module OW_parms_areal

use sep
use OW_type
use base_blkcal_mod

implicit none

character*6                                                           :: extrap
logical,public, save					:: encode,souwfld,every,encz,encw
logical,public, save					:: rand,encx,sgnpi,ench,encrsxz
logical,public, save					:: gold,enclin,hadamard,perfect
type (axis) ,public, save					:: kxm, kym, h, w, x, xm, z, shx, shy, ym
type (dat_grid)	 ,public, save				:: rec, sou
type (mod_grid)	 ,public, save				:: sog, recev
type (velo) ,public, save					:: v
type (block) ,public, save					:: wblock,zblock
integer	 ,public, save					:: nref,iwb,ikxa,tapertrx,tapertry
integer	 ,public, save					:: padtracesx,padtracesy,dref,isp,rndcte,subf,nreal,dipop
integer, public, save					:: imgxn, nphi, nx, ny, ngapx, ngapy, nwbott, shift, node
integer, public, save					:: nz,nzmax,nzmin,nhorz,nnx,nny
real ,public, save						:: xmin,xmax,ymin,ymax,zmin,zmax,pi,dwbott,trick
real, public, save						:: vsal,dphi,ophi,hmax,zmx,t0,eps,fac,maxang
!complex,allocatable,dimension(:),public,save			:: wavelet
! using lloyd
logical,public, save				:: lloyd,perc_start
real,public,save					:: min_r_pct,min_s_pct,slow_own
integer,public,save					:: niter_over



contains

subroutine owparam()

real							:: t1,wbott
integer							:: i,m,n,stat
integer							:: nzref,nhorz1
real							:: dzref,ozref

 call from_aux("sou_cig","n1",sog%xm%n);  call from_aux("sou_cig","d1",sog%xm%d);    call from_aux("sou_cig","o1",sog%xm%o)
 call from_aux("sou_cig","n2",sog%shx%n); call from_aux("sou_cig","d2",sog%shx%d);   call from_aux("sou_cig","o2",sog%shx%o)
 call from_aux("sou_cig","n3",sog%ym%n);  call from_aux("sou_cig","d3",sog%ym%d);    call from_aux("sou_cig","o3",sog%ym%o)
 call from_aux("sou_cig","n4",sog%shy%n); call from_aux("sou_cig","d4",sog%shy%d);   call from_aux("sou_cig","o4",sog%shy%o)
 call from_aux("sou_cig","n5",sog%z%n);   call from_aux("sou_cig","d5",sog%z%d);     call from_aux("sou_cig","o5",sog%z%o)  
 call from_aux("sou_cig","n6",nhorz1,1)

 call from_aux("rec_cig","n1",recev%xm%n);  call from_aux("rec_cig","d1",recev%xm%d);  call from_aux("rec_cig","o1",recev%xm%o)
 call from_aux("rec_cig","n2",recev%shx%n); call from_aux("rec_cig","d2",recev%shx%d); call from_aux("rec_cig","o2",recev%shx%o)
 call from_aux("rec_cig","n3",recev%ym%n);  call from_aux("rec_cig","d3",recev%ym%d);  call from_aux("rec_cig","o3",recev%ym%o)
 call from_aux("rec_cig","n4",recev%shy%n); call from_aux("rec_cig","d4",recev%shy%d); call from_aux("rec_cig","o4",recev%shy%o)
 call from_aux("rec_cig","n5",recev%z%n);   call from_aux("rec_cig","d5",recev%z%d);   call from_aux("rec_cig","o5",recev%z%o)  
 call from_aux("rec_cig","n6",nhorz,1)

 if (nhorz1.ne.nhorz) call seperr('#horz is different sou and rec')

 call from_aux("sou","n3",sou%w%n);	call from_aux("sou","d3",sou%w%d);	call from_aux("sou","o3",sou%w%o)

 call from_aux("vel","n1",v%xm%n);	call from_aux("vel","d1",v%xm%d);	call from_aux("vel","o1",v%xm%o)
 call from_aux("vel","n2",v%ym%n);	call from_aux("vel","d2",v%ym%d);	call from_aux("vel","o2",v%ym%o)
 call from_aux("vel","n3",v%z%n);	call from_aux("vel","d3",v%z%d);	call from_aux("vel","o3",v%z%o)
 call from_param("extrapolator",extrap,"ffd")
if (extrap(1:4).eq."pspi") then
   call from_param("lloyd",lloyd,.true.)
   call from_param("nrefvel",nref,1) 
   if (lloyd) then
      call from_param("min_region_pct",min_r_pct,1.)
      min_r_pct=min_r_pct/100.
     !MINIMUM ALLOWED DEVIATION BETWEEN CENTROIDS
      call from_param("min_slow_dev",min_s_pct,9.)
      min_s_pct=min_s_pct/100.
      call from_param("perc_start",perc_start,.false.)
      call from_param("vel_own",slow_own,-1.)
      slow_own=1./slow_own
     !NUMBER OF ITERATIONS
      call from_param("niter_lloyd",niter_over,-1)
   else
      call from_param("vsal",vsal)
   end if
else
   call from_param("trick",trick,0.14)
   call from_param("dipop",dipop,45)
end if
 if (sog%z%n .ne. recev%z%n) then
    call seperr("# source depth. .ne. # receiver depth.")
 end if
 if (sog%xm%n .ne. recev%xm%n) then
    call seperr("# source xlines .ne. # receiver xlines")
 end if
 if (sog%ym%n .ne. recev%ym%n) then
    call seperr("# source ilines .ne. # receiver ilines")
 end if
 if (sog%shx%n .ne. recev%shx%n) then
    call seperr("# x-source offsets. .ne. # x-receiver offsets")
 end if
 if (sog%shy%n .ne. recev%shy%n) then
    call seperr("# y-source offsets. .ne. # y-receiver offsets")
 end if
 if (sog%z%d .ne. recev%z%d) then
    call seperr("source depth.interval .ne. receiver depth.interval")
 end if
 if (sog%shx%d .ne. v%xm%d) then
    call seperr("source/receiv x-offset interval .ne. velocity trace x-interval")
 end if
 if (sog%shy%n.ne.1) then
    if (sog%shy%d .ne. v%ym%d) then
       call seperr("source/receiv y-offset interval .ne. velocity trace y-interval")
    end if
 end if
 if (sog%xm%d .ne. recev%xm%d) then
    call seperr("source cig interval .ne. receiver cig interval")
 end if
 if (sog%ym%d .ne. recev%ym%d) then
    call seperr("source cig y-interval .ne. receiver cig y-interval")
 end if
 if (sog%z%o .ne. recev%z%o) then
    call seperr("source depth origin .ne. receiver depth origin")
 end if
 if (sog%xm%o .ne. recev%xm%o) then
    call seperr("source cig x-origin .ne. receiver cig x-origin")
 end if
 if (sog%ym%o .ne. recev%ym%o) then
    call seperr("source cig y-origin .ne. receiver cig y-origin")
 end if
 if (sog%shx%o .ne. recev%shx%o) then
    call seperr("source offset x-origin .ne. receiver offset x-origin")
 end if
 if (sog%shy%o .ne. recev%shy%o) then
    call seperr("source offset y-origin .ne. receiver offset y-origin")
 end if

 if (sog%z%d .ne. v%z%d) then
    call seperr("source/receiv depth interval .ne. velocity depth interval")
 end if
 if (sog%xm%d .ne. sog%shx%d) then
    call seperr("ODCIG trace x-spacing .ne. offset x-spacing")
 end if
 if (sog%shy%n.ne.1) then
    if (sog%ym%d .ne. sog%shy%d) then
       call seperr("ODCIG trace y-spacing .ne. offset y-spacing")
    end if
 end if

pi=4*atan(1.)

 !!  get parameters

  call from_param("padtraces_x",padtracesx,0)
  call from_param("traces_in_taper_x",tapertrx,padtracesx/4)
  call from_param("padtraces_y",padtracesy,0)
  call from_param("traces_in_taper_y",tapertry,padtracesy/4)
  call from_param("wbott",wbott,0.)
  wbott=int(wbott/sog%z%d)*sog%z%d
  nwbott=int(wbott/sog%z%d)+1
  dwbott=wbott
 write(0,*) "wbott ",wbott,"dwbott ",dwbott,"nwbott ",nwbott
  call from_param("encode",encode,.false.)
  call from_param("random",rand,.false.)
  call from_param("gold",gold,.false.)
  call from_param("sgnpi",sgnpi,.false.)
  if (rand.and.gold.or.rand.and.sgnpi.or.gold.and.sgnpi) call seperr("Choose only one method: random or gold or sgnpi!")  

 subf=0
 call from_param("eps",eps,1.)
 call from_param("nreal",nreal,1)

if (encode) then
   call from_param("enclin",enclin,.false.)
   if (enclin) then
      t1=(2*pi/sou%w%d+2/(sou%w%n))/2
      call from_param("t0",t0,t1)
      write(0,*) "encode = ",encode," enclin = ",enclin," t1 = ",t1
   end if
   call from_param("encz",encz,.false.)
   if (encz) then
      if (rand.and.encz) then
         call from_param("rndcte",rndcte,1371)
         write(0,*) "encode = ",encode," encz = ",encz," eps = ",eps
      else if (gold) then
         call from_param("shift",shift,5)
         call from_param("rndcte",rndcte,137*shift)
         call from_param("fac",fac,1.)
         subf=1
         write(0,*) "encode = ",encode," gold = ",gold," eps = ",eps,"  shift = ",shift,"  fac = ",fac
      end if
    end if
    call from_param("encw",encw,.false.)
    if (encw) then
       if (rand) then
          call from_param("rndcte",rndcte,1371)
          write(0,*) "encode = ",encode," encw = ",encw," eps = ",eps
          subf=0
       else if (gold) then
          call from_param("shift",shift,5)
          call from_param("rndcte",rndcte,1373)
          call from_param("fac",fac,1.)
          subf=1
          write(0,*) "encode = ",encode," gold = ",gold," eps = ",eps,"  shift = ",shift
       end if
    end if
    call from_param("encx",encx,.false.)
    if (encx) then
       if (rand) then
          call from_param("rndcte",rndcte,1371)
          write(0,*) "encode = ",encode," encx = ",encx," eps = ",eps
          subf=0
       end if
    end if
    call from_param("ench",ench,.false.)
    if (ench) then
        call from_param("rndcte",rndcte,1371)
        write(0,*) "encode = ",encode," ench = ",ench," eps = ",eps
        subf=0
    end if
    call from_param("encrsxz",encrsxz,.false.)
    if (ench) then
        call from_param("rndcte",rndcte,1371)
        write(0,*) "encode = ",encode," encrsxz = ",encrsxz," eps = ",eps
        subf=0
    end if
    if (not(ench.or.encx.or.gold.or.encw.or.encz.or.enclin.or.encrsxz.or.exist_file("InCode"))) then
        write(0,*) 'How to encode the areal shots ? Modeling with no encode'
    end if
end if
  call from_param("ngapx",ngapx,1)
  call from_param("ngapy",ngapy,1)
  call from_param("node",node,9999)
  call from_param("maxang",maxang,80.)

 xmin=sog%xm%o-(padtracesx+tapertrx)*sog%xm%d+sog%shx%o
 xmax=((sog%xm%n-1)+(padtracesx+tapertrx))*sog%xm%d+sog%xm%o+(sog%shx%n-1)*sog%shx%d+sog%shx%o
 if (xmin < v%xm%o ) call seperr("Min.x-image < Min.x-velocity")
 if (xmax > v%xm%o+(v%xm%n-1)*v%xm%d ) call seperr("Max.x-image > Max.x-velocity")

 ymin=sog%ym%o-(padtracesy+tapertry)*sog%ym%d+sog%shy%o
 ymax=((sog%ym%n-1)+(padtracesy+tapertry))*sog%ym%d+sog%ym%o+(sog%shy%n-1)*sog%shy%d+sog%shy%o
 if (ymin < v%ym%o ) call seperr("Min.y-image < Min.y-velocity")
 if (ymax > v%ym%o+(v%ym%n-1)*v%ym%d ) call seperr("Max.y-image > Max.y-velocity")

 call from_param("zmin",zmin,sog%z%o)
 call from_param("zmax",zmax,(sog%z%n-1)*sog%z%d+sog%z%o)

 if (zmin < v%z%o ) call seperr("Min.z-image < Min.z-velocity")
 if (zmax > v%z%o+(v%z%n-1)*v%z%d ) call seperr("Max.z-image > Max.z-velocity")

 nx=(xmax-xmin)/sog%xm%d+1.5
 ny=(ymax-ymin)/sog%ym%d+1.5
 nz=(zmax-zmin)/sog%z%d+1.5
 nzmax=(zmax-v%z%o)/sog%z%d+1.5 
 nzmin=int((zmin-v%z%o)/v%z%d)+1.5
 write(0,*) "*************************"
 write(0,*) "* ",nz,nzmax,nzmin," *"
 write(0,*) "*************************"
 write(0,*) "* ",nx,ny,nz," *"
 nnx=int(sog%xm%n/ngapx)
 nny=int(sog%ym%n/ngapy)

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

!----------------------------------------------
!
! Wavefields
!
 call to_history("n1",(nx-2*(padtracesx+tapertrx)),"rec_out");call to_history("d1",sog%shx%d,"rec_out");call to_history("o1",xmin+(padtracesx+tapertrx)*sog%shx%d,"rec_out")  
 call to_history("n2",(ny-2*(padtracesy+tapertry)),"rec_out");call to_history("d2",sog%shy%d,"rec_out");call to_history("o2",ymin+(padtracesy+tapertry)*sog%shy%d,"rec_out")  
 call to_history("n3",sou%w%n,"rec_out");call to_history("d3",sou%w%d,"rec_out");call to_history("o3",sou%w%o,"rec_out")
 call to_history("n4",ngapx*ngapy,"rec_out");	call to_history("d4",1.,"rec_out");	call to_history("o4",1.,"rec_out")
 call to_history("n5",nreal,"rec_out");	call to_history("d5",1.,"rec_out");	call to_history("o5",1.,"rec_out") 
 call to_history("n6",1,"rec_out");	call to_history("d6",1.,"rec_out");	call to_history("o6",zmin,"rec_out") 
 call to_history("n7",1,"rec_out");	call to_history("d7",1.,"rec_out");	call to_history("o7",0.,"rec_out") 
 call to_history("esize",8,"rec_out")

 call to_history("n1",(nx-2*(padtracesx+tapertrx)),"sou_out");call to_history("d1",sog%shx%d,"sou_out");call to_history("o1",xmin+(padtracesx+tapertrx)*sog%shx%d,"sou_out")  
 call to_history("n2",(ny-2*(padtracesy+tapertry)),"sou_out");call to_history("d2",sog%shy%d,"sou_out");call to_history("o2",ymin+(padtracesy+tapertry)*sog%shy%d,"sou_out")  
 call to_history("n3",sou%w%n,"sou_out");	call to_history("d3",sou%w%d,"sou_out");	call to_history("o3",sou%w%o,"sou_out")
 call to_history("n4",ngapx*ngapy,"sou_out");	call to_history("d4",1.,"sou_out");	call to_history("o4",1.,"sou_out")
 call to_history("n5",nreal,"sou_out");	call to_history("d5",1.,"sou_out");	call to_history("o5",1.,"sou_out")  
 call to_history("n6",1,"sou_out");	call to_history("d6",1.,"sou_out");	call to_history("o6",zmin,"sou_out") 
 call to_history("n7",1,"sou_out");	call to_history("d7",1.,"sou_out");	call to_history("o7",0.,"sou_out") 
 call to_history("esize",8,"sou_out")

if(exist_file("velout")) then
 call to_history("n1",nx,"velout"); call to_history("d1",sog%shx%d,"velout"); call to_history("o1",xmin,"velout")
 call to_history("n2",ny,"velout"); call to_history("d2",sog%shy%d,"velout"); call to_history("o2",ymin,"velout")
 call to_history("n3",sog%z%n,"velout") ; call to_history("d3",sog%z%d,"velout") ; call to_history("o3",sog%z%o,"velout")
 call to_history("n4",1,"velout"); call to_history("d4",1.,"velout"); call to_history("o4",0.,"velout")
 call to_history("n5",1,"velout"); call to_history("d5",1.,"velout"); call to_history("o5",0.,"velout")  
 call to_history("esize",4,"velout")
end if
!--------------------------------------------------
! Frequency vector
! no zero frequency component
!
sou%w%o   = 2*pi * sou%w%o
sou%w%d   = 2*pi * sou%w%d
allocate(sou%w%all(sou%w%n))
do iwb = 1,sou%w%n
   sou%w%all(iwb) = sou%w%o + (iwb - 1) * sou%w%d
end do

!--------------------------------------------------
! Wavenumber vector - Kx
!
kxm%n = nx
allocate(kxm%all(kxm%n))
kxm%all=0.
if(kxm%n/=1) then
  m = kxm%n/2
  do i=1,m+1
     kxm%all(i)=i-1
  end do
  do i=m+2,kxm%n
     kxm%all(i)=i-kxm%n-1
  end do
else
  kxm%all=0. !! make sure k=0 if n=1
end if

kxm%all = kxm%all / (kxm%n*sog%xm%d)      !! 1/length
kxm%all = 2 * PI * kxm%all     !! radians

!--------------------------------------------------
! Wavenumber vector - Ky
!
kym%n = ny
allocate(kym%all(kym%n))
kym%all=0.
if(kym%n/=1) then
  m = kym%n/2
  do i=1,m+1
     kym%all(i)=i-1
  end do
  do i=m+2,kym%n
     kym%all(i)=i-kym%n-1
  end do
else
  kym%all=0. !! make sure k=0 if n=1
end if

kym%all = kym%all / (kym%n*sog%ym%d)      !! 1/length
kym%all = 2 * PI * kym%all     !! radians

!=====================
! Defining blocks
! w and z
call from_param("nwblk",wblock%ns,20)
call from_param("nzblk",zblock%ns,20)
call base_blkcal(wblock, sou%w)
call base_blkcal(zblock, sog%z)

write(0,*) "   PARM   "
write(0,*) "=========="
write(0,*) "PARAMETERS"
write(0,*) "=========="
write(0,*) //
if (lloyd) then
   write(0,*) "no.ref.velocities   = ",nref
   write(0,*) "max.limit velocity  = ",vsal
end if
write(0,*) "traces in taper x      = ",tapertrx
write(0,*) "padded traces   x      = ",padtracesx
write(0,*) "traces in taper y      = ",tapertry
write(0,*) "padded traces   y      = ",padtracesy
write(0,*) "phase encoding         = ",encode
write(0,*) "no.traces in gap x     = ",ngapx
write(0,*) "no.traces in gap x     = ",ngapy
write(0,*) //
write(0,*) "============"
write(0,*) "HEADER WORDS"
write(0,*) "============"
write(0,*) //
write(0,*) "Rotated CIG input"
write(0,*) "================="
write(0,*) "no.x-offsets      = ", sog%shx%n
write(0,*) "x-offset spacing  = ", sog%shx%d
write(0,*) "min.-xoffset      = ", sog%shx%o
write(0,*) "no.xlines       = ", sog%xm%n
write(0,*) "xline spacing   = ", sog%xm%d
write(0,*) "min.CIG x-coord = ", sog%xm%o
write(0,*) "no.y-offsets    = ", sog%shy%n
write(0,*) "y-offset spacing= ", sog%shy%d
write(0,*) "min.y-offset    = ", sog%shy%o
write(0,*) "no. ilines      = ", sog%ym%n
write(0,*) "iline spacing   = ", sog%ym%d
write(0,*) "min.CIG y-coord = ", sog%ym%o
write(0,*) "no.depth steps  = ", nz
write(0,*) "depth interval  = ", sog%z%d
write(0,*) "depth origin    = ", zmin
write(0,*) //
write(0,*) "VELOCITY MODEL"
write(0,*) "=============="
write(0,*) "no.xlines       = ", v%xm%n
write(0,*) "xline spacing   = ", v%xm%d
write(0,*) "min.trace-x     = ", v%xm%o
write(0,*) "no.ilines       = ", v%ym%n
write(0,*) "iline spacing   = ", v%ym%d
write(0,*) "min.trace-y     = ", v%ym%o
write(0,*) "no.depth samples= ", v%z%n
write(0,*) "depth interval  = ", v%z%d
write(0,*) "minimum depth   = ", v%z%o
write(0,*) //
write(0,*) "AREAL SPACE"
write(0,*) "==========="
write(0,*) "no.xlines       = ", nx
write(0,*) "xline spacing   = ", sog%shx%d
write(0,*) "min.trace-x     = ", xmin
write(0,*) "no.ilines       = ", ny
write(0,*) "iline spacing   = ", sog%shy%d
write(0,*) "min.trace-y     = ", ymin
write(0,*) "no.frequencies  = ", sou%w%n
write(0,*) "freq.interval   = ", sou%w%d/(2*pi)
write(0,*) "min.frequency   = ", sou%w%o/(2*pi)
write(0,*) //
write(0,*) "WAVENUMBER"
write(0,*) "=========="
write(0,*) "kxm%n           = ", kxm%n
write(0,*) "kym%n           = ", kym%n

end subroutine

end module
