How can I create a sismic trace of the third component of the velocity vector of the solid matrix of the porous media in matlab using inverse fourier transform?

Question

Marcia Azeredo 2013 年 3 月 14 日

0
リンク

この質問への直接リンク

https://jp.mathworks.com/matlabcentral/answers/67210-how-can-i-create-a-sismic-trace-of-the-third-component-of-the-velocity-vector-of-the-solid-matrix-of

Hi. I have to plot in matlab the trace of the third component of the velocity vector of the solid matrix of a porous media. In my case, I have two homogeneous mediuns and one discontinuos interface between this mediuns. The problem is that I work in frequency domain, and to plot the trace, I have to transform my solution to time domain. For this, I have to use the inverse fourier transform, but don't use the ifft function. In my code, I use two functions that are fsqrtcomplexo to calculate the square of one complex number and iht to calculate the hankel inverse transform. The problem is that I have to construct the trace in real time, showing the reflections registred for the receptors. The code and the functions that I am using are:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% SIMULAÇÃO PARA UM MODELO GEOLÓGICO COM 2 CAMADAS %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Estudo do caso mais simples: fonte na superfície e receptor na %%% superfície (somente uma interface)

%%%%%%%% MODELO GEOLÓGICO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %______________________________________________ z0 (0 m) %---------------------------------------------- zs(local da fonte) (zs=0) % CAMADA 1 % % % %_____________________________________________ z1 (150 m) % CAMADA 2 (infinita) % % %---------------------------------------------(300 m)(limite de propagação)

clear all close all clc format long %alterando a precisão dos números exibidos pelo Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% Dados de entrada %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Propriedades das Camadas %%%%%%%%%%%%

% Camada 1: Rho1=2.4e3;Rhof1= 9.84e2; RhoE1= 0.0; Kappa1= 1.0e-16; Eta1=1.0e-3; Lambda1=2.69e9; G1=3.46e9; C1=14.4e9; M1=12.4e9;

% Camada 2: Rho2=2.4e3;Rhof2= 9.84e2; RhoE2= 0.0; Kappa2= 1.0e-13; Eta2=1.0e-3; Lambda2=6.05e9; G2=7.78e9; C2=13.1e9; M2=12.7e9;

% frequencia dominante dada no artigo wd=15.0;

%frequencia de corte wc=3*sqrt(pi)*wd;

%%%%%%% Valores de beta%%%%%%%%%%%% Beta11=1/(C1^2-M1*(Lambda1+2*G1));%beta1 da camada 1 Beta12=1/(C2^2-M2*(Lambda2+2*G2));%beta 1 da camada 2

%%%%%%% Vetores de entrada %%%%%%%%%% KK1=linspace(0.1,0.5,5);%vetor com os valores de K1 (5 valores) KK2=linspace(0.1,0.5,5);%vetor com os valores de K2 (5 valores) raio = linspace(.0,1000,5);%raio=(0,250,500,750,1000)

freq=linspace(1,150,150);%vetor de frequências que irei trabalhar N=length(freq); t=linspace(0,1,N); %vetor de tempo

%%%%%%% Posição da interface, da fonte e do local da solução %%%%%%%%%%%%%% z=0; %analisando a solucao apenas na superficie zs=0; %analisando com a fonte apenas na superficie z1=500;%posição da interface de discontinuidade

%pulso de ricker no dominio do tempo for nt=1:length(t) H(nt)=(1-2*pi^(2)*wd^(2)*t(nt)^(2))*exp(1)^(-pi^(2)*wd^(2)*t(nt)^(2)); end

%plot(t,H) %transformada de fourier do pulso de ricker nfft=N; hhh=zeros(1,nfft); Sum=0; for g=1:nfft for jj=1:N Sum=Sum+H(jj)*exp(-2*pi*j*(jj-1)*(g-1)/nfft); end hhh(g)=Sum; Sum=0;% Reset end %plot(freq,hhh) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% INÍCIO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 a=0;%contador para montar o vetor que guarda os valores de K
 for nk1=1:length(KK1) %contador para número de termos de K1
    for nk2=1:length(KK2) %contador para número de termos de K2       
       a=a+1;      
        for nw=1:length(freq)%contador para número de termos de w
        w=freq(nw);%escolhendo w da posição nw do vetor ww
        K1 = KK1(nk1);%escolhendo K1 da posição nk1 do vetor KK1
        K2 = KK2(nk2);%escolhendo K2 da posição nk2 do vetor KK2
        K=sqrt(K1^(2)+K2^(2));
        KK(a)=K; %vetor KK é o vetor que vai guardando os valores de K 
        Gama=K/w;
       h=hhh(nw);
        %h=(2.0/sqrt(pi))*(w^(2)/wd^(3))*exp(-w^(2)/wd^(2));%espectro do momento sísmico da fonte de Ricker 
       hh(nw)=h;%vetor que guarda os valores de h
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %                        MATRIZES M1 e M2                                 %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %Ordem de leitura: Mijk - matriz i  do sistema j da camada k
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            M111=[-Beta11*M1, Beta11*Gama*(C1^2-Lambda1*M1), -Beta11*C1;Beta11*Gama*(C1^2-Lambda1*M1),Rho1+((1i*w*Rhof1^(2)*Kappa1)/(Eta1-1i*w*Kappa1*RhoE1))-4*Beta11*Gama^(2)*G1*(C1^(2)-M1*(Lambda1+G1)),2*Beta11*Gama*G1*C1-(1i*w*Rhof1*Gama*Kappa1)/(Eta1-1i*w*Kappa1*RhoE1);-Beta11*C1,2*Beta11*Gama*G1*C1-(1i*w*Rhof1*Gama*Kappa1)/(Eta1-1i*w*Kappa1*RhoE1),-Beta11*(Lambda1+2*G1)+(1i*w*Gama^(2)*Kappa1)/(Eta1-1i*w*Kappa1*RhoE1)];
            M211=[Rho1, Gama, -Rhof1;Gama, 1/G1, 0;-Rhof1, 0, -1*(Eta1-1i*w*Kappa1*RhoE1)/(1i*w*Kappa1)];
            %sistema 2
            M121=1/G1;
            M221=Rho1-G1*Gama^2+(1i*w*Rhof1^2*Kappa1)/(Eta1-1i*w*RhoE1*Kappa1);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            M112=[-Beta12*M2, Beta12*Gama*(C2^2-Lambda2*M2), -Beta12*C2;Beta12*Gama*(C2^2-Lambda2*M2),Rho2+(1i*w*Rhof2^(2)*Kappa2)/(Eta2-1i*w*Kappa2*RhoE2)-4*Beta12*Gama^(2)*G2*(C2^(2)-M2*(Lambda2+G2)),2*Beta12*Gama*G2*C2-(1i*w*Rhof2*Gama*Kappa2)/(Eta2-1i*w*Kappa2*RhoE2);-Beta12*C2,2*Beta12*Gama*G2*C2-(1i*w*Rhof2*Gama*Kappa2)/(Eta2-1i*w*Kappa2*RhoE2),-Beta12*(Lambda2+2*G2)+(1i*w*Gama^(2)*Kappa2)/(Eta2-1i*w*Kappa2*RhoE2)];
            M212=[Rho2, Gama, -Rhof2;Gama, 1/G2, 0;-Rhof2, 0, (-Eta2+1i*w*Kappa2*RhoE2)/(1i*w*Kappa2)];
            %sistema 2
            M122=1/G2;
            M222=Rho2-G2*Gama^2+(1i*w*Rhof2^2*Kappa2)/(Eta2-1i*w*RhoE2*Kappa2);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %                 AUTOVALORES E AUTOVETORES                               %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            %autovalores
            qquad1=fsqrtcomplexo((1i*(Lambda1+2*G1)*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))-M1*Rho1)^2-4*(M1*Rhof1-1i*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))*C1)*(C1*Rho1-(Lambda1+2*G1)*Rhof1));
            qquadrado111=-Gama^2+Beta11*(C1*Rhof1-(1/2)*M1*Rho1-0.5i*(Lambda1+2*G1)*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))+(0.5)*qquad1);
            qquadrado211=-Gama^2+Beta11*(C1*Rhof1-(1/2)*M1*Rho1-0.5i*(Lambda1+2*G1)*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))-0.5*qquad1);
            qquadrado311=-1*Gama^2+1/G1*(Rho1+1i*Rhof1^2*w*Kappa1/(Eta1-1i*w*RhoE1*Kappa1));
            D11=diag([fsqrtcomplexo(qquadrado111) fsqrtcomplexo(qquadrado211) fsqrtcomplexo(qquadrado311)]);
            %autovetores
            ksi111=(C1*Rho1-(Lambda1 +2*G1)*Rhof1)/((qquadrado111+Gama^2)/(Beta11)-C1*Rhof1+1i*(Lambda1+2*G1)*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w)));
            ksi211=(C1*Rho1-(Lambda1 +2*G1)*Rhof1)/((qquadrado211+Gama^2)/(Beta11)-C1*Rhof1+1i*(Lambda1+2*G1)*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w)));
            ksi311=0.0;
            abarra111=fsqrtcomplexo((D11(1,1))/(Rho1+2*Rhof1*ksi111+1i*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))*ksi111^(2)));
            abarra211=fsqrtcomplexo((D11(2,2))/(Rho1+2*Rhof1*ksi211+1i*((Eta1-1i*w*RhoE1*Kappa1)/(Kappa1*w))*ksi211^(2)));
            bbarra311=fsqrtcomplexo((D11(3,3))/(G1*(qquadrado311+Gama^2)));
            %matriz L111
            L111(1,1)=-abarra111;
            L111(2,1)=abarra111*2*Gama*G1;
            L111(3,1)=abarra111*ksi111;
            L111(1,2)=-abarra211;
            L111(2,2)=abarra211*2*Gama*G1;
            L111(3,2)=abarra211*ksi211;
            L111(1,3)=Gama*bbarra311/D11(3,3);
            L111(2,3)=-bbarra311/D11(3,3)*(G1*(Gama^2-qquadrado311));
            L111(3,3)=-bbarra311/D11(3,3)*1i*Gama*Rhof1*w*Kappa1/(Eta1-1i*w*RhoE1*Kappa1);
            %matriz L211
            L211(1,1)=abarra111/D11(1,1)*(2*Gama^2*G1-Rho1-Rhof1*ksi111);
            L211(2,1)=abarra111/D11(1,1)*Gama;
            L211(3,1)=abarra111/D11(1,1)*(Rhof1+1i*((Eta1-1i*w*RhoE1*Kappa1)/(w*Kappa1))*ksi111);
            L211(1,2)=abarra211/D11(2,2)*(2*Gama^2*G1-Rho1-Rhof1*ksi211);
            L211(2,2)=abarra211/D11(2,2)*Gama;
            L211(3,2)=abarra211/D11(2,2)*(Rhof1+1i*((Eta1-1i*w*RhoE1*Kappa1)/(w*Kappa1))*ksi211);
            L211(1,3)=bbarra311*2*Gama*G1;
            L211(2,3)=bbarra311;
            L211(3,3)=0.0;
            %sistema 2
            %autovalores
            qquadrado121=-1*Gama^2+1/G1*(Rho1+1i*Rhof1^2*w*Kappa1/(Eta1-1i*w*RhoE1*Kappa1));
            D12=diag(fsqrtcomplexo(qquadrado121));
            %autovetores
            L121(1,1)=1/fsqrtcomplexo(G1*fsqrtcomplexo(qquadrado121));
            L121=L121;
            L221(1,1)=fsqrtcomplexo(G1*fsqrtcomplexo(qquadrado121));
            L221=L221;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            %autovalores
            qquadrado112=-Gama^2+Beta12*(C2*Rhof2-0.5*M2*Rho2-0.5i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))+0.5*fsqrtcomplexo((1i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))-M2*Rho2)^2-4*(M2*Rhof2-1i*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))*C2)*(C2*Rho2-(Lambda2+2*G2)*Rhof2)));
            q112=fsqrtcomplexo(qquadrado112);
            qquadrado212=-Gama^(2)+Beta12*(C2*Rhof2-0.5*M2*Rho2-0.5i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))-0.5*fsqrtcomplexo((1i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))-M2*Rho2)^2-4*(M2*Rhof2-1i*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))*C2)*(C2*Rho2-(Lambda2+2*G2)*Rhof2)));
            q212=fsqrtcomplexo(qquadrado212);
            qquadrado312=-1*Gama^2+1/G2*(Rho2+(1i*Rhof2^2*w*Kappa2)/(Eta2-1i*w*RhoE2*Kappa2));
            q312=fsqrtcomplexo(qquadrado312);
            D21=diag([q112,q212,q312]);%Matriz diagonal com os autovalores da camada 2(Lambda maiusculo)
            %autovetores
            ksi112=(C2*Rho2-(Lambda2+2*G2)*Rhof2)/((qquadrado112+Gama^(2))/(Beta12)-C2*Rhof2+1i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2)));
            ksi212=(C2*Rho2-(Lambda2+2*G2)*Rhof2)/((qquadrado212+Gama^(2))/(Beta12)-C2*Rhof2+1i*(Lambda2+2*G2)*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2)));
            abarra112=fsqrtcomplexo(q112/(Rho2+2*Rhof2*ksi112+1i*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))*ksi112^2));
            abarra212=fsqrtcomplexo(q212/(Rho2+2*Rhof2*ksi212+1i*((Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2))*ksi212^2));
            bbarra312=fsqrtcomplexo(q312/(G2*(qquadrado312+Gama^2)));
            %matriz L112
            L112(1,1)=-abarra112;
            L112(2,1)=abarra112*2.0*Gama*G2;
            L112(3,1)=abarra112*ksi112;
            L112(1,2)=-abarra212;
            L112(2,2)=abarra212*2*Gama*G2;
            L112(3,2)=abarra212*ksi212;
            L112(1,3)=bbarra312/q312*Gama;
            L112(2,3)=-bbarra312/q312*G2*(Gama^2-qquadrado312);
            L112(3,3)=-bbarra312/q312*1i*Gama*Rhof2*w*Kappa2/(Eta2-1i*w*RhoE2*Kappa2);
            %matriz L212
            L212(1,1)=abarra112/q112*(2*Gama^2*G2-Rho2-Rhof2*ksi112);
            L212(2,1)=abarra112/q112*Gama;
            L212(3,1)=abarra112/q112*(Rhof2+1i*(Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2)*ksi112);
            L212(1,2)=abarra212/q212*(2*Gama^2*G2-Rho2-Rhof2*ksi212);
            L212(2,2)=abarra212/q212*Gama;
            L212(3,2)=abarra212/q212*(Rhof2+1i*(Eta2-1i*w*RhoE2*Kappa2)/(w*Kappa2)*ksi212);
            L212(1,3)=bbarra312*2.0*Gama*G2;
            L212(2,3)=bbarra312;
            L212(3,3)=0.0;%bbarra312*0.0;
            %sistema 2
            %autovalores
            qquadrado122=-Gama^2+1/G2*(Rho2+1i*Rhof2^2*((w*Kappa2)/(Eta2-1i*w*RhoE2*Kappa2)));
            D22=diag(fsqrtcomplexo(qquadrado122));%Matriz diagonal com os autovalores da camada 2(Lambda maiusculo)
            %autovetores
            L122(1,1)=1/fsqrtcomplexo(G2*fsqrtcomplexo(qquadrado122));
            L222(1,1)=fsqrtcomplexo(G2*fsqrtcomplexo(qquadrado122));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %                        MATRIZES DE SALTO J                              %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            JA11z=(1/2)*(L211.'*L111+L111.'*L211);%matriz JA do sistema 1 na camada 1 (interface ficticia z, onde z está na camada 1)
            %JA11z=eye(3);
            JB11z=(1/2)*(L211.'*L111-L111.'*L211);%matriz JB do sistema 1 na camada 1 (interface ficticia z, onde z está na camada 1)
            %JB11z=zeros(3);
            %sistema 2
            JA21z=0.5*(L221.'*L121+L121.'*L221);%matriz JA do sistema 2 na camada 1 (interface ficticia z,onde z está na camada 1)
            JB21z=0.5*(L221.'*L121-L121.'*L221);%matriz JB do sistema 2 na camada 1 (interface ficticia z, onde z está na camada 1)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 2 (interface z1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            JA12=0.5*(L212.'*L111+L112.'*L211);%Matriz JA do sistema 1 da camada 2 (interface z1)
            JB12=0.5*(L212.'*L111-L112.'*L211);%Matriz JB do sistema 1 da camada 2 (interface z1)
            %sistema 2
            JA22=0.5*(L222.'*L121+L122.'*L221);%Matriz JA do sistema 2 da camada 2 (interface z1)
            JB22=0.5*(L222.'*L121-L122.'*L221);%Matriz JB do sistema 2 da camada 2 (interface z1)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %                  MATRIZES DE TRANSMISSÃO E REFLEXÃO                     %
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 2 (interface z1-última interface) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            R12=(-JB12.')/(JA12.');%matriz de reflexão do sistema 1 da interfaze z1
            R12tils=(exp(1)^(1i*w*(z1-zs)*D11))*R12*(exp(1)^(1i*w*(z1-zs)*D11));%será usado para computar a matriz de reflexão na interface fictícia zs+, qdo zs estiver na camada 1
            %R12tils=[exp(1)^(i*w*(z1-zs)*D11(1,1)),0,0;0,exp(1)^(i*w*(z1-zs)*D11(2,2)),0;0,0,exp(1)^(i*w*(z1-zs)*D11(3,3))]*R12*[exp(1)^(i*w*(z1-zs)*D11(1,1)),0,0;0,exp(1)^(i*w*(z1-zs)*D11(2,2)),0;0,0,exp(1)^(i*w*(z1-zs)*D11(3,3))];
            %sistema 2
            R22=(-JB22.')/(JA22.');%matriz de reflexão do sistema 2 da interfaze z1
            R22tils=(exp(1)^(1i*w*(z1-zs)*D12))*R22*(exp(1)^(1i*w*(z1-zs)*D12));%será usado para computar a matriz de reflexão na interface fictícia zs+, qdo zs estiver na camada 1
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%CAMADA 1 (interface ficticia zs+)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            Rs1=(JA11z.'*R12tils-JB11z.')/(-1.0*JB11z.'*R12tils+JA11z.');%matriz de reflexão do sistema 1 na interface fictícia zs+, quando zs+ está na camada 1
            %sistema 2
            Rs2=(JA21z.'*R22tils-JB21z.')/(-1.0*JB21z.'*R22tils+JA21z.');%matriz de reflexão do sistema 2 na interface fictícia zs+, quando zs+ está na camada 1
            %%%%%%%%%%%%%%%%%%VETORES DE FONTE (SA e SB) %%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            SA1=i*Beta11*h*[w*(C1-M1);2*K*G1*(M1-C1);w*(Lambda1+2*G1-C1)];%SA do sistema 1 
            SB1=[0.0;0.0;0.0];%SB do sistema 1 
            %sistema 2
            SA2=0.0;%SA do sistema 2 
            SB2=0.0;%SB do sistema 2 
            %%%%%%%%%%%%%%MATRIZES GA e GB  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %sistema 1
            GA1=[1.0,0.0,0.0;0.0,0.0,0.0;0.0,1.0,0.0];%matriz GA do sistema 1
            GB1=[0.0,0.0,0.0;0.0,0.0,1.0;0.0,0.0,0.0];%matriz GB do sistema 1
            %sistema 2
            GA2=1.0;%matriz GA do sistema 2
            GB2=0.0;%matriz GB do sistema 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%% FONTE (zs) NA CAMADA 1 %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Phig1=inv((exp(1)^(i*w*D11*zs))*Rs1*(exp(1)^(i*w*D11*zs))*(L211.'*GA1-L111.'*GB1)-(L211.'*GA1+L111.'*GB1))*(exp(1)^(i*w*D11*zs))*(Rs1*(L211.'*SA1-L111.'*SB1)-(L211.'*SA1+L111.'*SB1));%soluçao Phig do sistema 1 camada 1 Phig2=inv((exp(1)^(i*w*D12*zs))*Rs2*(exp(1)^(i*w*D12*zs))*(L221.'*GA2-L121.'*GB2)-(L221.'*GA2+L121.'*GB2))*(exp(1)^(i*w*D12*zs))*(Rs2*(L221.'*SA2-L121.'*SB2)-(L221.'*SA2+L121.'*SB2));%soluçao Phig do sistema 2 camada 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%% SOLUÇÃO NA SUPERFÍCIE %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Phi1=[GA1*Phig1;GB1*Phig1];  %solução Phi para o sistema 1 em z=0+
 Phi2=[GA2*Phig2;GB2*Phig2]; %solução Phi para o sistema 2 em z=0+

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% Terceira componente da velocidade de fase sólida %%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

utilponto3(nw,a)=Phi1(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Rotação inversa na terceira componente da velocidade de fase sólida %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%saindo de ~u. para voltar para û. uchapeuponto3(nw,a) =utilponto3(nw,a);%matriz que guarda os valores de uchapeuponto3(=utilponto3) para cada frequencia e K

        end
end
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% Transformada de Hankel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %saindo de û. para u. (mudando da frequencia espacial para a distância %radial for nw=1:length(freq); uponto3(nw,:)=iht(uchapeuponto3(nw,:),KK,raio); end uponto3raio1=uponto3(:,1);%vetor que guarda os valores de uponto3 para raio igual a 1 e todas as frequencias uponto3raio2=uponto3(:,2); uponto3raio3=uponto3(:,3); uponto3raio4=uponto3(:,4); uponto3raio5=uponto3(:,5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%% Construindo o gráfico do pulso de Ricker %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Pulso de Ricker no domínio da frequência %%%%%%%%%% %plot(freq,real(hh)) %xlabel('frequência (Hz)') %ylabel('amplitude (m/s)') %title('Espectro do pulso sísmico de Ricker com frequência dominante de 15 Hz')

%%%% Pulso de Ricker no domínio do tempo %%%%%%%%%% %hhtempo=ifft(hh); %hhtempoc=ifftshift(hhtempo)*freqamostra; %plot(t,real(hhtempoc)) %xlabel('tempo (s)') %ylabel('amplitude (m/s)') %title('Traço do pulso sísmico de Ricker com frequência dominante de 15 Hz')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% Transformada Inversa de Fourier no tempo %%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nfft=N; uponto3temporaio3=zeros(1,nfft); Sum=0; for pp=1:nfft for jj=1:N Sum=Sum+uponto3raio3(jj)*exp(-2*pi*j*(jj-1)*(pp-1)/nfft);%transf. inversa de fourier de uponto3 para raio 1 end uponto3temporaio3(pp)=Sum; Sum=0;% Reset end uponto3temporaio3=uponto3temporaio3/N; figure, plot(t, real(uponto3temporaio3),'LineWidth',1.9), title(' Inverse Fourier Transform F^{-1}[F[H(x)]]=H(x)'), grid on %axis([-10 30 -1 2])

%xlabel('tempo (s)') %ylabel('amplitude (m/s)') %title('Traço da terceira componente da velocidade de fase sólida com afastamento fonte-receptor igual a 0 m e frequência dominante de 15 Hz')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sqrtcomplexo=fsqrtcomplexo(h)

sqrtcomplexo=sqrt(h); if imag(sqrtcomplexo)==0; if real(sqrtcomplexo)<0; sqrtcomplexo=-1*real(sqrtcomplexo); end else if imag(sqrtcomplexo)<0; sqrtcomplexo=-1*(sqrtcomplexo); end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%[h,J]=iht(H,k,r,J) %------------------ % %Inverse Hankel transform of order 0. % %Input: % H Spectrum K(k) % k Spatial frequencies [rad/m] {pi/numel(H)*(0:numel(H)-1)} % r Radial positions [m] {0:numel(H)-1} % J Integration kernel {default} % %Output: % h Signal h(r) % J Integration kernel % % % ) If the integration kernel is missing, it is % recomputed from the Bessel functions (slow). %

% Marcel Leutenegger June 2006 % function [h,J]=iht(H,k,r,J) if sum(size(H) > 1) > 1 error('Spectrum must be a vector.'); end if nargin < 2 | isempty(k) k=pi/numel(H)*(0:numel(H)-1).'; else [k,w]=sort(k(:).'); H=H(w); end if nargin < 3 | isempty(r) r=0:numel(H)-1; end if nargin < 4 | isempty(J) k=[(k(2:end) + k(1:end-1))/2 k(end)]; J=1/2/pi./r(:)*k.*besselj(1,r(:)*k); J(r == 0,:)=1/4/pi*k.*k; J=J - [zeros(numel(r),1) J(:,1:end-1)]; elseif exist('w','var') J=J(:,w); end h=reshape(J*H(:),size(r));