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?
1 view (last 30 days)
Show older comments
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));
0 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)