program OWmig2d_par

use sep 
use OW_pspi
use OW_fftw
use OW_type
use OW_exec_time
use OW_parms
use OW_rndseed
use OW_hadam1
use OW_mseq

implicit none

real						:: start, finish
integer						:: ixs,isp,i,stat,igap,n1,rxs,sxs,n_enc_sht,jgap
real, allocatable					:: velint(:,:),vel(:,:)
real, allocatable					:: rand1(:),phase(:,:),CHD(:,:)
complex, allocatable				:: rec1(:,:),rec2(:,:),source(:) !,GG(:,:,:)!,:)
complex, allocatable				:: encrec(:,:,:),encsou(:,:,:),soutemp1(:),soutemp2(:)

 call sep_init()
 call sep_begin_prog()

!---------------------------------------
! IO
!
 call owparam()

!----------------------------------------------------
! Memory Allocation
!
allocate(rec1(rec%h%n,rec%w%n),rec2(rec%h%n,rec%w%n),source(sou%w%n))
allocate(vel(img%z%n,kxm%n),velint(v%z%n,v%xm%n))

!--------------------------------------------------
! Reading source fctn and velocity
 call sreed("sou",source,sou%w%n*8)
 call sreed("vel",velint,v%z%n*v%xm%n*4)

n1 = floor((img%xm%o-v%xm%o)/v%xm%d)+1
vel=velint(1:img%z%n,n1:n1+kxm%n-1);

deallocate(velint)
!--------------------------------------------------
!
!
write(0,*) "DOWNWARD CONTINUES SOURCES AND RECEIVERS"
write(0,*) "========================================"

open(unit=2,file='lixo')
open(unit=3,file='lixo3')
open(unit=4,file='random.txt')
!--------------------------------------------------
! shot position in the acq.geom.
isp = abs(floor(rec%h%o/rec%h%d)) + 1

!--------------------------------------------------
! Encoding shots
!
!
n_enc_sht=int(rec%x%n/ngap*1.)

start=0.;finish=0. 
call owexec_time(start)
if (random) then
   allocate(phase(2*n_enc_sht,rec%w%n))
   do ixs=1,2*n_enc_sht
      allocate(rand1(rec%w%n))
      call init_random_seed(ixs,rec%w%n,rand1)
      phase(ixs,:)=rand1
      deallocate(rand1)
   end do
else if(mseq) then
   allocate(phase(2*n_enc_sht,rec%w%n))
   do ixs=1,2*n_enc_sht
      allocate(rand1(rec%w%n));rand1=0
      call owmseq1(2,8,rand1,ixs)
      phase(ixs,:)=rand1
      deallocate(rand1)
   end do
else if(golay) then
   allocate(phase(2*n_enc_sht,rec%w%n))
   allocate(CHD(rec%w%n,rec%w%n))
   call owhadamard(CHD)
   do i=1,n_enc_sht
      phase(i,:)=CHD((i-1)*int(rec%w%n/n_enc_sht*0.5)+1,:)
      phase(i+n_enc_sht,:)=CHD((i-1)*int(rec%w%n/n_enc_sht*0.5)+int(rec%w%n/2.)+1,:)
   end do
   deallocate(CHD)
else if(gold) then
   allocate(phase(2*n_enc_sht,rec%w%n))
   do ixs=1,2*n_enc_sht
      allocate(rand1(rec%w%n));rand1=0
      call owmseq1(2,8,rand1,ixs)
      phase(ixs,:)=rand1
      deallocate(rand1)
   end do
   allocate(CHD(2*n_enc_sht,rec%w%n));rand1=0
   do ixs=1,2*n_enc_sht
      if (ixs.le.n_enc_sht) then
         CHD(ixs,:)=ieor(int((phase(ixs,:)-1)/2.),int((phase(ixs+n_enc_sht,:)-1)/2.))
      else
         CHD(ixs,:)=ieor(int((phase(ixs,:)-1)/2.),int((phase(ixs-1,:)-1)/2.))
      end if
   end do   
   phase=CHD
   deallocate(CHD)
end if

call owexec_time(finish)
write(0,*) "Finished generating codes in ", finish-start," (s)"

allocate(encrec(2*ngap,img%xm%n,rec%w%n),encsou(2*ngap,img%xm%n,sou%w%n))
allocate(soutemp1(rec%w%n),soutemp2(rec%w%n))
encrec=0;encsou=0

start=0.;finish=0. 
call owexec_time(start)
do ixs=1,rec%x%n
   rec1=0;soutemp1=0;
   rec2=0;soutemp2=0;
   igap=ixs-int((ixs-1)/ngap)*ngap
   jgap=int((ixs-1)/ngap*1.)+1

   rxs = ceiling(((ixs-1)*rec%x%d+rec%x%o+rec%h%o-img%xm%o)/img%xm%d)
   sxs = ceiling(((ixs-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)

   call sreed("in",rec1,8*rec%h%n*rec%w%n)
   rec2=rec1
      do i=1,rec%h%n
         rec1(i,:)=rec1(i,:)*exp(cmplx(0,2*pi*phase(jgap,:)))
         rec2(i,:)=rec2(i,:)*exp(cmplx(0,2*pi*phase(jgap+n_enc_sht,:)))
         soutemp1=source*exp(cmplx(0,2*pi*phase(jgap,:)))
         soutemp2=source*exp(cmplx(0,2*pi*phase(jgap+n_enc_sht,:)))
      end do
!   else if(gold) then
!   end if
   
   encrec(igap,rxs:rxs+rec%h%n-1,:)=encrec(igap,rxs:rxs+rec%h%n-1,:)+rec1
   encrec(igap+ngap,rxs:rxs+rec%h%n-1,:)=encrec(igap+ngap,rxs:rxs+rec%h%n-1,:)+rec2

   encsou(igap,sxs,:)=encsou(igap,sxs,:)+soutemp1
   encsou(igap+ngap,sxs,:)=encsou(igap+ngap,sxs,:)+soutemp2
end do

call owexec_time(finish)
write(0,*) "Finished encoding shots in ", finish-start," (s)"

if(exist_file("encrec")) then
 call to_history("n1",ngap,"encrec"); call to_history("d1",1.,"encrec"); call to_history("o1",0.,"encrec")  
 call to_history("n2",img%xm%n,"encrec"); call to_history("d2",img%xm%d,"encrec"); call to_history("o2",img%xm%o,"encrec")
 call to_history("n3",rec%w%n,"encrec") ; call to_history("d3",rec%w%d/(2*pi),"encrec") ; call to_history("o3",rec%w%o/(2*pi),"encrec")
 call to_history("esize",8,"encrec")
 call srite("encrec",encrec,8*ngap*img%xm%n*rec%w%n)
end if

if(exist_file("encsou")) then
 call to_history("n1",ngap,"encsou"); call to_history("d1",1.,"encsou"); call to_history("o1",0.,"encsou")  
 call to_history("n2",img%xm%n,"encsou"); call to_history("d2",img%xm%d,"encsou"); call to_history("o2",img%xm%o,"encsou")
 call to_history("n3",rec%w%n,"encsou") ; call to_history("d3",rec%w%d/(2*pi),"encsou") ; call to_history("o3",rec%w%o/(2*pi),"encsou")
 call to_history("esize",8,"encsou")
 call srite("encsou",encsou,8*ngap*img%xm%n*rec%w%n)
end if

deallocate(rec1,rec2,soutemp1,soutemp2,source)

!--------------------------------------------------
! Migrating
!
!-----------------------------------------------------------------------------------------------      
shot:do ixs=1,2*ngap

    start=0.;finish=0. 
    call owexec_time(start)
    call owpspi_init(isp,ixs,wbott)
    call owpspi(encrec(ixs,:,:),encsou(ixs,:,:),vel)
    call owexec_time(finish)
    write(0,*) " DONE SHOT # ",ixs," OUT OF ",2*ngap," in ", (finish - start)/60.," (min)"

end do shot

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

call sep_end_prog()
end program OWmig2d_par
