ssf.m

分裂步程序matlab演示版

% Define
PI=3.14159165;
dt=0.004;
dx=30;
dz=14;

% Get Parameters:      
 %---------------
 nt=1024;
 nz=200;
 nx=257;

 %---------------
 ntfft =1040;
 nxfft =260;
 fw=0.0;
 fz=0;

 %---------------
 nw = ntfft/2+1;
 dw = 2.0*PI/(ntfft*dt);
 fw = 0;   % w ranges from 0 to Pi
 nk = nxfft;
 dk = 2.0*PI/(nxfft*dx);
 fk = -PI/dx;  % k rangs from -Pi to Pi, coz the Fourier Kernel is 1。

% clarification matrix
cresult=zeros(nx,nz);
CQ1=zeros(nk,nw);
 
% Read Data
 % Load traces into the zero-offset array
 % [D,H] = readsegy('spike.su');
 load data.mat

 % Load velicoty file
  % old mathod
   load v.mat;
%  FID = fopen('vfile');
%  status = fseek(FID,0,'bof');
%  v=fread(FID,[200,257],'float');
    
% Pad data with zeros and apply FFT
T=zeros(ntfft-nt,nx);
P=[D;T];
 CP=ntfft*ifft(P);
%----test-----      % not correct! Still cannot figure out '+1'&'-1'
%CP=fft(P);        % means what to the Fourier Kernel !!!
                     

CP=CP(1:nw,:);
clear T;
% clear D;
 
% loops over depth
profile on -detail 'builtin'
iz=1;
tz=fz;
while iz<=nz
   
    for ix=1:nx
        cresult(ix,iz)=0.0;
        for iw=1:nw
            cresult(ix,iz)=cresult(ix,iz)+real(CP(iw,ix))/ntfft;
        end
    end
    
    vmin=0;
    for ix=1:nx
        vmin=vmin+1.0/v(iz,ix)/nx;
    end
    vmin=1.0/vmin*0.5;
        
    % ???---------------------------------- 
     % data modified to fit SU-style FFT
     CQ=zeros(nw,nk);
     for ik=1:nx
         if mod(ik,2)==1
             CQ(:,ik)=CP(:,ik);
         else
             CQ(:,ik)=-CP(:,ik);
         end
     end 
    %   CQ(:,1:nx)=CP;
   % ???----------------------------------
        
    % FFT to W-K domain
    for iw=1:nw
        CQ(iw,:)=fft(CQ(iw,:));
     % CQ(iw,:)=fftshift(fft(CQ(iw,:)));  % No difference will be shown ?!
    end
     v1=vmin;
    
   %_________Phase Shift_________
    k=fk;
    ik=1;
    while ik<=nk        
       w=fw;
       iw=1;
       while iw<=nw
           %-------------------------------
           if w==0.0
               w=1.0e-10/dt;
           end
           kz1=1.0-(v1*k/w)^2;
           if kz1>0
               phase1=-w*sqrt(kz1)*dz/v1;
               cshift1=cos(phase1)+sqrt(-1)*sin(phase1);
               CQ1(ik,iw)=CQ(iw,ik)*cshift1;
           else 
               cshift1=0.0;
               CQ1(ik,iw)=CQ(iw,ik)*cshift1;
           end
           %--------------------------------
          iw=iw+1;
          w=w+dw;
       end                       
        ik=ik+1;
        k=k+dk;
    end
   %______________________________
    
   %-----FFT to W-X domain------
    CQ1=ifft(CQ1);
    
   %_________Time shift___________
   for ix=1:nx
      w=fw;
       for iw=1:nw           
           if mod(ix,2)==0
               CQ1(ix,iw)=-CQ1(ix,iw);
           else
               CQ1(ix,iw)=CQ1(ix,iw);
           end
           kz2=(1.0/v1-2.0/v(iz,ix))*w*dz;
           cshift2=cos(kz2)+sqrt(-1)*sin(kz2);
           CP(iw,ix)=CQ1(ix,iw)*cshift2;           
           w=w+dw;
       end           
   end
   %__________________________________
   
    tz=tz+dz;
    iz=iz+1
end

profile viewer
cresult=cresult';   

%---------------------------PLOT-------------------------------

x=(0:(nx-1))*dx;
z=linspace(0,(nt-1)*dt,nz);
[xx zz]=meshgrid(x,z);
clf;
pcolor(xx,zz,cresult);
set(gca, 'YDir', 'reverse','XAxisLocation', 'top')
shading interp
colormap winter