module OW_encoding

use sep 
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

contains

subroutine owencoding(encrec,encsou)

real						:: start, finish,a(2),b,c
integer						:: ixs,isp,i,j,stat,igap,n1,rxs,sxs,nsuper
integer						:: ni,n_enc_sht,jgap,ih,ib,ic,nth,indi
real, allocatable					:: velint(:,:),vel(:,:),codeX(:),codeall(:,:,:),vref(:,:,:)
real, allocatable					:: rand1(:),phase(:,:),CHD(:,:),hold(:,:),gol(:,:,:),code(:,:)
real, allocatable					:: coscode(:),sincode(:),randw(:,:),phase3(:,:,:)
complex, allocatable				:: rec1(:,:),rec2(:,:),source(:),tfld(:),source1(:,:)
complex, allocatable				:: 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()

!----------------------------------------------------
! 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))
allocate(vref(nref,2,v%z%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)
if (exist_file("vref")) call sreed("vref",vref,v%z%n*2*nref*4)

!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
if (noencode) then
   allocate(source1(kxm%n,rec%w%n))
   source1=0;source1(isp,:)=source
end if

!--------------------------------------------------
! Encoding shots
!
start=0.;finish=0. 
call owexec_time(start)

if (not(noencode)) then 
   if (random) then
      if (encall) then
         allocate(encrec(nreal,img%xm%n,rec%w%n),encsou(nreal,img%xm%n,rec%w%n))
         allocate(soutemp1(rec%w%n-subf),rec2(rec%h%n,rec%w%n))
         allocate(codeall(rec%w%n,rec%x%n,nreal))
         encrec=0;encsou=0
         call init_random_seed3(codeall)
         do ixs=1,rec%x%n
            rec1=0;soutemp1=0;rec2=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)
            do i=1,nreal
               do ih=1,rec%h%n
                  rec2(ih,:)=rec1(ih,:)*exp(cmplx(0,pi/2*codeall(:,ixs,i)*eps))
               end do
               soutemp1=source*exp(cmplx(0,pi/2*codeall(:,ixs,i)*eps))
               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
         nsuper=int(float(rec%x%n/nshots))
         allocate(encrec(nsuper,img%xm%n,rec%w%n),encsou(nsuper,img%xm%n,rec%w%n))
         allocate(soutemp1(rec%w%n),rec2(rec%h%n,rec%w%n))
         allocate(codeall(rec%w%n,rec%x%n,2))
         encrec=0;encsou=0
         call init_random_seed3(codeall)
         do ixs=1,rec%x%n
            i=int(float((ixs-1)/nshots))+1
            rec1=0;soutemp1=0;rec2=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)
            do ih=1,rec%h%n
               rec2(ih,:)=rec1(ih,:)*exp(cmplx(0,2*pi*codeall(:,ixs,1)*eps))
            end do
            soutemp1=source*exp(cmplx(0,2*pi*codeall(:,ixs,1)*eps))
            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
         if(exist_file("encrec")) then
            call to_history("n1",nsuper,"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",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*nsuper*img%xm%n*(rec%w%n-subf))
        end if

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

   else if(gold) then
      if (encall) then
         allocate(encrec(nreal,img%xm%n,rec%w%n),encsou(nreal,img%xm%n,rec%w%n))
         allocate(soutemp1(rec%w%n),rec2(rec%h%n,rec%w%n))
         allocate(codeall(rec%w%n,rec%x%n,nreal))
         encrec=0;encsou=0
         allocate(hold(48,2*rec%w%n-subf))
         allocate(gol(2*rec%w%n,rec%x%n,nreal))
         allocate(randw(rec%x%n,2*nreal))
         allocate(rand1(2*rec%w%n-subf))
   
         do i=1,48
            rand1=0
            call owmseq(2,9,rand1,2*rec%w%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,rec%x%n
            do iw=1,nreal
               b=(randw(ixs,iw)+1)/2;ib=int(b*48)
               c=(randw(ixs,nreal+iw)+1)/2;ic=int(c*48)
               if (ib==0) ib=1
               if (ic==0) ic=1
               if (b==c) then
                  gol(1:2*rec%w%n-subf,ixs,iw)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
               else
                  gol(1:2*rec%w%n-subf,ixs,iw)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib+iw)
               end if
            end do
         end do
         deallocate(hold);allocate(tfld(2*rec%w%n))

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

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

         do ixs=1,rec%x%n
!            i=int(float((ixs-1)/nshots))+1
!            j=ixs-int(float(ixs-1)/nshots)*nshots
            rec1=0;soutemp1=0;rec2=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)
               do iw=1,nreal
                  do ih=1,rec%h%n
                     rec2(ih,:)=rec1(ih,:)*exp(cmplx(0,codeall(:,ixs,iw)*eps))
                  end do
               soutemp1=source*exp(cmplx(0,codeall(:,ixs,iw)*eps))
               encrec(iw,rxs:rxs+rec%h%n-1,:)=encrec(iw,rxs:rxs+rec%h%n-1,:)+rec2
               encsou(iw,sxs,:)=encsou(iw,sxs,:)+soutemp1
            end do
         end do

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

        if(exist_file("encsou")) then
            call to_history("n1",nreal,"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",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*img%xm%n*(rec%w%n-subf))
         end if
      else
         nsuper=int(float(rec%x%n/nshots))
         allocate(encrec(nsuper,img%xm%n,rec%w%n),encsou(nsuper,img%xm%n,rec%w%n))
         allocate(soutemp1(rec%w%n),rec2(rec%h%n,rec%w%n))
         allocate(codeall(rec%w%n,rec%x%n,2))
         encrec=0;encsou=0
         allocate(hold(48,2*rec%w%n-subf))
         allocate(gol(2*rec%w%n,rec%x%n,2))
         allocate(randw(rec%x%n,2))
         allocate(rand1(2*rec%w%n-subf))
   
         do i=1,48
            rand1=0
            call owmseq(2,9,rand1,2*rec%w%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,rec%x%n
            b=(randw(ixs,1)+1)/2;ib=int(b*48)
            c=(randw(ixs,2)+1)/2;ic=int(c*48)
            if (ib==0) ib=1
            if (ic==0) ic=1
            if (b==c) then
               gol(1:2*rec%w%n-subf,ixs,1)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
            else
               gol(1:2*rec%w%n-subf,ixs,1)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib)
            end if
         end do
         deallocate(hold);allocate(tfld(2*rec%w%n))

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

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

         do ixs=1,rec%x%n
            i=int(float((ixs-1)/nshots))+1
            j=ixs-int(float(ixs-1)/nshots)*nshots
            rec1=0;soutemp1=0;rec2=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)
            do ih=1,rec%h%n
               rec2(ih,:)=rec1(ih,:)*exp(cmplx(0,codeall(:,ixs,1)*eps))
            end do
            soutemp1=source*exp(cmplx(0,codeall(:,ixs,1)*eps))
            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
         if(exist_file("encrec")) then
            call to_history("n1",nsuper,"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",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*nsuper*img%xm%n*(rec%w%n-subf))
        end if

        if(exist_file("encsou")) then
            call to_history("n1",nsuper,"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",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*nsuper*img%xm%n*(rec%w%n-subf))
         end if
      end if
   end if
end if
call owexec_time(finish)
write(0,*) "Finished generating codes in ", finish-start," (s)"
!
!if (exist_file("code").and.encall) 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

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