preffd_0227am_modeling.m

one wave prestack migration code using ffd method

function preffd_0227am_modeling
%modified by qinzhen /2009-02-26
nx=101;nz=101;nt=501; dx=10;dz=8;dt=0.002;
midfreq=50; fmax=0.5/dt;
df=1/nt/dt; dw=2*pi*df;
nw=round(3.5*midfreq/df); fw=0.001*dw;
if(nw>floor(nt/2)) nw=floor(nt/2);end
dwfl=zeros(nt,nx);%source data
uwfl=zeros(nt,nx);%receive data
ssp=zeros(nz,nx);%slowness data
rrr=zeros(nz-1,nx);%reflectivity data
ssp(1:nz,1:nx)=3500; ssp(1:20,:)=1500;ssp(21:40,:)=2000;
ssp(41:60,:)=2500; ssp(61:80,:)=3000;ssp(41:60,1:50)=2000;
ssp(1:nz,1:nx)=1./ssp;
rrr(19,:)=0.3; rrr(39,:)=0.3; rrr(59,:)=0.3; rrr(79,:)=0.3;
rrr(59,1:50)=0.1;
rrr(41:60,60)=0.1;
www=zeros(nz,nx,nw);  %wavefront waveshot

nzero=40;
nx1=nx+2*nzero;
dkx=2*pi/nx1/dx;
hamfilter=ones(1,nx1);%Hamming at the same time in x domain
htemp=blackman(2*nzero); %real size is 81
hamfilter(1:nzero)=htemp(1:nzero);
hamfilter(nx1-nzero+1:nx1)=htemp(nzero:-1:1);
slowness=ones(1,nx1);
% hamfilter(nx1-19:nx1)=0;hamfilter(1:20)=0;

for it=1:nt   %seismic source,ricker wave
    temp=pi*midfreq*(it-ceil((nt+1)/2))*dt;
    sd0(it)=(1-2*temp*temp)*exp(-temp*temp);
end
dwfl(:,20)=fft(ifftshift(sd0));

kx=ifftshift(([1:nx1]-ceil((nx1+1)/2))*dkx);
kx2=kx.*kx;

for iw=1:nw
    w=fw+(iw-1)*dw;%+0.001*i*dw;% stable??
    ww_s=zeros(1,nx1);
    ww_s(nzero+1:nzero+nx)=dwfl(iw,1:nx);
    www(1,1:nx,1)=ww_s(nzero+1:nzero+nx);
    ww_s=ww_s.*hamfilter;
    for iz=2:nz  
        slowness(1+nzero:nx+nzero)=ssp(iz,1:nx);
        slowness(1:nzero)=ssp(iz,1);
        slowness(nx+nzero+1:nx1)=ssp(iz,nx); 
        ww_s=qz_ffd(ww_s,kx2,w,dz,slowness);
        ww_s=ww_s.*hamfilter;
        www(iz,1:nx,iw)=conj(ww_s(nzero+1:nzero+nx));
        dwfl(iw,1:nx)=www(iz,1:nx,iw);
    end
    
    ww_r=zeros(1,nx1);
    for iz=nz-1:-1:7
         ww_r(nzero+1:nzero+nx)=ww_r(nzero+1:nzero+nx)+conj(www(iz+1,1:nx,iw)).*rrr(iz,1:nx);
%         ww_r=qz_ffd(ww_r,1,kx2,w,-dz,dx,dt,hamfilter,nx,nzero);
        slowness(1+nzero:nx+nzero)=ssp(iz,1:nx);
        slowness(1:nzero)=ssp(iz,1);
        slowness(nx+nzero+1:nx1)=ssp(iz,nx);
        ww_r=qz_sp(ww_r,kx2,w,dz,slowness);
        ww_r=ww_r.*hamfilter;
    end
    uwfl(iw,1:nx)=conj(ww_r(nzero+1:nzero+nx));
end

temp6=zeros(1,nt);  %compute certain time waveshot
out=zeros(nz,nx);
ntmp=ceil(nt/2);
for ix=1:nx
    for iz=1:nz
        temp6(1:nw)=www(iz,ix,1:nw);
        temp6(nt:-1:nt-ntmp+2)=conj(temp6(2:ntmp));
        temp7=ifft(temp6);
        out(iz,ix)=temp7(121)+temp7(80)+temp7(40);
    end
end

for ix=1:nx
    temp=zeros(1,nt);
    temp(1:nw)=uwfl(1:nw,ix);
    temp(nt:-1:nt-nw+2)=conj(temp(2:nw));
    uwfl(:,ix)=real(ifft(temp));
end

figure(3);imagesc(real(out));title('several time waveshot');
figure(4);imagesc(real(www(:,:,60)));title('single frequency waveshot');
figure(5);imagesc(uwfl);title('gather record');
end

%************end main program********************************************
%************************************************************************
function ww_r=qz_p(ww_r,kx2,w,dz,slowness)  
        wk=fft(ww_r);
        slow0=mean(slowness);
        kt=slow0*w;
        kt2=kt*kt;
        kz=kt2-kx2;
        if(dz<0) kz(kz<0)=0;end
        wk=wk.*exp(i*dz*sqrt(kz));
        ww_r=ifft(wk);       
end

function ww_r=qz_sp(ww_r,kx2,w,dz,slowness)  
        wk=fft(ww_r);
        slow0=max(slowness);
        kt=slow0*w;
        kt2=kt*kt;
        kz=kt2-kx2;
        if(dz<0) kz(kz<0)=0;end
        wk=wk.*exp(i*dz*sqrt(kz));
        ww_r=ifft(wk);    
        ww_r=ww_r.*exp(i*w*dz*(slowness-slow0));  %frequency shift   
      %  ww_r=fdmig(ww_r , 161, w, 2000,10,10,65,1);
end

function ww_r=qz_fd(ww_r,kx2,w,dz,slowness)  
%         wk=fft(ww_r);
%         slow0=1/1000;
%         kt=slow0*w;
%         kt2=kt*kt;
%         kz=sqrt(kt2-kx2);
%         wk=wk.*exp(i*dz*kz);
%         ww_r=ifft(wk);  
        ww_r=ww_r.*exp(-i*w*dz/2000);  %frequency shift
        ww_r=fdmig(ww_r , nx1, w, 1./sv,dz,10,45,1);  
        ww_r=ww_r.*hamfilter;
end


function ww_r=qz_ffd(ww_r,kx2,w,dz,slowness)
        slow0=min(slowness);
        nx1=length(kx2);
        %start w-x
        ww_s=fdmig( ww_r, nx1, w, 1./slowness,dz,10.0,65);
%         ww_r=fdmig(ww_r , nx1, w, ssp,dz,dx,65,-1);
        %end w-x
        ww_r=ww_r.*exp(i*w*dz*slow0);  %frequency shift
        %start w-k
        wk=fft(ww_r);
        slow0=min(slowness);
        kt=slow0*w;
        kt2=kt*kt;
        kz=kt2-kx2;
        if(dz<0) kz(kz<0)=0;end
        wk=wk.*exp(i*dz*sqrt(kz));
        ww_r=ifft(wk);    
        %end w-k        
end


%complex tridiagonal equation solver
function q=ctris(n,endl,a,c,b,endr,d)
e=zeros(1,n); f=zeros(1,n);
e(1)=-a(1)/endl;
f(1)=d(1)/endl;
for i=2:n-1
    den=b(i)+c(i)*e(i-1);
    e(i)=-a(i)/den;
    f(i)=(d(i)-c(i)*f(i-1))/den;
end
q(n)=(d(n)-c(n-1)*f(n-1))/(endr+c(n-1)*e(n-1));
for i=n-1:-1:1
    q(i)=e(i)*q(i+1)+f(i);
end
end

function [data]=fdmig( cp, nx, w, v,dz,dx,dip)
	trick=0.1;
	p=zeros(1,nx);
    s1=zeros(1,nx);
    s2=zeros(1,nx);
	data=zeros(1,nx);
	d=zeros(1,nx);
	a=zeros(1,nx);
	b=zeros(1,nx);
	c=zeros(1,nx);

	vc=min(v);
    for ix=1:nx
		p(ix)=vc/v(ix);
		p(ix)=(p(ix)*p(ix)+p(ix)+1.0);
    end

	if(dip~=65)
		coefa=0.5;coefb=0.25;
    else
		coefa=0.4784689;
		coefb=0.37607656;
    end
    v1=v(1);vn=v(nx);
    
	if(abs(w)<=1.0e-10) w=1.0e-10/dt; end
    s1(1:nx)=(v.*v).*p*coefb/(dx*dx*w*w)+trick;
    s2(1:nx)=-(1-vc./v).*v*dz*coefa/(w*dx*dx)*0.5;

    data=cp;
%------ready a1,a2 for boundary condition
    %conpute for a1,b1
	cp2=data(2);    cp3=data(3);
    a1=    2.0*(cp2*conj(cp3));
    b1= cp2*conj(cp2)+cp3*conj(cp3);
    if(b1==0.0+i*0)
		a1=exp(-i*w*dx*0.5/v1);
	else
		a1=a1/b1;
    end
	if(imag(a1)>0.0) a1=exp(-i*w*dx*0.5/v(1)); end
    %conpute for a2,b2
    cpnm1=data(nx-1);  cpnm2=data(nx-2);
	a2=2.0*cpnm1*conj(cpnm2);  
    b2=cpnm1*conj(cpnm1)+cpnm2*conj(cpnm2);	
    if(b2==0.0+i*0)
		a2=exp(-i*w*dx*0.5/vn);
	else
		a2=a2/b2;
    end
    if(imag(a2)>0.0) a2=exp(-i*w*dx*0.5/v(nx)); end
%--------------------
%---solve for (I+a*Lx)U(z+dz)=(I-a*Lx)U(z)
%---  a=s1+i*s2; b=I-2*a
	a=s1+s2*i;
    b=1-2*a;
	for ix=2:nx-1
		d(ix)=data(ix+1)*a(ix+1)+data(ix-1)*a(ix-1)+data(ix)*b(ix);
    end
    d(1)=(b(1)+a1*a(1))*data(1)+a(2)*data(2);
	d(nx)=(b(nx)+a2*a(nx))*data(nx)+a(nx-1)*data(nx-1);
	data=s1-i*s2;
	b=1-2*data;
	endl=b(1)+a1*data(1);  
	endr=b(nx)+a2*data(nx);  
	for ix=2:nx-1
		a(ix)=data(ix+1);
		c(ix)=data(ix-1);
    end
	a(1)=data(2);
	c(nx)=data(nx-1);
        
	data=ctris(nx,endl,a,c,b,endr,d);	
  end