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_rndseed1
use OW_hadam1
use OW_mseq

implicit none

real						:: start, finish,a(2),b,c
integer						:: ixs,isp,i,j,stat,igap,n1,rxs,sxs,n_enc_sht,jgap,ih,ib,ic,nth
real, allocatable					:: velint(:,:),vel(:,:),codeX(:),codeall(:,:,:)
real, allocatable					:: rand1(:),phase(:,:),CHD(:,:),hold(:,:),gol(:,:,:),code(:,:)
real, allocatable					:: coscode(:),sincode(:),randw(:,:),phase3(:,:,:)
complex, allocatable				:: rec1(:,:),rec2(:,:),source(:),tfld(:),tflds(:),tfldr(:)
complex, allocatable				:: encrec(:,:,:),encsou(:,:,:),soutemp1(:),soutemp2(:)

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

real,allocatable					:: image(:,:,:),GG(:,:,:)

 call sep_init()
 call sep_begin_prog()

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

!--------------------------------------------------
! Setting FFTW plans
if (node.eq.9999) then
   nth = omp_get_num_procs()
end if
!write(0,*) "**************",nth,"*****************"
!allocate(tflds(kxm%n));tflds=0
!call owfftws(tflds,nth)
!deallocate(tflds)
!
!allocate(tfldr(kxm%n));tfldr=0
!call owfftwr(tfldr,nth)
!deallocate(tfldr)

!----------------------------------------------------
! Memory Allocation
!
allocate(rec1(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.)

write(0,*) "subf ",subf

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

if (noencode) then 
   allocate(code(rec%x%n,sou%w%n-subf))
   code=0.
else
   if (random) then
     if (randall) then
        allocate(codeall(rec%w%n,rec%x%n,nreal))
        call init_random_seed3(codeall)
        allocate(phase3(rec%w%n,rec%x%n,nreal))
        phase3=cos(2*pi*codeall)
!########
!codeall=codeall*0.1
     else
        allocate(phase(nreal*n_enc_sht,rec%w%n-subf))
        allocate(code(nreal*n_enc_sht,rec%w%n-subf))
         do ixs=1,nreal*n_enc_sht
            allocate(rand1(rec%w%n-subf))
            call init_random_seed(ixs,rec%w%n-subf,rand1)
            code(ixs,:)=rand1
            deallocate(rand1)
         end do
!         phase=cos(2*pi*code)
     end if
!   else if(mseq) then
!
!      allocate(code(nreal*n_enc_sht,rec%w%n-subf))
!      allocate(phase(nreal*n_enc_sht,rec%w%n-subf))
!      allocate(tfld(rec%w%n-subf))
!      do ixs=1,nreal*n_enc_sht
!         allocate(rand1(rec%w%n-subf));rand1=0
!         call owmseq(2,8,rand1,2*(rec%w%n),ixs-int((ixs-1)/16)*16)
!         tfld=0.;tfld=cmplx(rand1,0)
!         rand1=2*rand1-1
!         phase(ixs,:)=rand1
!         call FFT1DT(tfld, 1)
!         code(ixs,:)=atan(imag(tfld)*real(tfld)/(real(tfld)*real(tfld)+tiny(1.)))
!         deallocate(rand1)
!      end do
!      deallocate(tfld)

   !else if(golay) then
   !  allocate(code(nreal*n_enc_sht,rec%w%n))
   !  allocate(phase(nreal*n_enc_sht,rec%w%n))
   !   allocate(CHD(rec%w%n,rec%w%n))
   !   call owhadamard(CHD)
   !   do i=1,nreal*n_enc_sht
   !      phase(i,:)=CHD((i-1)*int(rec%w%n/(nreal*n_enc_sht)*1.)+1,:)
   !      phase(i+n_enc_sht,:)=CHD((i-1)*int(rec%w%n/(nreal*n_enc_sht)*1.)+int(rec%w%n/2.)+1,:)
   !   end do
   !   deallocate(CHD)
   else if(gold) then
   
      allocate(codeall(rec%w%n,rec%x%n,nreal))
      allocate(hold(16,rec%w%n-subf))
      allocate(gol(rec%w%n,rec%x%n,nreal))
      allocate(phase(nreal*rec%x%n,rec%w%n))
      allocate(randw(rec%x%n,2))
      allocate(rand1(rec%w%n-subf))

      do i=1,16
         rand1=0
         call owmseq(2,8,rand1,2*(rec%w%n-subf),i-int((i-1)/16)*16)
         hold(i,:)=rand1
      end do
      deallocate(rand1)

      hold=hold*2-1

   do i=1,nreal
      call init_random_seed2(randw)
      do ixs=1,rec%x%n
         b=(randw(ixs,1)+1)/2;ib=int(b*16)
         c=(randw(ixs,2)+1)/2;ic=int(c*16)
         if (ib==0) ib=1
         if (ic==0) ic=1
         write(2,*) "nreal = ",i," tiro = ",ixs,"ib = ",ib,"  ic = ",ic
         do iw=1,rec%w%n
            if (b==c) then
               gol(1:rec%w%n-subf,ixs,i)=hold(ib,:)*cshift(hold(ic+1,:),iw+shift)
!               if (ixs==1) write(2,*) "*******"; write(2,*) gol
            else
               gol(1:rec%w%n-subf,ixs,i)=hold(ib,:)*cshift(hold(ic,:),iw+shift+ib)
!               if (ixs==1) write(2,*) "*******"; write(2,*) gol
            end if
         end do
      end do
   end do
   deallocate(hold);allocate(tfld(rec%w%n))

   do i=1,nreal
      do ixs=1,rec%x%n
         if (gol(1,ixs,i).eq.1) then
            gol(rec%w%n,ixs,i)=-1
         else
            gol(1,ixs,i)=1;gol(rec%w%n,ixs,i)=1
         end if
      end do
   end do

   do i=1,nreal
      do j=1,rec%x%n
         tfld=0
         tfld=cmplx(gol(:,j,i),0)
         call FFT1DT(tfld, 1)
         codeall(:,j,i)=atan(imag(tfld)*real(tfld)/(real(tfld)*real(tfld)+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   

   allocate(phase3(rec%w%n,rec%x%n,nreal))
   phase3=cos(2*pi*codeall)

   else if(xiao) then
      allocate(codeX(nreal*rec%x%n));codeX=0
      allocate(coscode(n_enc_sht),sincode(n_enc_sht))
      do ixs=1,ngap*nreal
         call init_random_seed(ixs,size(a),a)
         codeX((ixs-1)*n_enc_sht+1)=a(1)
         coscode=0;sincode=0;coscode(1)=cos(2*pi*a(1));sincode(1)=sin(2*pi*a(1))
         do i=2,n_enc_sht
            coscode(i)=cos(2*pi*codex((ixs-1)*n_enc_sht+i-1))
            sincode(i)=sin(2*pi*codex((ixs-1)*n_enc_sht+i-1))
!            codeX((ixs-1)*n_enc_sht+i)=atan(-(sum(coscode*sincode/(sincode**2+tiny(1.)))))
            codeX((ixs-1)*n_enc_sht+i)=atan(-sum(coscode)/(sum(sincode)+tiny(1.)))
         end do
      end do
   end if
end if
call owexec_time(finish)
write(0,*) "Finished generating codes in ", finish-start," (s)"

if (exist_file("code").and.randall) then
 call to_history("n1",rec%w%n,"code"); call to_history("d1",rec%w%d/(2*pi),"code"); call to_history("o1",rec%w%o/(2*pi),"code")  
 call to_history("n2",rec%x%n,"code"); call to_history("d2",rec%x%d,"code"); call to_history("o2",rec%x%o,"code")
 call to_history("n3",nreal,"code") ; call to_history("d3",1.,"code") ; call to_history("o3",1.,"code")
 call to_history("esize",4,"code")
 call srite("code",phase3,4*nreal*rec%x%n*rec%w%n)
 goto 10000
end if


start=0.;finish=0. 
call owexec_time(start)
if (xiao) then
   allocate(encrec(nreal*ngap,img%xm%n,rec%w%n),encsou(nreal*ngap,img%xm%n,rec%w%n))
   allocate(soutemp1(rec%w%n))
   encrec=0;encsou=0
   do ixs=1,ngap
      do i=1,n_enc_sht
         rxs = ceiling((((ixs-1)*n_enc_sht+i-1)*rec%x%d+rec%x%o+rec%h%o-img%xm%o)/img%xm%d)
         sxs = ceiling((((ixs-1)*n_enc_sht+i-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)
         write(0,*) ixs,i,rxs,rxs+rec%h%n-1,sxs
         rec1=0;soutemp1=0;
         call sreed("in",rec1,8*rec%h%n*rec%w%n)
         rec1=rec1*exp(cmplx(0,2*pi*codeX((ixs-1)*n_enc_sht+i)))
         soutemp1=source*exp(cmplx(0,2*pi*codeX((ixs-1)*n_enc_sht+i)))
         encrec(ixs,rxs:rxs+rec%h%n-1,:)=encrec(ixs,rxs:rxs+rec%h%n-1,:)+rec1
         encsou(ixs,sxs,:)=encsou(ixs,sxs,:)+soutemp1
      end do
   end do
else
   if ((random.and.randall).or.(gold.and.randall)) then
      allocate(encrec(nreal,img%xm%n,rec%w%n-subf),encsou(nreal,img%xm%n,rec%w%n-subf))
      allocate(soutemp1(rec%w%n-subf),rec2(rec%h%n,rec%w%n))! ,soutemp2(rec%w%n-subf))
      encrec=0;encsou=0
      do ixs=1,rec%x%n
         rec1=0;soutemp1=0;
         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)
         write(0,*) "Encoded shot # ",ixs
         do i=1,nreal
            do ih=1,rec%h%n
               rec2(ih,:)=rec1(ih,:)*exp(cmplx(0,2*pi*codeall(:,ixs,i)))
            end do
            soutemp1=source*exp(cmplx(0,2*pi*codeall(:,ixs,i)))
            encrec(i,rxs:rxs+rec%h%n-1,:)=encrec(i,rxs:rxs+rec%h%n-1,:)+rec2
            encsou(i,sxs,:)=encsou(i,sxs,:)+soutemp1
         end do
      end do
   else
      allocate(encrec(nreal*ngap,img%xm%n,rec%w%n-subf),encsou(nreal*ngap,img%xm%n,rec%w%n-subf))
      allocate(soutemp1(rec%w%n-subf))! ,soutemp2(rec%w%n-subf))
      encrec=0;encsou=0
      do ixs=1,rec%x%n
         rec1=0;soutemp1=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)
         do i=1,rec%h%n
            rec1(i,1:rec%w%n-subf)=rec1(i,1:rec%w%n-subf)*exp(cmplx(0,2*pi*code(jgap,1:rec%w%n-subf)))
         end do
         soutemp1=source(1:rec%w%n-subf)*exp(cmplx(0,2*pi*code(jgap,1:rec%w%n-subf)))
         encrec(igap,rxs:rxs+rec%h%n-1,:)=encrec(igap,rxs:rxs+rec%h%n-1,:)+rec1(:,1:rec%w%n-subf)
         encsou(igap,sxs,:)=encsou(igap,sxs,:)+soutemp1
      end do
   end if
end if
call owexec_time(finish)
write(0,*) "Finished encoding shots in ", finish-start," (s)"

if(exist_file("encrec")) then
 call to_history("n1",nreal*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-subf,"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*nreal*ngap*img%xm%n*(rec%w%n-subf))
end if

if(exist_file("encsou")) then
 call to_history("n1",nreal*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-subf,"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*nreal*ngap*img%xm%n*(rec%w%n-subf))
end if

deallocate(rec1,soutemp1,source)

!--------------------------------------------------
! Migrating
!
!-----------------------------------------------------------------------------------------------      
if (randall) then
   do ixs=1,nreal
       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 ",nreal," in ", (finish - start)/60.," (min)"
   end do
else
   do ixs=1,nreal*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 ",nreal*ngap," in ", (finish - start)/60.," (min)"
   end do
end if
close(2)
close(3)
close(4)

10000 call sep_end_prog()
end program OWmig2d_par