program OWmig
!
! x      => CMP coordinates(m)     vector
! xs     => shot coordinates(m)    vector
! h      => full offsets(m)        vector
! sh     => subsurf.off.(m)        vector
! z      => depth(m)               vector
! t      => time(ms)               vector
! v      => int.velocity(m/s)      matrix nz X nx
! source => source function        vector
! shots  => shot gathers           matrix nt X nh X nxs
! ixpad  => traces to pad          scalar
!

use sep 
!use OW_io

 call sep_init()

!
! IO
!

 call from_history("n1",rec%xmr%n);	call from_history("d1",rec%xmr%d);	call from_history("o1",rec%xmr%o)
 call from_history("n2",rec%ymr%n);	call from_history("d2",rec%ymr%d);	call from_history("o2",rec%ymr%o)  
 call from_history("n3",rec%wr%n);	call from_history("d3",rec%wr%d);	call from_history("o3",rec%wr%o)
 call from_history("n4",rec%xr%n);	call from_history("d4",rec%xr%d);	call from_history("o4",rec%xr%o)  
 call from_history("n5",rec%yr%n);	call from_history("d5",rec%yr%d);	call from_history("o5",rec%yr%o)

 call from_aux("sou","n1",wav%xms%n);	call from_aux("sou","d1",wav%xms%d);	call from_aux("sou","o1",wav%xms%o)
 call from_aux("sou","n2",wav%yms%n);	call from_aux("sou","d2",wav%yms%d);	call from_aux("sou","o2",wav%yms%o)  
 call from_aux("sou","n3",wav%ws%n);	call from_aux("sou","d3",wav%ws%d);	call from_aux("sou","o3",wav%ws%o)
 call from_aux("sou","n4",wav%xs%n);	call from_aux("sou","d4",wav%xs%d);	call from_aux("sou","o4",wav%xs%o)  
 call from_aux("sou","n5",wav%ys%n);	call from_aux("sou","d5",wav%ys%d);	call from_aux("sou","o5",wav%ys%o)

 call from_aux("vel","n1",vel%xm%n);	call from_aux("vel","d1",vel%xm%d);	call from_aux("vel","o1",vel%xm%o)
 call from_aux("vel","n2",vel%ym%n);	call from_aux("vel","d2",vel%ym%d);	call from_aux("vel","o2",vel%ym%o)  
 call from_aux("vel","n3",vel%z%n);	call from_aux("vel","d3",vel%z%d);	call from_aux("vel","o3",vel%z%o)

 if (wav%ws%n != rec%wr%n) call seperr("# source freq. != # receiver freq.")
 if (wav%ws%d != rec%wr%d) call seperr("source freq.interval != receiver freq.interval")
 if (wav%ws%o != rec%wr%o) call seperr("source freq.origin != receiver freq.origin")

 !!  get image space 
  call from_param("xmin",xmin) 
  call from_param("xmax",xmax)
  call from_param("ymin",ymin,0)
  call from_param("ymax",ymax,0)

 if (xmin < vel%xm%o ) call seperr("Min.x-image space < Min.x-velocity")
 if (xmax > vel%xm%o+(vel%xm%n-1)*vel%xm%d ) call seperr("Max.x-image space > Max.x-velocity")
 if (ymin < vel%ym%o ) call seperr("Min.y-image space < Min.y-velocity")
 if (ymax > vel%ym%o+(vel%ym%n-1)*vel%ym%d ) call seperr("Max.y-image space > Max.y-velocity")

 if (rec%xmr%o < xmin ) call seperr("Min.x-receiver < Min.x-image space")
 if (rec%xmr%o+(rec%xmr%n-1)*rec%xmr%d < xmax ) call seperr("Max.x-receiver < Max.x-image space")
 if (rec%ymr%o < ymin ) call seperr("Min.y-receiver < Min.y-image space")
 if (rec%ymr%o+(rec%ymr%n-1)*rec%ymr%d < ymax ) call seperr("Max.y-receiver < Max.y-image space")

 if (sou%xms%o < xmin ) call seperr("Min.x-source < Min.x-image space")
 if (sou%xms%o+(sou%xms%n-1)*sou%xms%d < xmax ) call seperr("Max.x-source < Max.x-image space")
 if (sou%yms%o < ymin ) call seperr("Min.y-source < Min.y-image space")
 if (sou%yms%o+(sou%yms%n-1)*sou%yms%d < ymax ) call seperr("Max.y-source < Max.y-image space")

 if (rec%xmr%d != vel%xm%d) call seperr("receiver x-interval != velocity x-interval")
 if (rec%ymr%d != vel%ym%d) call seperr("receiver y-interval != velocity y-interval")


  call from_param("emt4",emt(4))
  call from_param("emt5",emt(5))
  call from_param("emt6",emt(6))
  call from_param("emt7",emt(7))

sign=1;

isp = ixpad - floor(oh/dh) + 1;                                                                     ! shot position in the acq.geom.
!
! Frequency vector
!
nw = floor(nt/2);
dw = 1./(nt*dt);
! ! no zero frequency component
!ow = dw;        
ow=0;   
wall=zeros(nw,1);
for iw = 1:nw;
	wall(iw) = ow  + (iw - 1) * dw;
end;
wall = 2*pi * wall;
!
! Wavenumber vector
!
nkx = nh+2*ixpad;
dkx = 1./(nkx*dx);
okx = 1./(2*dx)-dkx;
kxall=zeros(nkx,1);
for ik=1:nkx;
	kxall(ik) = rem(okx+(ik-1)*dkx,1./dx)-okx;
end;
kxall = 2*pi * kxall;
!
! expand velocity if needed, assuming:
! icmp_vel = icmp_image
! iz_vel   = iz_image
!
[nzv,nxv]=size(v);
if nxv<nx;!+2*ixpad;
    vel=v;
    v=zeros(nzv,nx+2*ixpad);
    for iv=1:ixpad;
        v(:,iv)=vel(:,1);
        v(:,ixpad+nxv+iv)=vel(:,length(vel(1,:)));
    end;
    v(:,ixpad+1:ixpad+nxv)=vel;
end;
if nzv<nz;
    disp('entrei')
    vel=v;
    v=zeros(nz,length(v(1,:)));
    for iv=nzv+1:nz;
        v(iv,:)=vel(nzv,:);
    end;
    v(1:nzv,:)=vel;
end;
clear vel;
!
! Pad and position wavelet 
!
wavelet=zeros(nt,1);
if (rem(length(source),2) == 0);
    nsour2=floor(length(source)/2);
    wavelet(nw-nsour2+1:nw+nsour2)=source;
else;
	nsour2=floor(length(source)/2)+1;
    wavelet(nw-nsour2+1:nw+nsour2-1)=source;
end;
wavelet=half(-1,wavelet,dt);
wavelet=fft(wavelet);
wave=zeros(nw,1);wave=wavelet(1:nw);
wave=wave.*exp(complex(0.,wall*nw*dt));
!
!
!
vref=zeros(nref,nz,nsh);
for ixs=1:nxs;                                                                                          ! shot loop    
!
! Pre-compute vref
!
    owrk=(ixs-1)*dxs+oxs+oh-ixpad*dx;                                                                  ! workspace origin 
	if owrk < x(1) || owrk > x(end);
        error(['WORKING SPACE GREATER THAN IMAGE SPACE FOR SHOT #',num2str(ixs)]);
	end;
! grab velocity in the workspace       
    ni=floor((owrk-ox)/dx)+1;
	vel=zeros(nz,nkx);
    vel(:,:)=v(:,ni:ni+nkx-1);
    for iz=1:nz;                                                                                        ! depth loop
! !        minv=min(vel(iz,:));
! !        maxv=max(vel(iz,:));



! !        dvel=(maxv-minv)/(nref-1);
       vref(:,iz,ixs)=(minv:dvel:maxv);
        for iref=1:nref;
! !             vref(iref,iz,ixs)=minv+(iref-1)*dvel;
        k0=0;v0=0;
!        if iz >=120 && iz <=130;disp([num2str(iz),'  ',num2str(minv),'  ',num2str(maxv),'  ',num2str(min(v(iz,:))),'  ',num2str(max(v(iz,:)))]);end;
        if minv~=maxv;
            if maxv >= vsal;
                vref(nref,iz,ixs)=maxv;
                for ix=1:nkx;
                    if vel(iz,ix) < vsal;
                        v0=v0 + vel(iz,ix);
                        k0=k0+1;
                    end;
                end;
                v0=v0/k0;
                if nref == 2;
                    vref(1:nref-1,iz,ixs)=v0;
                else;
                    if v0 == minv;
                        vref(1:nref-1,iz,ixs)=v0;
                    else;
                        vref(1:nref-1,iz,ixs)=(minv:2*(v0-minv)/(nref-2):2*v0-minv);
                    end;
                end;
            else;
                vref(:,iz,ixs)=(minv:(maxv-minv)/(nref-1):maxv);
            end;
        else;
            vref(1:nref,iz,ixs)=minv;
!             vref(1:nref,iz,ixs)=mean(vel(iz,:));
!             if mean(vel(iz,:)) == 0;
!                 disp([])
! !         end;
    end;
end;
! figure;imagesc(permute(vref(1:3:nref,:,1),[2 1 3]));colormap('jet');pause;
!
!
!
disp('DOWNWARD CONTINUES SOURCES AND RECEIVERS')
disp('======== ========= ======= === =========')
!
image=zeros(nz,nsh,nx);sws=zeros(nw,nz,nkx);swr=sws;
for ixs=1:nxs;                                                                                          ! shot loop    
!
! FFT shots
!
    tempshots=zeros(nt,nh+2*ixpad);
    tempshots(:,ixpad+1:ixpad+nh)=shots(:,:,ixs);
    tempshots = fft(tempshots);
    GG=zeros(nz,nsh,nkx);
	owrk=(ixs-1)*dxs+oxs+oh-ixpad*dx;                                                                  ! workspace origin 
	if owrk < x(1) || owrk > x(end);
        error(['WORKING SPACE GREATER THAN IMAGE SPACE FOR SHOT #',num2str(ixs)]);
	end;
    
! grab velocity in the workspace       
    ni=floor((owrk-ox)/dx)+1;
	vel=zeros(nz,nkx);
    vel(:,:)=v(:,ni:ni+nkx-1);
!      
! exec.time
!
    cputa=cputime;
    cli1a=date; cli2a=clock;

	ws=zeros(nz,nkx);wr=ws;
    gs=zeros(nkx,nref);gs0=zeros(nkx,1);       
	gr=zeros(nkx,nref);gr0=zeros(nkx,1);               

    for iw = 1:nw;                                                                                      ! frequency loop
        ws(1,isp)=wave(iw);
        wr(1,:)=tempshots(iw,:);
        for iz=2:nz;                                                                                        ! depth loop
            ws(iz,:)=fft(ws(iz-1,:));
            wr(iz,:)=fft(wr(iz-1,:));
            minv=min(vel(iz,:)); maxv=max(vel(iz,:));

!      Operations in the wavenumber domain
           if minv == maxv;                                                                                ! v(x.z) condition
                gs0(:)=ws(iz,:);  gr0(:)=wr(iz,:);
                for ik = 1:nkx;                                                                         ! wavenumber loop
!      phase shift
                   sqa0 = wall(iw)^2 / minv^2 - kxall(ik)^2;
                   if sqa0 > 0;                                                                         ! evanescent condition v0
                       gs0(ik)=gs0(ik)*exp(complex(0,-sign*dz*sqrt(sqa0)));
                       gr0(ik)=gr0(ik)*exp(complex(0,sign*dz*sqrt(sqa0)));
                   else;
                       gs0(ik)=gs0(ik)*exp(-sqrt(-sqa0)*dz);
                       gr0(ik)=gr0(ik)*exp(-sqrt(-sqa0)*dz);
                   end;
                end;                                                                                    ! end wavenumber loop
! Operations in x            
                gs0=ifft(gs0);  gr0=ifft(gr0); 
!      split-step
                gs0 = gs0 .* exp(complex(0,-sign*dz*wall(iw)).*(1./vel(iz,:)-1./minv));;              ! spl.step
                gr0 = gr0 .* exp(complex(0,sign*dz*wall(iw)).*(1./vel(iz,:)-1./minv));;              ! spl.step
                ws(iz,:)=gs0;   wr(iz,:)=gr0;               
            else;                                                                                       ! v(z) condition
                for iref=1:nref;                                                                        ! reference velocities                 
                    gs(:,iref)=ws(iz,:);  gr(:,iref)=wr(iz,:); 
                    for ik = 1:nkx;                                                                         ! wavenumber loop
!      phase shift
                        sqa0 = wall(iw)^2 / vref(iref,iz,ixs)^2 - kxall(ik)^2;
                        if sqa0 > 0;                                                                        ! evanescent condition v0
                            gs(ik,iref)=gs(ik,iref)*exp(complex(0,-sign*dz*sqrt(sqa0)));
                            gr(ik,iref)=gr(ik,iref)*exp(complex(0,sign*dz*sqrt(sqa0)));
                        else;
                            gs(ik,iref)=0.;
                            gr(ik,iref)=0.;
                        end;
                    end;                                                                                    ! end wavenumber loop
! Operations in x            
                    gs(:,iref)=ifft(gs(:,iref)); 
                    gr(:,iref)=ifft(gr(:,iref)); 
!      split step
                    gs(:,iref) = gs(:,iref) .* exp(complex(0,-sign*dz*wall(iw)).*(1./vel(iz,:)-1./vref(iref,iz,ixs)));
                    gr(:,iref) = gr(:,iref) .* exp(complex(0,sign*dz*wall(iw)).*(1./vel(iz,:)-1./vref(iref,iz,ixs)));
                end;                                                                                       ! end reference velocities

                for ix=1:nkx;
                    for iref=1:nref-1;
                        if vref(iref,iz,ixs) ~= vref(iref+1,iz,ixs);
                            if vel(iz,ix) >= vref(iref,iz,ixs) && vel(iz,ix) < vref(iref+1,iz,ixs);
                                f=(vel(iz,ix) - vref(iref,iz,ixs))/vel(iz,ix); g=1-f;
                                ws(iz,ix)=g*gs(ix,iref)+f*gs(ix,iref+1);
                                wr(iz,ix)=g*gr(ix,iref)+f*gr(ix,iref+1);
!                                ws(iz,ix)=gs(ix,iref);
!                                wr(iz,ix)=gr(ix,iref);
                            elseif vel(iz,ix) >= vref(nref,iz,ixs);
                                ws(iz,ix)=gs(ix,nref);
                                wr(iz,ix)=gr(ix,nref);
                            end;
                        else;
                            ws(iz,ix)=gs(ix,iref);
                            wr(iz,ix)=gr(ix,iref);
                        end;                            
                    end;
                end;
           end;                                                                                    ! end v(x,z) condition
!
! Imaging Condition
!
            for ish=1:nsh;
                ihso=floor(((ish-1)*dsh+osh)/(dx));      
                for ixw=1:nkx;
                    if ixw-ihso < 1 || ixw-ihso > nkx || ixw+ihso < 1 || ixw+ihso > nkx;
                        continue;
                    else;
                        GG(iz,ish,ixw)=GG(iz,ish,ixw)+conj(ws(iz,ixw-ihso))*(wr(iz,ixw+ihso));             ! correlation to image
                    end;
                end;
            end;
        end;                                                                                            ! end depth loop
        disp([' DONE FREQUENCY #',num2str(iw),' OF ',num2str(nw),' / SHOT #',num2str(ixs)]);        
        sws(iw,:,:)=ws(:,:);
        swr(iw,:,:)=wr(:,:);
    end;                                                                                                ! end frequency loop

    image(:,:,ni:ni+nkx-1)=image(:,:,ni:ni+nkx-1)+GG;
    cputb=cputime;
    disp([' DONE SHOT # ',num2str(ixs),' OF ',num2str(nxs),' in ',num2str((cputb-cputa)/60),' min']);
end;                                                                                                    ! end shot loop

image=real(image);
sws(nw+1:nt,:,:)=0;
swr(nw+1:nt,:,:)=0;
sws=real(ifft(conj(sws)));
swr=real(ifft((swr)));
a=zeros(floor(nt/5)+1,nz,nkx);b=a;
for i=1:5:nt;
    b((i-1)/5+1,:,:)=swr(i,:,:);
    a((i-1)/5+1,:,:)=sws(i,:,:);
end;
clear swr sws;
swr=b;sws=a;
