program OW_arealsht

use sep 
use OW_ctrillwei
use OW_type
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_fftw_areal
use OW_picksh
use OW_rndseed1
use OW_mseq
use OW_gcd

implicit none

real						:: start, finish,start1, finish1,a(2),b,c
real						:: rw(2)=0.,rz(2,2)=0.
integer						:: i,nv,ixs,j,ni,nsy,iw,ih,icomb,ix,stat,iz,ihor,n_enc_sog,iq,nq
integer						:: ind_trac,i_pos,i_ini,i_fin,igap,jgap,jxs,nnn,hij,k,ib,ic
integer,allocatable					:: window(:,:),HD(:,:),q(:),prefer(:,:,:)
real, allocatable					:: velint(:,:),taper(:),vel(:,:),rectmp(:,:),soutmp(:,:)
real,allocatable					:: reccig(:,:),soucig(:,:),soucig_rot(:,:),source1(:),temp(:)
real,allocatable					:: tempr(:,:),temps(:,:),randz(:,:,:),randw(:,:),rand1(:),vref(:,:,:)
real,allocatable					:: hold(:,:),gol(:,:,:),randq(:,:),coscode(:,:),sincode(:,:)
complex,allocatable					:: source(:),GGR(:,:),GGS(:,:)
complex, allocatable				:: GGG(:,:,:),wave_R(:,:,:),wave_S(:,:,:),tfld(:)

integer, external 					:: omp_get_thread_num,omp_get_num_threads,omp_get_num_procs,omp_set_num_threads

 call sep_init()
 call sep_begin_prog()

open(unit=2,file="logfile.txt")
open(unit=3,file="logfile1.txt")
open(unit=4,file="logfile2.txt")
open(unit=5,file="logfile3.txt")

start=0.;finish=0.
call owexec_time(start)

!---------------------------------------
! IO
!
 call owparam(.false.)

!--------------------------------------------------

allocate(source(sou%w%n),source1(sou%w%n))
call sreed("sou",source,sou%w%n*8)
source1=cabs(source)
deallocate(source)

!---------------------------------------------------
!
! Read velocities
!
allocate(velint(v%z%n,v%xm%n))
call sreed("vel",velint,v%z%n*v%xm%n*4)
if (exist_file("vref")) then
   allocate(vref(nref,2,v%z%n))
   call sreed("vref",vref,v%z%n*2*nref*4)
!   do i=1,v%z%n
!      do j=1,2
!     write(0,*) int(vref(:,j,i))
!       end do
!   end do
end if

nv=int((xmin-v%xm%o)/v%xm%d)+1.5
allocate(vel(sog%z%n,nx))
vel(1:sog%z%n,:)=velint(1:sog%z%n,nv:nv+nx-1)
deallocate(velint)

!--------------------------------------------------
! Picks
!
if (exist_file("hor1")) then
   call owpick()
end if
!--------------------------------------------------
!
write(0,*) "GENERATING AREAL DATA"
write(0,*) "====================="

allocate(taper(kxm%n))
call owtaper(int(padtraces/4),taper)

!--------------------------------------------------
!
! Read rotated CIGs
!

if (allocated(reccig)) deallocate(reccig);allocate(reccig(sog%z%n,nx));reccig=0
if (allocated(soucig)) deallocate(soucig);allocate(soucig(sog%z%n,nx));soucig=0

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

if (encode) then
if(rand.and.encz) then 
   allocate(randz(sog%z%n,sou%w%n,sog%xm%n))
   call init_random_seed3(randz)
else if (rand.and.encw) then
   allocate(randw(sou%w%n,sog%xm%n))
   call init_random_seed2(randw)
else if(xiao) then
   n_enc_sog=int(sog%xm%n/ngap)
   allocate(randw(sou%w%n,sog%xm%n));randw=0
   allocate(rand1(sou%w%n));rand1=0
!   allocate(coscode(sou%w%n,sog%xm%n),sincode(sou%w%n,sog%xm%n))
   coscode=0;sincode=0

!   call init_random_seed1(rand1)
!   coscode(:,1)=cos(2*pi*rand1);sincode(:,1)=sin(2*pi*rand1)
!   randw(:,1)=rand1
!
!   do ixs=2,sog%xm%n
!      coscode(:,ixs)=cos(2*pi*coscode(:,ixs-1))
!      sincode(:,ixs)=sin(2*pi*sincode(:,ixs-1))
!!            codeX((ixs-1)*n_enc_sht+i)=atan(-(sum(coscode*sincode/(sincode**2+tiny(1.)))))
!      do iz=1,sou%w%n
!         randw(iz,ixs)=atan(-sum(coscode(iz,1:ixs))/(sum(sincode(iz,1:ixs))+tiny(1.)))
!      end do
!   end do

   allocate(coscode(sou%w%n,n_enc_sog),sincode(sou%w%n,n_enc_sog))
   do ixs=1,ngap
      call init_random_seed1(rand1,ixs)
      coscode=0;sincode=0
      coscode(:,1)=cos(2*pi*rand1);sincode(:,1)=sin(2*pi*rand1)
      randw(:,(ixs-1)*n_enc_sog+1)=rand1
      do i=2,n_enc_sog
         coscode(:,i)=cos(2*pi*coscode(:,i-1))
         sincode(:,i)=sin(2*pi*sincode(:,i-1))
!               codeX((ixs-1)*n_enc_sht+i)=atan(-(sum(coscode*sincode/(sincode**2+tiny(1.)))))
         do iz=1,sou%w%n
            randw(iz,(ixs-1)*n_enc_sog+i)=atan(-sum(coscode(iz,1:i))/(sum(sincode(iz,1:i))+tiny(1.)))
         end do
      end do
   end do

do i=1,sog%xm%n
write(2,1) randw(1:100,i),'***'
end do
1 format(100f10.5,a3)
else if (gold.and.encz.and.not(perfect)) then
   allocate(randz(sog%z%n,sou%w%n,sog%xm%n))
   allocate(hold(50,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,50
      rand1=0
      call owmseq(2,9,rand1,2*sog%z%n,i-int((i-1)/50)*50)
      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))

!--------------------------------------------------
! Setting FFTW plans
!
   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)   
else if (gold.and.encz.and.perfect) then
   allocate(randz(sog%z%n,sou%w%n,sog%xm%n))
   allocate(hold(50,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)

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

   allocate(prefer(50,nq,2*sog%z%n-subf))

!   ms4(1)=ms2(1);
!   a=ms2;
!   for i=2:511; 
!       a=circshift(a,51);
!       ms4(i)=a(1);
!   end;

   do iq=1,nq
      do i=1,50
         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
      
if (exist_file("code1").and.gold.and.perfect) then
 call to_history("n1",50,"code1"); call to_history("d1",1.,"code1"); call to_history("o1",1.,"code1")  
 call to_history("n2",nq,"code1"); call to_history("d2",1.,"code1"); call to_history("o2",1.,"code1")
 call to_history("n3",2*sog%z%n-subf,"code1") ; call to_history("d3",sog%z%d,"code1") ; call to_history("o3",sog%z%o,"code1")
 call to_history("esize",4,"code1")
 call srite("code1",float(prefer),4*(2*sog%z%n-subf)*50*nq)
end if

   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
!         ic=nq-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))

!--------------------------------------------------
! Setting FFTW plans
!
   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)   

if (exist_file("code").and.gold.and.perfect) then
 call to_history("n1",sog%z%n,"code"); call to_history("d1",sog%z%d,"code"); call to_history("o1",sog%z%o,"code")  
 call to_history("n2",sou%w%n,"code"); call to_history("d2",sou%w%d/(2*pi),"code"); call to_history("o2",sou%w%o/(2*pi),"code")
 call to_history("n3",sog%xm%n,"code") ; call to_history("d3",sog%xm%d,"code") ; call to_history("o3",sog%xm%o,"code")
 call to_history("esize",4,"code")
 call srite("code",randz,4*sog%z%n*sou%w%n*sog%xm%n)
end if



else if(gold.and.encw) then
   allocate(randz(sog%z%n,sou%w%n,sog%xm%n))
   allocate(hold(50,2*sou%w%n-subf))
   allocate(gol(sog%z%n,2*sou%w%n,sog%xm%n))
   allocate(randw(sog%xm%n*sog%z%n,2))
   allocate(rand1(2*sou%w%n-subf))
   do i=1,50
      rand1=0
      call owmseq(2,9,rand1,2*sou%w%n,i-int((i-1)/50)*50)
      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,sog%z%n
         b=(randw((ixs-1)*sog%xm%n+iw,1)+1)/2;ib=int(b*50)
         c=(randw((ixs-1)*sog%xm%n+iw,2)+1)/2;ic=int(c*50)
         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,sog%z%n
         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))

!--------------------------------------------------
! Setting FFTW plans
!
   call owfftw(tfld)

   do ixs=1,sog%xm%n
      do iw=1,sog%z%n
         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 if
call owexec_time(finish)
write(0,*) "Finished generating encode functions in ",finish-start,"(s)"
end if

   if (allocated(wave_S)) deallocate(wave_S);allocate(wave_S(sou%w%n,nx,ngap));wave_S=0
   if (allocated(wave_R)) deallocate(wave_R);allocate(wave_R(sou%w%n,nx,ngap));wave_R=0
   if (allocated(GGS)) deallocate(GGS);allocate(GGS(sou%w%n,kxm%n));GGS=0
   if (allocated(GGR)) deallocate(GGR);allocate(GGR(sou%w%n,kxm%n));GGR=0

   do ixs=1,sog%xm%n

      allocate(soucig_rot(sog%z%n,sog%sh%n));soucig_rot=0.
      allocate(rectmp(sog%z%n,sog%sh%n));rectmp=0.
      allocate(soutmp(sog%z%n,sog%sh%n));soutmp=0.
      allocate(tempr(sog%z%n,sog%sh%n));tempr=0.
      allocate(temps(sog%z%n,sog%sh%n));temps=0.
   
      reccig=0.;soucig=0.
      start=0.;finish=0. 
      GGS=0.;GGR=0.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Selection scheme: adjacent cigs, in total ngap groups
! The old one is better
      if (encode.and.xiao) then
         igap=int((ixs-1)/n_enc_sog)+1
         jgap=ixs-int(ixs/n_enc_sog)*n_enc_sog
      else
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Selection scheme: cigs separated by ngap
         igap=ixs-int((ixs-1)/ngap)*ngap
         jgap=floor((ixs-1)/ngap*1.)+1
      end if

      jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
      call sreed("rec_cig",rectmp,4*sog%z%n*sog%sh%n)
      call sreed("sou_cig",soutmp,4*sog%z%n*sog%sh%n)

      do ih=1,sog%sh%n
         soucig_rot(:,ih)=soutmp(:,sog%sh%n-ih+1)
      end do

      deallocate(soutmp)

      !--------------------------------------------------
      !
      ! Upward continuing CIGs
      !
      !--------------------------------------------------     

      if (nhorz .ne. 0) then
!         allocate(window(sog%xm%n,2*nhorz))
         do ihor=1,nhorz
!            window(ixs,2*ihor-1)=minval(picks(ihor,ixs,:))
!            window(ixs,2*ihor)=maxval(picks(ihor,ixs,:))
!            np=minval(pick(ihor,:), mask=pick(ihor,:).ne.0)
!            if (np
            do iz=1,sog%z%n
               do ih=1,sog%sh%n
!               if (iz .ge. picks(ihor,ixs,ih)-wavz .and. iz .le. picks(ihor,ixs,ih)+wavz) then
               if (iz .ge. picks(ihor,ixs)-wavz .and. iz .le. picks(ihor,ixs)+wavz) then
                   tempr(iz,ih)=rectmp(iz,ih)
                   temps(iz,ih)=soucig_rot(iz,ih)
               end if
               end do
            end do
         end do
         reccig(:,jxs:jxs+sog%sh%n-1)=tempr
         soucig(:,jxs:jxs+sog%sh%n-1)=temps
      else
         reccig(:,jxs:jxs+sog%sh%n-1)=rectmp
         soucig(:,jxs:jxs+sog%sh%n-1)=soucig_rot
      end if


      deallocate(rectmp,soucig_rot,tempr,temps)

    if(exist_file("soucig")) then
      call srite("soucig",soucig,4*sog%z%n*nx)
    end if

      call owexec_time(start)
      call owctrlwei_init(ixs)

      if (exist_file("vref")) then
         if ((rand.and.encz).or.(gold.and.encz)) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,rw,randz(:,:,ixs),vref)
         else if (rand.and.encw) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,randw(:,ixs),rz,vref)
         else if (gold.and.encw) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,rw,randz(:,:,ixs),vref)
         else
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,rw,rz,vref)
         end if
      else
         if ((rand.and.encz).or.(gold.and.encz)) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,rw,randz(:,:,ixs))
         else if (rand.and.encw) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,randw(:,ixs))
         else if (gold.and.encw) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,rw,randz(:,:,ixs))
         else if (xiao) then
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS,randw(:,ixs))
         else
            call owctrlwei(taper,vsal,vel,reccig,soucig,GGR,GGS)
         end if
      end if

      call owexec_time(finish)

      write(0,3) "CIGs UPWARD PROPAGATED AT ",(ixs-1)*sog%xm%d+sog%xm%o, &
                 " # ",ixs," OUT OF ",sog%xm%n,' FROM IGAP ',igap,' in ',finish-start,' s'
3     format(A26,f10.2,A3,i3,A8,i5,a11,i4,a4,f10.2,a2)

      wave_R(:,:,igap)=wave_R(:,:,igap)+GGR
      wave_S(:,:,igap)=wave_S(:,:,igap)+GGS
   end do

   call srite("rec_out",wave_R,8*sou%w%n*nx*ngap)
   call srite("sou_out",wave_S,8*sou%w%n*nx*ngap)

close(2)
close(3)
close(4)
close(5)

deallocate(GGR,GGS,wave_R,wave_S)
deallocate(source1,vel)

call sep_end_prog()
end program OW_arealsht
