module goldcode

use fftw
use OW_mseq
use OW_gcd
USE simulated_anneal

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 (x,nhorz,encz)

INTEGER, PARAMETER            :: dp = SELECTED_REAL_KIND(14, 60)
logical			:: encz
real(dp)			:: x(:)
real,allocatable		:: randz(:,:,:,:)
integer			:: nhorz,nodes,whichSeq
integer, external 		:: omp_get_thread_num,omp_get_num_threads

allocate(randz(2*nw-1,nx,nhorz,nreal))
if (nodes.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
end if

if (pwval==2) whichSeq=1
if (pwval==3) whichSeq=2
if (pwval==4) whichSeq=2
if (pwval==5) whichSeq=6
if (pwval==6) whichSeq=6
if (pwval==7) whichSeq=18
if (pwval==8) whichSeq=16
if (pwval==9) whichSeq=48
if (pwval==10) whichSeq=60

   allocate(hold(whichSeq,2*nw-1))
   allocate(rand1(2*nw-1))
   allocate(q(2*nw-1));q=0
   call owgcd(pwval,nq,q,indq)

   do i=1,whichSeq
      rand1=0
      call owmseq(2,pwval,rand1,2*nw,i-int((i-1)/whichSeq)*whichSeq)
      hold(i,:)=rand1
   end do
   rand1=0
   hold=hold*2-1
   allocate(prefer(whichSeq,indq(1),2*nw-1))
   do iq=1,indq(1)
      do i=1,whichSeq
         prefer(i,iq,1)=hold(i,1)
         rand1=hold(i,:)      
         do iz=2,2*nw-1
            rand1=cshift(rand1,q(iq))
            prefer(i,iq,iz)=rand1(1)
         end do
      end do
   end do
   deallocate(rand1)
   allocate(gol(whichSeq,indq(1),2*nw-1,2*nw-1));gol=0;

!$OMP parallel do private (ixs,&
                         ihorz,&
                             i) schedule(dynamic,1) num_threads(nodes)
   do ixs=1,2*nw-1
      do ihorz=1,indq(1)
         do i=1,whichSeq
            gol(i,ihorz,:,ixs)=hold(i,:)*cshift(prefer(i,ihorz,:),ixs)
         end do
      end do
   end do
!$OMP END PARALLEL DO
   deallocate(hold)
   if (allocated(tfld)) deallocate(tfld);allocate(tfld(2*nw-1))
   call owfftw(tfld)
   if (encz) then
!      gap=int((2*nw-1)/(nhorz-1))
      gap=int((2*nw-1)/(nhorz))
   else
      gap=int((2*nw-1)/(nx/ngap))
   end if

   if (encz) then
      do ihorz=1,nhorz
         do ixs=1,nx
!!            jgap=int((ihorz-1)/(nhorz))*gap+1
!!            randz(:,ixs,ihorz,1)=gol(1,ihorz,jgap,:)
            randz(:,ixs,ihorz,1)=gol(1,ihorz,ixs,:)
!            tfld=0
!            tfld=cmplx(gol(1,ihorz,jgap,:),0)
!            call owfftexe(tfld)
!            randz(:,ixs,ihorz,1)=atan(imag(tfld(1:nw))*real(tfld(1:nw))/ &
!                      (real(tfld(1:nw))*real(tfld(1:nw))+tiny(1.)))
         end do
      end do
   else
       do ihorz=1,nhorz
         do ixs=1,nx
!!            jgap=int((ixs-1)/(ngap)+1)*gap
!!            if (jgap > 2*nw-1) jgap=jgap/2
!!            randz(:,ixs,ihorz,1)=gol(1,ihorz,jgap,:)
            randz(:,ixs,ihorz,1)=gol(1,ihorz,ixs,:)
!            tfld=0
!            tfld=cmplx(gol(1,ihorz,jgap,:),0)
!            call owfftexe(tfld)
!            randz(:,ixs,ihorz,1)=atan(imag(tfld(1:nw))*real(tfld(1:nw))/ &
!                      (real(tfld(1:nw))*real(tfld(1:nw))+tiny(1.)))
         end do
      end do
   end if
!   randz=randz*2
   do ixs=1,nx
      x((ixs-1)*(2*nw-1)+1:ixs*(2*nw-1))=dble(randz(:,ixs,nhorz,1))
   end do
   call owfftwdestroy()
   deallocate(tfld,gol,randz)   
write(0,*) "Finished generating Gold codes"
end subroutine
end module
