module OW_encoding_areal

use sep 
use OW_type
use OW_parms
use OW_fftw
use OW_rndseed1
use OW_mseq
use OW_gcd
use OW_exec_time

implicit none

real,private			:: start,finish
integer,private			:: i,ixs,iw,iz,iq,nq
integer,private			:: igap,jgap,indq(1),ihorz,gap,irand
integer,allocatable,private		:: q(:),prefer(:,:,:)
real,allocatable,private		:: rand1(:),hold(:,:),gol(:,:,:,:),correl(:,:)
complex,allocatable,private		:: tfld(:),fun1(:)

contains 

!==================================================
subroutine owgoldpw (source,encrec,encsou)
real			:: randz(rec%w%n,rec%x%n)
integer			:: nodes
integer, external 		:: omp_get_thread_num,omp_get_num_threads
!-------------------------
integer			:: rxs,sxs,ih
complex			:: rec1(rec%h%n,rec%w%n),soutemp1(rec%w%n),rec2(rec%h%n,rec%w%n)
complex			:: encrec(:,:,:),encsou(:,:,:),source(:)

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

   allocate(hold(48,2*sou%w%n-subf))
   allocate(rand1(2*sou%w%n-subf))
   allocate(q(2*sou%w%n-subf));q=0
   call owexec_time(start)
   call owgcd(9,nq,q,indq)
   call owexec_time(finish)

   write(0,*) "gcd ",(finish-start)/60,"(min)"

   call owexec_time(start)
   do i=1,48
      rand1=0
      call owmseq(2,9,rand1,2*sou%w%n,i-int((i-1)/48)*48)
      hold(i,:)=rand1
   end do
   rand1=0
   hold=hold*2-1
   call owexec_time(finish)

   write(0,*) "mseq ",(finish-start)/60,"(min)"

   allocate(prefer(48,indq(1),2*sou%w%n-subf))

   call owexec_time(start)
   do iq=1,indq(1)
      do i=1,48
         prefer(i,iq,1)=hold(i,1)
         rand1=hold(i,:)      
         do iz=2,2*sou%w%n-subf
            rand1=cshift(rand1,q(iq))
            prefer(i,iq,iz)=rand1(1)
         end do
      end do
   end do
   call owexec_time(finish)

   write(0,*) "prefer ",(finish-start)/60,"(min)"

   deallocate(rand1)
   allocate(gol(48,indq(1),2*sou%w%n-subf,2*sou%w%n-subf));gol=0;

!
! ii can not be greater than 48
!
   call owexec_time(start)
!$OMP parallel do private (ixs,&
                         ihorz,&
                             i) schedule(dynamic,1) num_threads(nodes)
   do ixs=1,2*sou%w%n-subf
      do ihorz=1,indq(1)
         do i=1,48
            gol(i,ihorz,:,ixs)=hold(i,:)*cshift(prefer(i,ihorz,:),ixs)
         end do
      end do
   end do
!$OMP END PARALLEL DO

   call owexec_time(finish)
   write(0,*) "gold ",(finish-start)/60,"(min)"

   deallocate(hold)

   if (allocated(tfld)) deallocate(tfld);allocate(tfld(2*sou%w%n-subf))

   call owfftw(tfld)

   call owexec_time(start)
   gap=int((2*sou%w%n-subf)/nshots)
   do ixs=1,rec%x%n
      jgap=int((ixs-1)/(rec%x%n/nshots)+1.)*gap
!      jgap=int((ixs-1)/(rec%x%n/nshots)+1.)*gap
!      jgap=(ixs-1-int((ixs-1)/(rec%x%n/nshots))*(rec%x%n/nshots))*gap+1

      if (jgap > 2*sou%w%n-subf) jgap=jgap/2
      tfld=0
      tfld=cmplx(gol(1,1,jgap,:),0)
      call sfftw_execute(planF)
      randz(:,ixs)=atan(imag(tfld(1:sou%w%n))*real(tfld(1:sou%w%n))/ &
                (real(tfld(1:sou%w%n))*real(tfld(1:sou%w%n))+tiny(1.)))
   end do
   randz=randz*2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   rec1=0;
   do ixs=1,rec%x%n
      i=ixs-int(float((ixs-1)/(rec%x%n/nshots)))*(rec%x%n/nshots)
      soutemp1=0;rec2=0.
      rxs = (((ixs-1)*rec%x%d+rec%x%o+rec%h%o-img%xm%o)/img%xm%d)+1.5
      sxs = (((ixs-1)*rec%x%d+rec%x%o-img%xm%o)/img%xm%d)+1.5
      write(0,*) ixs, i, rxs, sxs
      call sreed("in",rec1,8*rec%h%n*rec%w%n)
! !$OMP parallel do private (iw,ih) schedule(dynamic,1) num_threads(nodes)
      do iw=1,rec%w%n
         do ih=1,rec%h%n
            rec2(ih,iw)=rec1(ih,iw)*exp(cmplx(0,randz(iw,ixs)*eps))
            encrec(rxs+ih-1,iw,i)=encrec(rxs+ih-1,iw,i)+rec2(ih,iw)
         end do
         soutemp1(iw)=source(iw)*exp(cmplx(0,randz(iw,ixs)*eps))
         encsou(sxs,iw,i)=encsou(sxs,iw,i)+soutemp1(iw)
      end do
! !$OMP END PARALLEL DO
      if (mod(ixs,20).eq.1) write(0,*) "encoded shot ",ixs
   end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   call owexec_time(finish)
   write(0,*) "encode ",(finish-start)/60,"(min)"

   call owfftdestroy1()

!if (exist_file("encrec")) then
! call to_history("n1",size(encrec,1),"encrec"); call to_history("d1",1.,"encrec"); call to_history("o1",1.,"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",1,"encrec"); call to_history("d3",1.,"encrec"); call to_history("o3",0.,"encrec")
! call to_history("n4",1,"encrec"); call to_history("d4",1.,"encrec"); call to_history("o4",0.,"encrec")
! call to_history("n5",rec%w%n,"encrec") ; call to_history("d5",rec%w%d/(2*pi),"encrec") ; call to_history("o5",rec%w%o/!(2*pi),"encrec")
! call to_history("esize",8,"encrec")
! call srite("encrec",encrec,8*size(encrec,1)*img%xm%n*rec%w%n)
!end if
!
!if (exist_file("encsou")) then
! call to_history("n1",size(encsou,1),"encsou"); call to_history("d1",1.,"encsou"); call to_history("o1",1.,"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",1,"encsou"); call to_history("d3",1.,"encsou"); call to_history("o3",0.,"encsou")
! call to_history("n4",1,"encsou"); call to_history("d4",1.,"encsou"); call to_history("o4",0.,"encsou")
! call to_history("n5",rec%w%n,"encsou") ; call to_history("d5",rec%w%d/(2*pi),"encsou") ; call to_history("o5",rec%w%o/!(2*pi),"encsou")
! call to_history("esize",8,"encsou")
! call srite("encsou",encsou,8*size(encsou,1)*img%xm%n*rec%w%n)
!end if

   deallocate(tfld,gol)   

end subroutine
!========================================================

end module
!      if (mod(ixs,nshots).eq.0) then
!         j=int(rec%x%n/nshots)*(nshots-1)+int((ixs-1)/nshots+1)
!      else    
!         j=int(rec%x%n/nshots)*(mod(ixs,nshots)-1)+int((ixs-1)/nshots+1)
!      end if
