program greenxs4

use sep 
use OW_type
use igrad1

integer						:: time1(8)
real						:: start_time, finish_time

integer						:: ixs,iw,iz,ikx,iref,ixw,ish,ino
integer						:: nref,ni,nv,isp,inode,nnode,nr
real						:: xmin,xmax,ymin,ymax,zmin,zmax
real						:: sign, pi
real						:: minv, maxv, dvel, owrk
real, allocatable				:: vref(:,:,:),vel(:,:),velint(:,:),sqa0(:),reflect(:,:)
complex, allocatable				:: g0(:),gx0(:,:)
complex, allocatable				:: source(:,:,:),GG(:,:,:),G(:,:,:)
complex, allocatable				:: tfield(:)

type (axis) 					:: kxm, h, w, x, xm, z, sh
type (dat_grid)					:: rec, sou
type (mod_grid)					:: img
type (velo)					:: v

integer, external 				:: omp_get_thread_num
integer, external 				:: omp_get_num_threads

integer, parameter   				:: FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
integer, parameter   				:: FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
integer, parameter   				:: FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
integer, parameter   				:: FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
integer, parameter   				:: FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
integer, parameter   				:: FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
integer, parameter   				:: FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
integer, parameter   				:: FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536
 
!integer*8, pointer,private   			:: planF,planI
!integer*8, private	   			:: planF,planI
integer*8		   			:: planF,planI  

 call sep_init()
 call sep_begin_prog()

start_time=0.;finish_time=0.
call date_and_time(values=time1)
start_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)

!
! IO
!

 call from_history("n1",rec%h%n);	call from_history("d1",rec%h%d);	call from_history("o1",rec%h%o)
 call from_history("n3",rec%w%n);	call from_history("d3",rec%w%d);	call from_history("o3",rec%w%o)
 call from_history("n4",rec%x%n);	call from_history("d4",rec%x%d);	call from_history("o4",rec%x%o)  

 call from_aux("sou","n1",sou%h%n);	call from_aux("sou","d1",sou%h%d);	call from_aux("sou","o1",sou%h%o)
 call from_aux("sou","n3",sou%w%n);	call from_aux("sou","d3",sou%w%d);	call from_aux("sou","o3",sou%w%o)
 call from_aux("sou","n4",sou%x%n);	call from_aux("sou","d4",sou%x%d);	call from_aux("sou","o4",sou%x%o)  

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

 if (sou%w%n .ne. rec%w%n) then
    call seperr("# source freq. .ne. # receiver freq.")
 end if
 if (sou%w%d .ne. rec%w%d) then
    call seperr("source freq.interval .ne. receiver freq.interval")
 end if
 if (sou%w%o .ne. rec%w%o) then
    call seperr("source freq.origin .ne. receiver freq.origin")
 end if

 !!  get parameters
  call from_param("nrefvel",nref) 
  call from_param("xmin",xmin) 
  call from_param("xmax",xmax)
  call from_param("ymin",ymin,0.)
  call from_param("ymax",ymax,0.)
  call from_param("zmin",zmin,0.)
  call from_param("zmax",zmax)

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

 if (zmin < v%xm%o ) call seperr("Min.z-image space < Min.z-velocity")
 if (zmax > v%z%o+(v%z%n-1)*v%z%d ) call seperr("Max.z-image space > Max.z-velocity")

 if (rec%h%o + rec%x%o < xmin ) call seperr("Min.x-receiver < Min.x-image space")
 if (rec%h%o+(rec%h%n-1)*rec%h%d + rec%x%o+(rec%x%n-1)*rec%x%d > xmax ) call seperr("Max.x-receiver < Max.x-image space")

 if (sou%h%o + sou%x%o < xmin ) call seperr("Min.x-source < Min.x-image space")
 if (sou%h%o+(sou%h%n-1)*sou%h%d + sou%x%o+(sou%x%n-1)*sou%x%d > xmax ) call seperr("Max.x-source < Max.x-image space")

 if (rec%h%d .ne. v%xm%d) call seperr("x-offset interval .ne. velocity x-interval")

 img%xm%n=floor((xmax-xmin)/v%xm%d)+1
 img%xm%d=v%xm%d
 img%xm%o=xmin

 img%z%n=floor((zmax-zmin)/v%z%d)+1
 img%z%o=zmin
 img%z%d=v%z%d

 call to_history("n1",rec%w%n) ; call to_history("d1",rec%w%d) ; call to_history("o1",rec%w%o)
 call to_history("n2",rec%h%n); call to_history("d2",rec%h%d); call to_history("o2",rec%h%o)
 call to_history("n3",rec%x%n); call to_history("d3",rec%x%d); call to_history("o3",rec%x%o)  
 call to_history("n4",1); call to_history("d4",1); call to_history("o4",1)
 call to_history("n5",1); call to_history("d5",1); call to_history("o5",1)  
 call to_history("esize",8)

nnode=omp_get_num_threads()

 allocate(source(sou%h%n,sou%w%n,sou%x%n))
 allocate(velint(v%z%n,v%xm%n))
 allocate(reflect(img%z%n,img%xm%n))
 allocate(vref(nref,img%z%n,img%xm%n))
 allocate(vel(v%z%n,rec%h%n))

sign=-1;
pi=4*atan(1.)
open(unit=2,file="lixo")
!--------------------------------------------------
! Frequency vector
!
! no zero frequency component
rec%w%o   = 2*pi * rec%w%o
rec%w%d   = 2*pi * rec%w%d
allocate(rec%w%all(rec%w%n))
do iw = 1,rec%w%n
   rec%w%all(iw) = rec%w%o + (iw - 1) * rec%w%d
end do
!--------------------------------------------------
! Wavenumber vector
!
kxm%n = rec%h%n
kxm%d = 1./(kxm%n*rec%h%d)
kxm%o = 1./(2*rec%h%d)-kxm%d
allocate(kxm%all(kxm%n))
do ikx=1,kxm%n
   kxm%all(ikx) = mod(kxm%o+(ikx-1)*kxm%d,1./rec%h%d)-kxm%o
end do
kxm%all = 2*pi * kxm%all

!--------------------------------------------------
! Setting FFTW plans
!

allocate(tfield(kxm%n))
call sfftw_plan_dft(planF,1,kxm%n,tfield,tfield,fftW_FORWARD,fftW_PATIENT)
call sfftw_plan_dft(planI,1,kxm%n,tfield,tfield,fftw_BACKWARD,fftW_PATIENT)

!allocate(trace(
!call sfftw_plan_dft(planI,1,kxm%n,tfield,tfield,fftw_BACKWARD,fftW_PATIENT)
!--------------------------------------------------
! Reading source fctn, velocity

 call sreed("sou",source,sou%h%n*sou%w%n*sou%x%n*8)
 call sreed("vel",velint,v%z%n*v%xm%n*4)

nr=floor((img%xm%o-v%xm%o)/img%xm%d)+1
do ix=1,img%xm%n
   stat = igrad1_lop (.true.,.false.,reflect(:,ix),velint(:,nr+ix-1))
end do


!--------------------------------------------------
! Pre-compute vref
!

v%z%n=img%z%n
vref=0
write(0,*) "owrk_max ",(rec%h%n+floor(rec%x%d/rec%h%d)*rec%x%n)*rec%h%d+rec%x%o+rec%h%o,"nvmax ",floor(((img%xm%n-1)*rec%h%d+rec%x%o+rec%h%o-v%xm%o)/v%xm%d)+1,&
            "  ",floor(((img%xm%n-1)*rec%h%d+rec%x%o+rec%h%o-v%xm%o)/v%xm%d)+rec%h%n

nsh=rec%h%n+floor(rec%x%d/rec%h%d)*rec%x%n+1
do ixs=1,nsh
    owrk=(ixs-1)*rec%h%d+rec%x%o+rec%h%o
    if (owrk < img%xm%o .or. owrk + (rec%h%n-1)*rec%h%d+rec%h%o > img%xm%o+img%xm%d*(img%xm%n-1)) then
       write(0,*) " SHOT # ",ixs ; call seperr(' OUT OF WORKING SPACE')
    end if
! grab velocity in the workspace       
    nv=floor((owrk-v%xm%o)/v%xm%d)+1;
!    write(0,*) ixs,nv,nv+rec%h%n-1,img%xm%n,v%xm%n
    vel=0;
    vel=velint(:,nv:nv+rec%h%n-1);
!
    do iz=1,v%z%n
        if (nref .eq. 1) then
           vref(nref,iz,ixs)=sum(vel(iz,:))/rec%h%n
        else
           minv=minval(vel(iz,:))
           maxv=maxval(vel(iz,:))
           dvel=(maxv-minv)/(nref-1) 
           do iref=1,nref
              vref(iref,iz,ixs)=minv+(iref-1)*dvel
           end do 
        end if
    end do 
end do 

call date_and_time(values=time1)
finish_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)
write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish_time - start_time)/60./60.,' (hrs)'

isp = abs(floor(rec%h%o/rec%h%d)) + 1

sign=-1;

!
!
!
    write(0,*) "COMPUTE GREENS FUNCTIONS AND GENERATE DATA"
    write(0,*) "======= ====== ========= === ======== ===="

	if (allocated(GG)) deallocate(GG);allocate(GG(rec%w%n,rec%h%n,rec%x%n));GG=0.
    
    do iw = 1,rec%w%n;                                                                                          ! frequency loop
! exec.time
    call date_and_time(values=time1)
    start_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)

    if (allocated(G)) deallocate(G);		allocate(G(v%z%n,rec%h%n,nsh));	G=0.
    if (allocated(g0)) deallocate(g0);		allocate(g0(kxm%n));		g0=0.
    if (allocated(gx0)) deallocate(gx0);	allocate(gx0(kxm%n,nref));		gx0=0.
    if (allocated(tfield)) deallocate(tfield);	allocate(tfield(kxm%n));		tfield=0.
    if (allocated(sqa0)) deallocate(sqa0);	allocate(sqa0(kxm%n));		sqa0=0.

       do ixs=1,nsh                                                                                       ! shot loop

! workspace origin 
            owrk=(ixs-1)*rec%h%d+rec%x%o+rec%h%o
! grab velocity in the workspace    
   
            nv = floor((owrk-v%xm%o)/v%xm%d)+1
            ni = nv-floor((img%xm%o-v%xm%o)/img%xm%d)+1

            vel=0;
            vel(:,:)=velint(:,nv:nv+rec%h%n-1);

            ixp=floor(((rec%x%o-v%xm%o)+rec%x%d*(ixs-1))/v%xm%d)+1;
!            write(0,*) ixp, v%xm%n
            
            if (ixp>v%xm%n) then
                write(0,*) " SHOT # ",ixs
                call seperr(" OUT OF BOUNDS")
            end if

            G(1,isp,ixs)=source(1,iw,1)
            g0(isp)=G(1,isp,ixs)
            write(0,*) maxval(cabs(g0))

            do iz=2,img%z%n;                                                                                    ! depth loop

                minv=minval(vel(iz,:))
                maxv=maxval(vel(iz,:))

                vxz:if (minv .eq. maxv) then
!      Operations in the wavenumber domain
!      phase shift
                   sqa0 = rec%w%all(iw)**2 / minv**2 - kxm%all**2
                   evanescent1:where (sqa0 > 0)
                       g0=g0*exp(cmplx(0,-sign*v%z%d*sqrt(sqa0)))
                   elsewhere
                       g0=g0*exp(-sqrt(-sqa0)*v%z%d)
                   end where evanescent1
! Operations in x            

                   tfield=g0; call sfftw_execute(planI); g0=tfield/sqrt(1.*kxm%n)
                   G(iz,:,ixs)=g0

                else  

                   velref: do iref=1,nref
!      Operations in the wavenumber domain
!      phase shift

                        sqa0 = rec%w%all(iw)**2 / vref(iref,iz,ixs)**2 - kxm%all**2;
                        evanescent2:where (sqa0 > 0)
                            gx0(:,iref)=gx0(:,iref)*exp(cmplx(0,-sign*v%z%d*sqrt(sqa0)))
                        elsewhere
                            gx0(:,iref)=gx0(:,iref)*exp(-sqrt(-sqa0)*v%z%d)
                        end where evanescent2
! Operations in x            

                        tfield=gx0(:,iref); call sfftw_execute(planI); gx0(:,iref)=tfield/sqrt(1.*kxm%n)

!      split step
                        gx0(:,iref) = gx0(:,iref) * exp(cmplx(0,-sign*v%z%d*rec%w%all(iw))*(1./vel(iz,:)-1./vref(iref,iz,ixs)));
                   end do velref

                do iref=1,nref-1
                interpwvfld: where (vel(iz,:) >= vref(iref,iz,ixs) .and. vel(iz,:) < vref(iref+1,iz,ixs))
                     G(iz,:,ixs)=gx0(:,iref)
                     end where interpwvfld
                end do
                where (vel(iz,:) >= vref(nref,iz,ixs))
                      G(iz,:,ixs)=gx0(:,nref)
                end where 

!                do ix=1,kxm%n
!                    do iref=1,nref-1
!                        if (vref(iref,iz,ixs) .ne. vref(iref+1,iz,ixs)) then
!                            if (vel(iz,ix) >= vref(iref,iz,ixs) .and. vel(iz,ix) < vref(iref+1,iz,ixs)) then
!                                f=(vel(iz,ix) - vref(iref,iz,ixs))/vel(iz,ix)
!                                g=1-f
!                                G(iz,ix,ixs)=g*gx0(ix,iref)+f*gx0(ix,iref+1)
!!                                G(iz,ix,ixs)=gx0(ix,iref)
!                            elseif (vel(iz,ix) >= vref(nref,iz,ixs)) then
!                                G(iz,ix,ixs)=gx0(ix,nref)
!                            end if
!                        else
!                            G(iz,ix,ixs)=gx0(ix,iref)
!                        end if                            
!                    end do
!                end do
           end if vxz
       end do
  end do

  write(2,*) iz,ixs,maxval(cabs(G(iz,:,ixs)))
!
! Convolutions to generate data
!

       do ixs=1,floor(rec%x%d/rec%h%d),rec%x%n;
           owrk=(ixs-1)*rec%h%d+rec%x%o+rec%h%o
! grab reflectivity in the workspace       
           nv = floor((owrk-v%xm%o)/v%xm%d)+1
           ni = nv-floor((img%xm%o-v%xm%o)/img%xm%d)+1
           do jxs=1,rec%h%n
               do ix=ni,ni+rec%h%n-1;
                   do iz=1,img%z%n;
                       GG(iw,jxs,ixs)=GG(iw,jxs,ixs)+G(iz,jxs,ixs)*G(iz,jxs,ixs+jxs)*reflect(iz,ix);      ! convolve to generate data
                   end do
               end do
           end do
       end do

    call date_and_time(values=time1)
    finish_time = time1 (5) * 3600 + time1 (6) * 60 + time1 (7) + 0.001 * time1 (8)
    write(0,*) " DONE FREQUENCY # ",iw," OF ",rec%w%n," in ",(finish_time - start_time)/60," (min)"
5   format(a18,2(i3,a4),f6.2,a6)

  end do                                                                                                    ! end frequency loop

 call srite("out",GG,8*rec%w%n*rec%h%n*rec%x%n)
 close(2)
 call sep_end_prog()
 end program greenxs4
