program OWplane_areal

use sep 
use OW_type
use OW_exec_time

implicit none

real						:: start,finish,start1,finish1,v0,pi,oxplane,owrk,dx,alpha
integer						:: iang,iw,idip,ih,iz,nxplane,iowrk,node
complex, allocatable				:: recv(:,:),source(:,:),source1(:),phase1(:)
complex, allocatable				:: recv_pl(:,:,:),source_pl(:,:,:)
complex                                                     :: phase
type (axis) 					:: ang,h,w,xd
type (dat_grid)	 				:: rec,sou
integer, external 					:: omp_get_thread_num,omp_get_num_threads

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)

 call from_param("ango", ang%o, -30.); call from_param("angd", ang%d, 0.5); call from_param("angn", ang%n, int(2*abs(ang%o)/ang%d)+1 )
 call from_param("v0",v0,1500.)
 call from_param("node",node,9999)

if (node.eq.9999) then
  !$OMP PARALLEL
  node = omp_get_num_threads()
  !$OMP END PARALLEL
write(0,*) " NUM.THREADS = ",node
end if

 call from_aux("rec","n1",rec%h%n);	call from_aux("rec","d1",rec%h%d);call from_aux("rec","o1",rec%h%o)
 call from_aux("rec","n3",rec%w%n);	call from_aux("rec","d3",rec%w%d);call from_aux("rec","o3",rec%w%o)
 call from_aux("rec","n4",rec%xd%n);call from_aux("rec","d4",rec%xd%d);call from_aux("rec","o4",rec%xd%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)

 nxplane=rec%h%n+int((rec%xd%n-1)*(rec%xd%d/rec%h%d))
 oxplane=rec%h%o+rec%xd%o
 pi=4*atan(1.)

 call to_history("n1", nxplane,"rec_pl");call to_history("d1", rec%h%d,"rec_pl");call to_history("o1", oxplane,"rec_pl")
 call to_history("n2",       1,"rec_pl");call to_history("d2",      1.,"rec_pl");call to_history("o2",      0.,"rec_pl")
 call to_history("n3", rec%w%n,"rec_pl");call to_history("d3", rec%w%d,"rec_pl");call to_history("o3", rec%w%o,"rec_pl")
 call to_history("n4",   ang%n,"rec_pl");call to_history("d4",   ang%d,"rec_pl");call to_history("o4",   ang%o,"rec_pl")  
 call to_history("esize",8,"rec_pl")

 call to_history("n1", nxplane,"sou_pl");call to_history("d1", rec%h%d,"sou_pl");call to_history("o1", oxplane,"sou_pl")
 call to_history("n2",       1,"sou_pl");call to_history("d2",      1.,"sou_pl");call to_history("o2",      0.,"sou_pl")
 call to_history("n3", rec%w%n,"sou_pl");call to_history("d3", rec%w%d,"sou_pl");call to_history("o3", rec%w%o,"sou_pl")
 call to_history("n4",   ang%n,"sou_pl");call to_history("d4",   ang%d,"sou_pl");call to_history("o4",   ang%o,"sou_pl")  
 call to_history("esize",8,"sou_pl")

 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
!--------------------------------------------------
! Frequency vector
!
rec%w%o   = 2*pi * rec%w%o
rec%w%d   = 2*pi * rec%w%d
allocate(rec%w%all(rec%w%n),phase1(rec%w%n))
do iw = 1,rec%w%n
   rec%w%all(iw) = rec%w%o + (iw - 1) * rec%w%d
end do
 call owexec_time(start1)
 if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
 source1=0;
 call sreed("sou",source1,8*sou%w%n)
 if (allocated(source_pl)) deallocate(source_pl); allocate(source_pl(nxplane,sou%w%n,ang%n))
 source_pl=0
 phase1=cmplx(cos(pi/4.),-sin(pi/4))*sqrt(rec%w%all)
 do idip=1,nxplane
    dx=(idip-int(nxplane/2.+1))*rec%h%d
!$OMP parallel do private (iw,iang,alpha,phase) schedule(dynamic,1) num_threads(node)
    do iw = 1,rec%w%n
       do iang=1,ang%n
          alpha=sind((iang-1)*ang%d+ang%o)
          phase=exp(cmplx(0,(dx*rec%w%all(iw)*alpha/v0)))
          source_pl(idip,iw,iang)=source_pl(idip,iw,iang)+source1(iw)*phase*phase1(iw)
       end do
    end do
!$OMP END PARALLEL DO
    if (mod(idip,20) ==1) write(0,*) "Encoded source", idip
 end do

 call srite("sou_pl",source_pl,8*nxplane*sou%w%n*ang%n)
 deallocate(source_pl,source1)
 call owexec_time(finish1)

 write(0,*) "Encoded sources in ",(finish1-start1)/60," min"
 start1=0.;finish1=0.;

 call owexec_time(start1)

 if (allocated(recv_pl  )) deallocate(recv_pl  ); allocate(  recv_pl(nxplane,sou%w%n,ang%n))
 recv_pl=0;
 do idip=1,rec%xd%n
    dx=(idip-int(rec%xd%n/2.+1))*rec%xd%d
    owrk=(idip-1)*rec%xd%d+rec%xd%o+rec%h%o
    iowrk=(owrk-oxplane)/rec%h%d+1.5
    call owexec_time(start)
    if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
    recv=0
    call sreed("rec",recv,8*rec%h%n*rec%w%n)
!$OMP parallel do private (iw,iang,alpha,phase) schedule(dynamic,1) num_threads(node)
    do iw = 1,rec%w%n
       do iang=1,ang%n
          alpha=sind((iang-1)*ang%d+ang%o)
          phase=exp(cmplx(0,(dx*rec%w%all(iw)*alpha/v0)))
          recv_pl(iowrk:iowrk-1+rec%h%n,iw,iang)=recv_pl(iowrk:iowrk-1+rec%h%n,iw,iang)+recv(:,iw)*phase*phase1(iw)
       end do
    end do
!$OMP END PARALLEL DO
    if (mod(idip,20) ==1) write(0,*) "Encoded shot", idip
 end do

 call srite("rec_pl",recv_pl,8*nxplane*sou%w%n*ang%n)
 deallocate(recv_pl,recv)
 call owexec_time(finish1)
 write(0,*) "Encoded receivers in ",(finish1-start1)/60," min"

 call sep_end_prog()

end program 
