module OW_encoding

use sep 
use OW_type
use OW_parms_areal
use OW_fftw_areal
use OW_rndseed1
use OW_mseq
use OW_gcd

implicit none

real,private			:: a(2),b,c
integer,private			:: i,ixs,iw,iz,iq,nq,j
integer,private			:: igap,jgap,ib,ic,ireal,indq(1),ihorz,iih
integer,allocatable,private		:: q(:),prefer(:,:,:)
real,allocatable,private		:: randw(:,:),rand4(:,:,:,:),rand1(:),randx(:)
real,allocatable,private		:: hold(:,:),gol(:,:,:,:),randq(:,:),coscode(:,:),sincode(:,:)
complex,allocatable,private		:: tfld(:)
contains 

!====================================================
subroutine owgoldz (randz)
real			:: randz(:,:,:)

   allocate(hold(48,2*sog%z%n-subf))
   allocate(gol(2*sog%z%n,sou%w%n,sog%xm%n))
   allocate(randw(sog%xm%n*sou%w%n,2))
   allocate(rand1(2*sog%z%n-subf))
   do i=1,48
      rand1=0
      call owmseq(2,9,rand1,2*sog%z%n,i-int((i-1)/48)*48)
      hold(i,:)=rand1
   end do
   deallocate(rand1)
   hold=hold*2-1
   call init_random_seed2(randw)
   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         b=(randw((ixs-1)*sog%xm%n+iw,1)+1)/2;ib=int(b*49)
         c=(randw((ixs-1)*sog%xm%n+iw,2)+1)/2;ic=int(c*49)
         if (ib==0) ib=1
         if (ic==0) ic=1
         if (b==c) then
            gol(1:2*sog%z%n-subf,iw,ixs)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
         else
            gol(1:2*sog%z%n-subf,iw,ixs)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib)
         end if
      end do
   end do
   deallocate(hold)

   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         if (gol(1,iw,ixs).eq.1) then
            gol(2*sog%z%n,iw,ixs)=-1
         else
            gol(1,iw,ixs)=1;gol(2*sog%z%n,iw,ixs)=1
         end if
      end do
   end do

   allocate(tfld(2*sog%z%n))

   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         tfld=0
         tfld=cmplx(gol(:,iw,ixs),0)
         call sfftw_execute(planF)
         randz(:,iw,ixs)=atan(imag(tfld(sog%z%n+1:2*sog%z%n))*real(tfld(sog%z%n+1:2*sog%z%n))/ &
                            (real(tfld(sog%z%n+1:2*sog%z%n))*real(tfld(sog%z%n+1:2*sog%z%n))+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   

end subroutine

!=======================================================
subroutine owgoldpz (randz)
real			:: randz(:,:,:)

   allocate(hold(48,2*sog%z%n-subf))
   allocate(gol(2*sog%z%n,sou%w%n,sog%xm%n))
   allocate(randw(sog%xm%n*sou%w%n,2))
   allocate(rand1(2*sog%z%n-subf))
   allocate(q(2*sog%z%n-subf))

   call owgcd(9,nq,q,indq)

   do i=1,48
      rand1=0
      call owmseq(2,9,rand1,2*sog%z%n,i-int((i-1)/48)*48)
      hold(i,:)=rand1
   end do
   rand1=0
   hold=hold*2-1

   allocate(prefer(48,indq(1),2*sog%z%n-subf))

   do iq=1,indq(1)
      do i=1,48
         prefer(i,iq,1)=hold(i,1)
         rand1=hold(i,:)      
         do iz=2,2*sog%z%n-subf
            rand1=cshift(rand1,-q(iq))
            prefer(i,iq,iz)=rand1(1)
         end do
      end do
   end do
      
   deallocate(rand1)

   call init_random_seed2(randw)

   do ixs=1,sog%xm%n
      igap=ixs-int((ixs-1)/ngap)*ngap
      do iw=1,sou%w%n
         b=(randw((ixs-1)*sog%xm%n+iw,1)+1)/2;ib=int(b*nq)
         if (ib==0) ib=1
         c=(randw((ixs-1)*sog%xm%n+iw,2)+1)/2;ic=int(c*nq)
         if (ib==0) ib=1
         if (ic==0) ic=1
         if (ib==ic) then
            gol(1:2*sog%z%n-subf,iw,ixs)=prefer(igap,ic,:)*cshift(prefer(igap,ib,:),ixs)
         else
            gol(1:2*sog%z%n-subf,iw,ixs)=prefer(igap,ic,:)*cshift(prefer(igap,ib,:),ixs+ib)
         end if
      end do
   end do
   deallocate(hold)

   allocate(tfld(2*sog%z%n-subf))

   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,sou%w%n
         tfld=0
         tfld=cmplx(gol(:,iw,ixs),0)
         call sfftw_execute(planF)
         randz(:,iw,ixs)=atan(imag(tfld(1:sog%z%n))*real(tfld(1:sog%z%n))/ &
                            (real(tfld(1:sog%z%n))*real(tfld(1:sog%z%n))+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   

end subroutine
!========================================================
subroutine owgoldw (randz)
real			:: randz(:,:,:)

   allocate(hold(48,2*sou%w%n-subf))
   allocate(gol(2,2*sou%w%n,sog%xm%n))
   allocate(randw(sog%xm%n*sou%w%n,2))
   allocate(rand1(2*sou%w%n-subf))
   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

   hold=hold*2-1

   allocate(prefer(48,nq,2*sou%w%n-subf))
   rand1=0

   do iq=1,nq
      do i=1,48
         prefer(i,iq,1)=hold(i,1)
         rand1=hold(i,:)      
         do iz=2,2*sog%z%n-subf
            rand1=cshift(rand1,-q(iq))
            prefer(i,iq,iz)=rand1(1)
         end do
      end do
   end do
      
   deallocate(rand1)

   call init_random_seed2(randw)

   do ixs=1,sog%xm%n
      do iw=1,2
         b=(randw((ixs-1)*sog%xm%n+iw,1)+1)/2;ib=int(b*48)
         c=(randw((ixs-1)*sog%xm%n+iw,2)+1)/2;ic=int(c*48)
         if (ib==0) ib=1
         if (ic==0) ic=1
         if (b==c) then
            gol(iw,1:2*sou%w%n-subf,ixs)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
         else
            gol(iw,1:2*sou%w%n-subf,ixs)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib)
         end if
      end do
   end do
   deallocate(hold)

   do ixs=1,sog%xm%n
      do iw=1,2
         if (gol(iw,1,ixs).eq.1) then
            gol(iw,2*sou%w%n,ixs)=-1
         else
            gol(iw,1,ixs)=1;gol(iw,2*sou%w%n,ixs)=1
         end if
      end do
   end do

   allocate(tfld(2*sou%w%n))

   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,2
         tfld=0
         tfld=cmplx(gol(iw,:,ixs),0)
         call sfftw_execute(planF)
         randz(iw,:,ixs)=atan(imag(tfld(sou%w%n+1:2*sou%w%n))*real(tfld(sou%w%n+1:2*sou%w%n))/ &
                            (real(tfld(sou%w%n+1:2*sou%w%n))*real(tfld(sou%w%n+1:2*sou%w%n))+tiny(1.)))
      end do
   end do
   deallocate(tfld,gol)   

end subroutine
!==================================================
subroutine owgoldpw (randz)
real			:: randz(:,:,:,:)

   allocate(hold(48,2*sou%w%n-subf))
   allocate(gol(2*sou%w%n-subf,sog%xm%n,nhorz,nreal));gol=0;
   allocate(rand1(2*sou%w%n-subf))
   allocate(q(2*sou%w%n-subf));q=0

   call owgcd(9,nq,q,indq)

   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

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

   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

   deallocate(rand1)

!
! ii can not be greater than 48
!
   do i=1,nreal
      do ihorz=1,nhorz
         iih=(i-1)*nhorz+ihorz
         do ixs=1,sog%xm%n
            jgap=ngap*(ixs-1)+1+int((ixs-1)/ngap)
            iq=ixs-int((ixs-1)/indq(1))*indq(1)
            gol(:,ixs,ihorz,i)=prefer(ii,1,:)*cshift(prefer(ii,iq,:),jgap)
         end do
      end do
   end do

   do ixs=1,sog%xm%n
      do i=1,2*sou%w%n-subf
         write(3,*) i,ixs, gol(i,ixs,1,1)
      end do
   end do

   deallocate(hold)

   allocate(tfld(2*sou%w%n-subf))

   call owfftw(tfld)

   do i=1,nreal
      do ihorz=1,nhorz
         do ixs=1,sog%xm%n
            tfld=0
            tfld=cmplx(gol(:,ixs,ihorz,i),0)
            call sfftw_execute(planF)
            randz(:,ixs,ihorz,i)=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
      end do
   end do
   randz=randz*2
   deallocate(tfld,gol)   

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

end module
