module oneway_mig
use sep
use myfft_mod
use ctaper_mod
use sort_mod
use vel_mod
use findsr_mod
!use global_defi_mod
!use global_mod

implicit none

contains
!==============================================
subroutine oneway(data,x, xs, xr, z, t, v, vnum, source, nxpad, imin_z, imax_z, jz, jx, jxr, jxs, reflect,restart,tape,opt,what)

complex , dimension(:,:,:)      , pointer :: data
complex , dimension(:,:,:)   ,allocatable :: fdata
complex , dimension(:,:)        , pointer :: mig
complex , dimension(:,:)     ,allocatable :: img,tmp_mig
real    , dimension(:,:)        , pointer :: v
real    , dimension(:)       ,allocatable :: vtmp
real    , dimension(:,:)        , pointer :: reflect
real    , dimension(:,:)     ,allocatable :: vel
complex , dimension(:,:,:)   ,allocatable :: S
complex , dimension(:,:  )   ,allocatable :: R
complex , dimension(:,:)     ,allocatable :: green_s, green_r
complex , dimension(:)          , pointer :: source
complex , dimension(:)       ,allocatable :: wavelet,wave,fwavelet
real    , dimension(:)          , pointer :: x, z, xs, xr, t
real    , dimension(:)       ,allocatable :: tmp 
real    , dimension(:)          , pointer :: wall,kxall
real    , dimension(:)       ,allocatable :: vref
integer , dimension(:)       ,allocatable :: ind_vref        
complex , dimension(:)       ,allocatable :: tmp_wfld        
complex , dimension(:,:,:,:) ,allocatable :: wfld_d, wfld_u        
real    , dimension(:)       ,allocatable :: sqa        
real    , dimension(:,:,:)   ,allocatable :: data2        
real    , dimension(:,:,:)   ,allocatable :: data3        

real                    :: dx, dy, dz, ox, oy, oz  
real                    :: dxr, dxs, oxr, oxs 
real                    :: dt, ot 
real                    :: dw, ow, dkx, okx  
real                    :: tape

integer                 :: ix, iz, iw, it, itt, ik, iv
integer                 :: nx, ny, nz
integer                 :: nxr, nxs
integer                 :: nw, nkx, nt
integer                 :: nxv
integer                 :: vnum, nref
integer                 :: nxpad, imax_z, imin_z
integer                 :: jz, jx, jxr, jxs
integer                 :: restart, what, opt
integer                 :: nsour2
integer                 :: iww, niw
integer                 :: i, j, k, ii, jj, kk
integer                 :: ir, is
integer                 :: ngpts
integer                 :: TID, OMP_GET_THREAD_NUM

integer                 :: ixp, ixg, ixs, ixr, iixs, iixr
integer                 :: sign
real                    :: pi
character(len=256)      :: vel_="vel2"
INTEGER, EXTERNAL       :: omp_get_num_procs 
INTEGER                 :: inode, nnode, i1, i2, thread
real                    :: xxs, xxr 
character(len=256)      :: sou1_="sou1", sou2_="sou2", fdata_="fdata", mig_="mig"
logical                 :: invert
!=============================================================
pi=3.1415937; sign=1;    
!=============================================================
!=thread stuff
nnode=omp_get_num_procs()
thread=omp_get_thread_num()
write(*,*)'nnonde,thread'
write(*,*)nnode,thread

call omp_set_num_threads(floor(nnode*.98))
!call omp_set_num_threads(floor(1.0))

sign=1;
!=============================================================
!=================  Initial computations  ====================
!=============================================================
!=model grid....vel grid...image grid...
        nx=size(x)  ;   nt=size(t)  ;   nz=size(z)  ;     ny=1;  
        dx=x(2)-x(1);   dt=t(2)-t(1);   dz=z(2)-z(1);     dy=dx;
        ox=x(1)     ;   ot=t(1)     ;   oz=z(1)     ;     oy=0;
!=============================================================
!=source/receiver axis
        nxs=size(xs) ; nxr=size(xr) ;
        if (nxs > 1)then; dxs=xs(2)-xs(1); oxs=xs(1); else; dxs=1; oxs=0.; endif;
        if (nxr > 1)then; dxr=xr(2)-xr(1); oxr=xr(1); else; dxr=1; oxr=0.; endif;
!=============================================================
!=Frequency axes
!=temporal freq.
        nw = floor(nt/2.); dw = 1./(nt*dt); ow=0;  allocate(wall(nw));
        do iw = 1, nw; wall(iw) = ow  + (iw - 1) * dw; enddo; wall = 2*pi * wall;
!==
!=spatial freq.
        nkx = nx+2*nxpad; dkx = 1./(nkx*dx); okx = 1./(2*dx)-dkx; allocate(kxall(nkx));
        do ik=1,nkx; kxall(ik) = mod(okx+(ik-1)*dkx,1./dx)-okx; enddo; kxall = 2*pi * kxall;
!=============================================================
!=propagation grid...
        nxv=size(v(1,:));
        if (nxv<nkx)then 
            allocate (vel(nz,nx+2*nxpad)) 
            do iv=1,nxpad; 
                vel(:,iv)=v(:,1);
                vel(:,nxpad+nxv+iv)=v(:,size(v(1,:))); 
            enddo; 
            vel(:,nxpad+1:nxpad+nxv)=v;
        else
            allocate (vel(nz,nx)) 
            vel=v;
        endif;
!=============================================================
!=pad and position wavelet velocities
    allocate(wavelet(nt)); wavelet=source;
    call FFT1DT(wavelet, 1); allocate(wave(nw)); wave=wavelet(1:nw);
    call ctaper1d(wave,wave,dw, nw, 0, 0.1, 2);
!call auxinout(sou2_)
!call auxputch("n1", "i", nw , sou2_);  
!call auxputch("d1", "r", dw , sou2_); 
!call auxputch("o1", "r", ow , sou2_); 
!call auxputch("label1", "s", "w" , sou2_)
!call auxputch("esize" , "i",  8  , sou2_)
!call srite("sou2",wave , 8*nw)
    
    wave=wave*exp(cmplx(0.,wall*nw*dt))

!call auxinout(fdata_)
!call auxputch("n1", "i", nw , fdata_);  
!call auxputch("d1", "r", dw , fdata_); 
!call auxputch("o1", "r", ow , fdata_); 
!call auxputch("label1", "s", "w" , fdata_)
!call auxputch("esize" , "i",  8  , fdata_)
!call srite("fdata",wave , 8*nw)
!=compute the useful freq. 
    iww=0;
    do niw=floor(nw*0.1),nw
        if (abs(wave(niw))<=maxval(abs(wave))*.005)then; iww=iww+1; if (iww>=floor(nw*0.02)) exit; endif
    enddo
!=============================================================
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt))
if (allocated(fdata)) deallocate(fdata); allocate(fdata(nw,nxr,nxs))
tmp_wfld=.0;
do iixs=1,nxs
        do iixr=1,nxr
        write(*,*)iixs,iixr
                tmp_wfld=data(:,iixr,iixs); 
                call FFT1DT(tmp_wfld, 1);
                !fdata(:,iixr,iixs)= conjg(tmp_wfld(1:nw))
                fdata(:,iixr,iixs)= conjg(tmp_wfld(1:nw))
        enddo
enddo
call auxinout(fdata_)
call auxputch("n1", "i", nw , fdata_); call auxputch("n2", "i", nxr , fdata_); call auxputch("n3", "i", nxs , fdata_); 
call auxputch("d1", "r", dw , fdata_); call auxputch("d2", "r", dxr , fdata_); call auxputch("d3", "r", dxs , fdata_); 
call auxputch("o1", "r", ow , fdata_); call auxputch("o2", "r", oxr , fdata_); call auxputch("o3", "r", oxs , fdata_); 
call auxputch("label1", "s", "w" , fdata_)
call auxputch("esize" , "i",  8  , fdata_)
call srite("fdata",fdata, 8*nw*nxr*nxs)

!=============================================================
!%initial:w

allocate(mig(nz,nkx),tmp_mig(nz,nkx)); tmp_mig=0.

allocate(wfld_d(nz,nx,nt,1),wfld_u(nz,nx,nt,1),tmp(nxr)); wfld_d=0; tmp=xr

!$OMP parallel do private(inode,iw,ixp,ix,iz,iv,ik,nref,sqa,tmp_wfld,green_s,green_r,S,R,vtmp,vref,ind_vref,ixs,ixr,xxs,xxr)  schedule(dynamic,1) ordered 
do iw = 1,niw;      ! frequency loop
        inode=omp_get_thread_num();
        write(*,*)'freq. num.: ', iw ,' of ',niw,' inode ', inode;         
        if(abs(wave(iw))<=maxval(abs(wave))*.0001)cycle 

        if (allocated(S)) deallocate(S); allocate(S(nz,nkx,nxs)); S=0
        if (allocated(R)) deallocate(R); allocate(R(nz,nkx    )); R=0
        do ixs=1, nxs;                                                      !shot loop
            ixp=floor(((oxs-ox)+dxs*(ixs-1.))/dx)+1;
            if (ixp>nkx) then; call seperr('source out of velocity model'); endif;
            S(1,ixp+nxpad,ixs)=wave(iw);
            R(1,  1+nxpad:nxr+nxpad )=fdata(iw,:,ixs) ;

!            R(1,ixp+nxpad)=wave(iw) ;

!            write(*,*)S(1,ixp+nxpad,ixs), R(1,ixp+nxpad)
            !R(1,  1+nxpad:nxr+nxpad)=fdata(iw,:,ixs);
            do iz=2,nz;                                                 !depth loop
               if (allocated(vtmp)) deallocate(vtmp); allocate(vtmp(nkx))
               if (allocated(ind_vref)) deallocate(ind_vref); allocate(ind_vref(nkx))
               if (allocated(vref)) deallocate(vref); allocate(vref(vnum))
               vtmp=vel(iz,:); 
               call  velocity(vref,vtmp, ind_vref, vnum, nref); 
               if (allocated(green_s )) deallocate(green_s); allocate(green_s(nref,nkx))
               if (allocated(green_r )) deallocate(green_r); allocate(green_r(nref,nkx))
               if (allocated(sqa     )) deallocate(sqa    ); allocate(sqa    (nref))
               if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nkx))

               tmp_wfld=S(iz-1,:,ixs); call FFT1DX(tmp_wfld, 1); S(iz,:,ixs)=tmp_wfld;
               tmp_wfld=R(iz-1,:    ); call FFT1DX(tmp_wfld, 1); R(iz,:    )=tmp_wfld;
               do iv=1, nref; green_s(iv,:)=S(iz,:,ixs);enddo
               do iv=1, nref; green_r(iv,:)=R(iz,:    );enddo

               do iv=1,nref                     !v(x,z)
                  do ik = 1,nkx;                !wavenumber loop
                     sqa(iv) = wall(iw)**2 / vref(iv)**2 - kxall(ik)**2;
                     if (sqa(iv) > 0)then       !evanescent condition v0
                        green_s(iv,ik)=green_s(iv,ik)*exp(cmplx(0,   sign*dz*sqrt(sqa(iv))));
                        green_r(iv,ik)=green_r(iv,ik)*exp(cmplx(0,   sign*dz*sqrt(sqa(iv))));
                     else
                        green_s(iv,ik)=0.;
                        green_r(iv,ik)=0.;
                     endif
                  enddo                          !end wavenumber loop                        

! split step  
                  tmp_wfld=green_s(iv,:); call FFT1DX(tmp_wfld, -1); green_s(iv,:)=tmp_wfld
                  tmp_wfld=green_r(iv,:); call FFT1DX(tmp_wfld, -1); green_r(iv,:)=tmp_wfld

                  green_s(iv,:)=green_s(iv,:) * exp(cmplx(0,   sign*dz*wall(iw))*(1./vel(iz,:)-1./vref(iv)));    
                  green_r(iv,:)=green_r(iv,:) * exp(cmplx(0,   sign*dz*wall(iw))*(1./vel(iz,:)-1./vref(iv)));   

                  tmp_wfld=green_s(iv,:); call ctaper1d(tmp_wfld,tmp_wfld, dx, nkx, floor(nxpad*0.), tape, opt); green_s(iv,:)=tmp_wfld
                  tmp_wfld=green_r(iv,:); call ctaper1d(tmp_wfld,tmp_wfld, dx, nkx, floor(nxpad*0.), tape, opt); green_r(iv,:)=tmp_wfld
               enddo
               !write(*,*)maxval(real(green_r)),maxval(real(green_s));       

               xxs = oxs + (1-1)*dxs; call findsr(iixs,xxs,tmp);

               do ix=1,nkx
                  S(iz,ix,ixs)=green_s(ind_vref(ix),ix);
                  R(iz,ix    )=green_r(ind_vref(ix),ix);
               enddo
!              write(*,*)size(green_s(:,1+nxpad:nx+nxpad)),size(wfld_d(:,:,iw,1))
!              wfld_d(:,:,iw,1)=S(:,1+nxpad:nx+nxpad,iixs);
                if (ixs==2)then
               wfld_d(:,:,iw,1)=S(:,1+nxpad:nx+nxpad, ixs);
               wfld_u(:,:,iw,1)=R(:,1+nxpad:nx+nxpad);
                endif
               tmp_mig(iz,:)=tmp_mig(iz,:)+real(R(iz,:)*(S(iz,:,ixs)));
               !write(*,*)maxval(real(wfld_d)),maxval(real(wfld_u));       
               !write(*,*)sum(real(wfld_d)),sum(real(wfld_u));      
               !write(*,*)maxval(real(S)),maxval(real(R)),maxval(real(tmp_mig));       
           enddo                !end depth loop
        enddo                   !end shot loop
enddo                       !end frequency loop
!$OMP end parallel do   

if(allocated(data2)) deallocate(data2); allocate(data2(nx,nz,nt))
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
!call omp_set_num_threads(floor(1.0))
!$OMP parallel do private(iz,ix,tmp_wfld)   
do iz=1,nz
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
inode=omp_get_thread_num();
write(*,*)'iz, inode:',iz,wfld_d(iz,1,54,1),  wfld_u(iz,1,54,1),inode
        do ix=1,nx
                tmp_wfld=wfld_d(iz,ix,:,1);
                call FFT1DT(tmp_wfld, -1);
                data2(iz,ix,:)=tmp_wfld 
        enddo
enddo
!$OMP end parallel do

if (allocated(data3)) deallocate(data3); allocate(data3(nz,nx,ceiling(nt/10.)))

itt=1
do it=1,nt,10
       data3(:,:,itt)=data2(:,:,it) 
        itt=itt+1
        write(*,*)it,itt,nt
end do

call auxinout(sou1_)
call auxputch("n1", "i", nz , sou1_); call auxputch("n2", "i", nx , sou1_); call auxputch("n3", "i", ceiling(nt/10.), sou1_); 
call auxputch("d1", "r", dz , sou1_); call auxputch("d2", "r", dx , sou1_); call auxputch("d3", "r", dt*10 , sou1_);
call auxputch("o1", "r", oz , sou1_); call auxputch("o2", "r", ox , sou1_); call auxputch("o3", "r", ot , sou1_);
call auxputch("label1", "s", "nz" , sou1_)
call auxputch("label2", "s", "nx" , sou1_)
call auxputch("label3", "s", "t" , sou1_)
call auxputch("esize" , "i",  4  , sou1_)
call srite("sou1",data3 , 4*ceiling(nt/10.)*nz*nx)
!deallocate(wfld_d,data2,data3)

write(*,*)'here------------------>1'

if(allocated(data2)) deallocate(data2); allocate(data2(nx,nz,nt))
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
!call omp_set_num_threads(floor(1.0))
!$OMP parallel do private(iz,ix,tmp_wfld)   
do iz=1,nz
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
inode=omp_get_thread_num();
write(*,*)'iz, inode:',iz,wfld_d(iz,1,54,1),  wfld_u(iz,1,54,1),inode
        do ix=1,nx
                tmp_wfld=wfld_u(iz,ix,:,1);
                call FFT1DT(tmp_wfld, -1);
                data2(iz,ix,:)=tmp_wfld 
        enddo
enddo
!$OMP end parallel do

if (allocated(data3)) deallocate(data3); allocate(data3(nz,nx,ceiling(nt/10.)))

itt=1
do it=1,nt,10
       data3(:,:,itt)=data2(:,:,it) 
        itt=itt+1
        write(*,*)it,itt,nt
end do

call auxinout(sou2_)
call auxputch("n1", "i", nz , sou2_); call auxputch("n2", "i", nx , sou2_); call auxputch("n3", "i", ceiling(nt/10.), sou2_); 
call auxputch("d1", "r", dz , sou2_); call auxputch("d2", "r", dx , sou2_); call auxputch("d3", "r", dt*10 , sou2_);
call auxputch("o1", "r", oz , sou2_); call auxputch("o2", "r", ox , sou2_); call auxputch("o3", "r", ot , sou2_);
call auxputch("label1", "s", "nz" , sou2_)
call auxputch("label2", "s", "nx" , sou2_)
call auxputch("label3", "s", "t" , sou2_)
call auxputch("esize" , "i",  4  , sou2_)
call srite("sou2",data3 , 4*ceiling(nt/10.)*nz*nx)
deallocate(wfld_d,data2,data3)





write(*,*)'here------------------>2'


write(*,*)maxval(real(S)),maxval(real(tmp_mig)),maxval(real(mig));       
allocate(mig(nz,nx))
write(*,*)maxval(real(S)),maxval(real(tmp_mig)),maxval(real(mig));       


mig=tmp_mig(:,1+nxpad:nx+nxpad);
write(*,*)'writing img.....'
call auxinout(mig_)
call auxputch("n1", "i", nz, mig_); call auxputch("n2", "i", nz, mig_); call auxputch("n3", "i", ny, mig_); 
call auxputch("d1", "r", dx, mig_); call auxputch("d2", "r", dx, mig_); call auxputch("d3", "r", dy, mig_); 
call auxputch("o1", "r", oy, mig_); call auxputch("o2", "r", oy, mig_); call auxputch("o3", "r", oy, mig_); 
call auxputch("label1", "s", "z" , mig_)
call auxputch("label2", "s", "x" , mig_)
call auxputch("label3", "s", "y" , mig_)
call auxputch("esize" , "i",  4  , mig_)
call srite(mig_, real(mig), 4*nz*nx*ny)

    
    if (invert/=1)then
    else
        write(*,*)'.......Inverting ur migrated image size: ',size(x), ' by ', size(z);
    endif
    
    
end subroutine
end module
