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

implicit none

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(:,:,:)
real, allocatable					:: rand1(:),phase(:,:),CHD(:,:),hold(:,:),gol(:,:,:),code(:,:)
real, allocatable					:: coscode(:),sincode(:),randw(:,:),phase3(:,:,:)
complex, allocatable				:: rec1(:,:),rec2(:,:),source(:),tfld(:),source1(:,:)
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()

!----------------------------------------------------
! 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)

!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-subf),encsou(nreal,img%xm%n,rec%w%n-subf))
         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
            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)
!write(0,*) ixs,i,j
            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
!         write(0,*) "Encoded shot #", ixs
         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
      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*16)
            c=(randw(ixs,2)+1)/2;ic=int(c*16)
            if (ib==0) ib=1
            if (ic==0) ic=1
!            do iw=1,rec%w%n
               if (b==c) then
                  gol(1:2*rec%w%n-subf,ixs,1)=hold(ib,:)*cshift(hold(ic+1,:),ixs+shift)
   !!               if (ixs==1) write(2,*) "*******"; write(2,*) gol
               else
                  gol(1:2*rec%w%n-subf,ixs,1)=hold(ib,:)*cshift(hold(ic,:),ixs+shift+ib)
   !!               if (ixs==1) write(2,*) "*******"; write(2,*) gol
               end if
!            end do
         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)
!write(0,*) ixs,i,j
            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
!         write(0,*) "Encoded shot #", ixs
         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)"

!--------------------------------------------------
! Migrating
!
!-----------------------------------------------------------------------------------------------  
allocate(GG(kxm%n,img%sh%n,img%z%n))    
if (noencode) then
   allocate(image(img%xm%n,img%sh%n,img%z%n))
   do ixs=1,rec%x%n
       start=0.;finish=0.;GG=0;image=0;
       call owexec_time(start)
       call sreed("in",rec1,8*rec%h%n*rec%w%n)
       call owgrab_vel(ixs,velint,vel,ni)
       call owpspi_init(isp,ixs,wbott)
       call owpspi(rec1,source1,vel,GG)
       call owexec_time(finish)
       image(ni:ni+kxm%n-1,:,:)=GG
       call owimage_io(image,ixs)
       write(0,*) " DONE SHOT # ",ixs," OUT OF ",rec%x%n," in ", (finish - start)/60.," (min)"
   end do
!   call srite("image",image,4*img%z%n*img%sh%n*img%xm%n)
   deallocate(image,GG)
else if ((random.and.encall).or.(gold.and.encall)) then
   n1 = floor((img%xm%o-v%xm%o)/v%xm%d)+1
   vel=velint(1:img%z%n,n1:n1+kxm%n-1);
   do ixs=1,nreal
       start=0.;finish=0.;GG=0
       call owexec_time(start)
       call owpspi_init(isp,ixs,wbott)
       call owpspi(encrec(ixs,:,:),encsou(ixs,:,:),vel,GG)
       call owexec_time(finish)
       call owimage_io(GG,ixs)
       write(0,*) " DONE SHOT # ",ixs," OUT OF ",nreal," in ", (finish - start)/60.," (min)"
   end do
   deallocate(GG)
else if ((random.and.not(encall)).or.(gold.and.not(encall))) then
   n1 = floor((img%xm%o-v%xm%o)/v%xm%d)+1
   vel=velint(1:img%z%n,n1:n1+kxm%n-1);

   do ixs=1,nsuper
       start=0.;finish=0.;GG=0
       call owexec_time(start)
       call owpspi_init(isp,ixs,wbott)
       call owpspi(encrec(ixs,:,:),encsou(ixs,:,:),vel,GG)
       call owexec_time(finish)
       call owimage_io(GG,ixs)
       write(0,*) " DONE SHOT # ",ixs," OUT OF ",nsuper," in ", (finish - start)/60.," (min)"
   end do
   deallocate(GG)
end if
close(2)
close(3)
close(4)

10000 call sep_end_prog()
end program OWmig2d_par
