program OW_arealsht

use sep 
use OW_ctrill
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
!use OW_datum
use OW_taper_areal
use OW_rndseed


implicit none

real					:: start, finish
integer					:: node,i,nv,ixs,j,ni,nsy,iw,ih,iphi,ix,izphi
real					:: owrk,phi,spec
real, allocatable				:: velint(:,:),taper(:),vel(:,:),rnd_phase(:)
real,allocatable				:: reflect(:,:),rand1(:)
complex, allocatable			:: recv(:,:,:),source(:,:,:),synshot(:,:)
complex,allocatable				:: G1(:,:),refl(:,:),tfld(:),G2(:,:)

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()


node=omp_get_num_threads()

write(0,*) 'node = ',node

start=0.;finish=0.
call owexec_time(start)
!---------------------------------------
! IO
!
 call owparam(.false.)

!----------------------------------------------------
! Memory Allocation
!
allocate(source(sou%h%n,sou%w%n,sou%x%n))
allocate(velint(v%z%n,v%xm%n))
allocate(vel(img%z%n,img%xm%n))
allocate(taper(img%xm%n)) 

call owexec_time(finish)
write(0,'(8x,a,f10.2,a)')' Init.phase - Total elapsed time:', (finish - start),' (sec)'

call sreed("sou",source,sou%h%n*sou%w%n*sou%x%n*8)

!---------------------------------------------------
!
! Read and extend velocities
!
call sreed("vel",velint,v%z%n*v%xm%n*4)
nv = floor((img%xm%o-v%xm%o)/v%xm%d)+1
vel(:,:)=velint(1:img%z%n,nv+1:nv+img%xm%n)

if (exist_file("velout")) then
   call to_history("n1",img%z%n,"velout");call to_history("o1",v%z%o,"velout");call to_history("d1",v%z%d,"velout"); 
   call to_history("n2",img%xm%n,"velout");call to_history("o2",img%xm%o,"velout");call to_history("d2",img%xm%d,"velout");
   call srite("velout",vel,4*v%z%n*img%xm%n)
end if
deallocate(velint)

if(exist_file("shot_in")) then
   allocate(recv(rec%h%n,rec%w%n,rec%x%n))
   call sreed("shot_in",recv,rec%w%n*rec%h%n*rec%x%n*8)
end if

!--------------------------------------------------
!
!
write(0,*) "UPWARD CONTINUES THE PREDEFINED WAVEFIELD"
write(0,*) "========================================="

open(unit=2,file='lixo')

call owtaper(tapertr,taper)

!--------------------------------------------------
!
! Generating the Synthesis Operator
!
!-----------------------------------------------------------------------------------------------      
allocate(reflect(img%z%n,img%xm%n))

do iphi=1,nphi
   reflect=0 
   phi=(iphi-1)*dphi+ophi
   if (exist_file("synth_op")) then
       allocate(G1(rec%w%n,img%xm%n))
       call sreed("synth_op",G1,img%xm%n*rec%w%n*8)
   else
       do ix=1,img%xm%n
          izphi=tan(phi*pi/180)*(ix-img%xm%n/2)
          if (dref-izphi > img%z%n) then
              reflect(img%z%n,ix)=0.
          elseif (dref-izphi < 1) then
              reflect(1,ix)=0.
          else
              reflect(dref-izphi,ix)=1.
          end if
       end do
       allocate(G1(rec%w%n,img%xm%n)); G1=0
       start=0.;finish=0. 
       call owexec_time(start)
       call owctrl_init(wbott,node)
!       call owctrl(rnd_phase,taper,reflect,source,vsal,vel,G1)
       call owctrl(taper,reflect,source,vsal,vel,G1)
       call owexec_time(finish)
!    allocate(G2(rec%w%n,img%xm%n));G2=G1
!    write(0,*) img%xm%n,imgxn,imgxn+padtraces
!    G1=G2(:,padtraces+1:imgxn+padtraces)
!    call srite("out",G1,8*rec%w%n*imgxn)
!       call srite("out",G1,8*rec%w%n*img%xm%n)
!    deallocate(G1);allocate(G1(rec%w%n,img%xm%n))
!    G1=G2;deallocate(G2)
   end if

!---------------------------------------------------
!
! Synthesizing areal shots
!
!---------------------------------------------------
   if(exist_file("shot_in")) then
   write(0,*) "Synthesizing areal shot for wvfld.dip = ", (iphi-1)*dphi+ophi
   write(0,*) "================================================="
   allocate(synshot(rec%w%n,nsh_img));synshot=0
   allocate(rand1(rec%w%n)); rand1=0.

!       spec=1/(cabs(source(1,:,1))+100*minval(cabs(source(1,:,1))))
       do ixs=1,rec%x%n
          owrk=(ixs-1)*rec%x%d+rec%x%o
          ni = (ixs-1)*floor(rec%x%d/img%xm%d) + 1
          nsy=floor(((ixs-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)+1
          write(0,*) nsy
          do iw=1,rec%w%n
             if (encode) then 
                call init_random_seed(ixs,rec%w%n,rand1,iw)
                G1(iw,nsy)=exp(cmplx(0,(2*pi*rand1(iw))))*G1(iw,nsy)
             end if
             do ih = 1,rec%h%n
                synshot(iw,ni+ih-1)=synshot(iw,ni+ih-1)+recv(ih,iw,ixs)*(G1(iw,nsy))!*spec
             end do
          end do
       end do
       call srite("synsh",synshot,8*rec%w%n*nsh_img)
       call srite("out",G1,8*rec%w%n*img%xm%n)
   end if
   deallocate(G1,synshot,rand1)
end do

deallocate(source,vel)
close(2)
call sep_end_prog()
end program OW_arealsht
