function preffd_0227am_modeling%modified by qinzhen /2009-02-26nx=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);enddwfl=zeros(nt,nx);%source datauwfl=zeros(nt,nx);%receive datassp=zeros(nz,nx);%slowness datarrr=zeros(nz-1,nx);%reflectivity datassp(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 waveshotnzero=40;nx1=nx+2*nzero;dkx=2*pi/nx1/dx;hamfilter=ones(1,nx1);%Hamming at the same time in x domainhtemp=blackman(2*nzero); %real size is 81hamfilter(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);enddwfl(:,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));endtemp6=zeros(1,nt);  %compute certain time waveshotout=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);    endendfor 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));endfigure(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);       endfunction 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);endfunction 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;endfunction 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 solverfunction 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;endq(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);endendfunction [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