clear;
nt=5000; ng=280; 
dx=10; dt=0.001; 
dz=dx; nz=170;
zg=0; ns=1;

load den; imagesc(den); title('Model');
figure;

load s;

kx=faxis(2*zg+ng,dx); f=faxis(nt,dt); 
[gkx,gf]=meshgrid(kx,f);
signgf=sign(gf);
gf2=gf.^2;
gkx2=gkx.^2;

Z=zeros(nt,zg);
for is=650:10:740
    is
    fn=sprintf('d%d.mat',is); 
    load(fn); data=d; 
    fn=sprintf('t%d.mat',is);
    load(fn);  
    t=round(t/dt); t(t>nt)=nt; t(t<1)=1;
    
    mig=zeros(nz,ng); 
    wu=[Z,data,Z];

    for iz=2:nz
        wu(:,1:zg)=Z;
        wu(:,ng+1+zg:end)=Z;
        fk=fft2(wu);
% PSPI        
        mins=min(s(iz,:));
        maxs=max(s(iz,:));
        weight=s(iz,:);
        for ix=1:ng
           weight(ix)=(1/s(ix)-1/maxs)/(1/mins-1/maxs); 
        end
        
        [kz,ind]=eikon(signgf,gf2,gkx2,mins^2,1);
        f1=fk.*ind.*exp(i*kz*dz);
        f1=real(ifft2(f1));
        
        [kz,ind]=eikon(signgf,gf2,gkx2,maxs^2,1);
        f2=fk.*ind.*exp(i*kz*dz);
        f2=real(ifft2(f2));
        
        for ix=zg+1:zg+ng
           wu(:,ix)=weight(ix-zg)*f1(:,ix)+(1-weight(ix-zg))*f2(:,ix); 
        end
        
        tt=t(iz,:);
        ind=sub2ind([nt,ng+2*zg],tt,[1:ng]+zg);
        mig(iz,:)=wu(ind);
        imagesc(wu); hold on; plot([1:ng]+zg,tt); hold off; pause(0.1)
    end
    imagesc(fliplr(mig'))
    fn=sprintf('mig%d.mat',is);
    save(fn,'mig');
end

