program OW_arealsht

use sep 
use OW_ctrillwei
use OW_type
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_fftw_areal
use OW_picksh
use OW_rndseed1
use OW_encoding

implicit none

real						:: start, finish,start1, finish1,a(2),b,c
real						:: rw(2)=0.,rz(2,2)=0.
integer						:: i,nv,ixs,j,ni,nsy,iw,ih,icomb,ix,stat,iz,ihor,n_enc_sog,iq,nq
integer						:: ind_trac,i_pos,i_ini,i_fin,igap
integer						:: jgap,jxs,nnn,hij,k,ib,ic,ireal,igap1
integer,allocatable					:: window(:,:),HD(:,:),q(:),prefer(:,:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:)
real,allocatable					:: randz(:,:,:,:)
real, allocatable					:: randw(:,:),rand4(:,:,:,:),rand1(:),vref(:,:,:),randx(:)
real,allocatable					:: hold(:,:),gol(:,:,:),randq(:,:),correl(:,:)
complex,allocatable					:: source(:),GGR(:,:),GGS(:,:),fun1(:)
complex, allocatable				:: GGG(:,:,:),wave_R(:,:,:),wave_S(:,:,:),tfld(:)
complex						:: phase

integer, external 					:: omp_get_thread_num,omp_get_num_threads,omp_get_num_procs,omp_set_num_threads

 call sep_init()
 call sep_begin_prog()

open(unit=2,file="logfile.txt")
open(unit=3,file="logfile1.txt")
open(unit=4,file="logfile2.txt")
open(unit=5,file="logfile3.txt")

start=0.;finish=0.
call owexec_time(start)

!---------------------------------------
! IO
!
 call owparam(.false.)

!--------------------------------------------------
! Encoding functions
!
start=0.;finish=0. 
call owexec_time(start)

!!!!!!!!!!!!!!!
!
! To be deleted 
! when working 
! with horizons
!
nhorz=2
!
!
if (encode) then
   allocate(randz(sou%w%n,sog%xm%n,nhorz,nreal))
   if(rand) then
      call init_random_seed4(randz)
      if (exist_file("correlations")) then
         if (allocated(tfld)) deallocate(tfld);allocate(tfld(sou%w%n))
         call owfftw(tfld)
         allocate(fun1(sou%w%n))
         allocate(correl(sog%xm%n,sou%w%n))
         tfld=cmplx(randz(:,1,1,1),0.)
         call sfftw_execute(planF)
         fun1=conjg(tfld/sqrt(1.*(sou%w%n)))
         do ixs=1,sog%xm%n
            tfld=cmplx(randz(:,ixs,1,1),0.)
            call sfftw_execute(planF)
            tfld=tfld/sqrt(1.*(sou%w%n))*fun1
            call sfftw_execute(planI)
            correl(ixs,:)=tfld
            correl(ixs,:)=cshift(correl(ixs,:),int((sou%w%n)/2))
         end do
         correl=correl/maxval(correl(1,:))*2*sou%w%n
         call to_history("n1",(sog%xm%n),"correlations"); call to_history("d1",1.,"correlations"); 
         call to_history("o1",1.,"correlations")
         call to_history("n2",sou%w%n,"correlations") ; call to_history("d2",2.,"correlations") ;
         call to_history("o2",-float(sou%w%n),"correlations")
         call to_history("n3",1,"correlations"); call to_history("d3",1.,"correlations");
         call to_history("o3",1.,"correlations")  
         call to_history("n4",1,"correlations"); call to_history("d4",1.,"correlations");
         call to_history("o4",0.,"correlations")
         call to_history("n5",1,"correlations"); call to_history("d5",1.,"correlations");
         call to_history("o5",0.,"correlations")  
         call to_history("esize",4,"correlations")
         call srite("correlations",correl,4*sog%xm%n*sou%w%n)
         call owfftwdestroy()
         deallocate(tfld)
      end if
!   else if (gold.and.not(perfect)) then
!      call owgoldw (randz)
   else if (gold.and.perfect) then
      call owgoldpw (randz,nhorz)
   else if (exist_file("InCode")) then
nhorz=1
      if (allocated(randz)) deallocate(randz); allocate(randz(sou%w%n,sog%xm%n,nhorz,nreal))
      call sreed("InCode",randz,4*sou%w%n*sog%xm%n*nreal*nhorz)
   end if
   if(exist_file("OutCode")) then
     call srite("OutCode",randz,4*sou%w%n*sog%xm%n*nreal*nhorz)
   end if
   call owexec_time(finish)
   write(0,*) "Finished generating encode functions in ",finish-start,"(s)"
end if

!---------------------------------------------------
!
! Read velocities
!
allocate(velint(v%z%n,v%xm%n))
call sreed("vel",velint,v%z%n*v%xm%n*4)
if (exist_file("vref")) then
   allocate(vref(nref,2,v%z%n))
   call sreed("vref",vref,v%z%n*2*nref*4)
end if

nv=int((xmin-v%xm%o)/v%xm%d)+1.5
allocate(vel(sog%z%n,nx))
vel(1:sog%z%n,:)=velint(1:sog%z%n,nv:nv+nx-1)
deallocate(velint)

!--------------------------------------------------
! Picks
!
if (exist_file("hor1")) then
   call owpick()
end if

!--------------------------------------------------
! Computing taper
!
allocate(taper(kxm%n))
call owtaper(int(padtraces/4),taper)

allocate(wave_S(sou%w%n,nx,nreal),wave_R(sou%w%n,nx,nreal));wave_S=0;wave_R=0
allocate(GGS(sou%w%n,kxm%n),GGR(sou%w%n,kxm%n))

call owctrlwei_init()

do ireal=1,nreal
   GGS=0;GGR=0
   if (exist_file("vref")) then
      call owexec_time(start)
      call owctrlwei(taper,vsal,vel,GGR,+1,randz(:,:,:,ireal),vref)
      call owexec_time(finish)
      write(0,3) "Receiver wavefield - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
      call owexec_time(start)
      call owctrlwei(taper,vsal,vel,GGS,-1,randz(:,:,:,ireal),vref)
      call owexec_time(finish)
      write(0,3) "Source wavefield   - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
   else
      call owexec_time(start)
      call owctrlwei(taper,vsal,vel,GGR,+1,randz(:,:,:,ireal))
      call owexec_time(finish)
      write(0,3) "Receiver wavefield - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
      call owexec_time(start)
      call owctrlwei(taper,vsal,vel,GGS,-1,randz(:,:,:,ireal))
      call owexec_time(finish)
      write(0,3) "Source wavefield   - realization # ",ireal,' FROM ',nreal,' in ',finish-start,' s'
   end if
   wave_R(:,:,ireal)=GGR
   wave_S(:,:,ireal)=GGS
end do
3  format(A35,i4,A6,i4,A4,f6.2,a2)

!do iw=1,rec%w%n
!   phase = cmplx(cos(pi/4.), -sin(pi/4.)) / sqrt(rec%w%all(iw))
!   wave_R(iw,:,:)=wave_R(iw,:,:)*(phase)
!   wave_S(iw,:,:)=wave_S(iw,:,:)*conjg(phase)
!end do
!wave_S=wave_S*phase
!wave_S=wave_S*conjg(phase)
!wave_S=wave_S*phase
!
!       recv sour
! 2.T = conj ----
! 3.T = conj conj
! 4.T = ---- ----
! 5.T = ---- conj

call srite("rec_out",wave_R,8*sou%w%n*nx*nreal)
call srite("sou_out",wave_S,8*sou%w%n*nx*nreal)

close(2)
close(3)
close(4)
close(5)

deallocate(GGR,GGS,wave_R,wave_S)
deallocate(vel)

call sep_end_prog()
end program OW_arealsht
