module OW_parms_areal

use sep
use OW_type

implicit none

character*3,public, save                                              :: gat
logical,public, save					:: encode,souwfld,every,encz,encw,rand,encx,sgnpi,ench,encrsxz
logical,public, save					:: gold,enclin,hadamard,perfect
type (axis) ,public, save					:: kxm, kym, w, x, y, xm, ym, z, shx, shy
type (dat_grid)	 ,public, save				:: rec, sou
type (mod_grid)	 ,public, save				:: sog, img, recev
type (velo) ,public, save					:: v
integer	 ,public, save					:: nref,iwb,ikxa,x_tap,y_tap,x_pad,y_pad,dref,isp,rndcte,subf,nreal
integer, public, save					:: imgxn, nphi, nx, ny, nx_gap, ny_gap, nwbott, shift, node
integer, public, save					:: nz,nzmax,nzmin,nhorz
real ,public, save						:: xmin,xmax,ymin,ymax,zmin,zmax,pi,dwbott
real, public, save						:: vsal,dphi,ophi,hmax,zmx,t0,eps,fac,maxang

! using lloyd
logical,public, save					:: lloyd,perc_start
real,public,save						:: min_r_pct,min_s_pct,slow_own
integer,public,save						:: niter_over

! using blocking strategy
integer,public,save						:: wbn,zbn
contains

subroutine owparam()

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

 call from_param("rand_phase",rand_phase,.false.)
if (rand_phase) then
   call from_param("nreal",nreal,1)
   call from_history("n1",sog%xm%n);  call from_history("d1",sog%xm%d);  call from_history("o1",sog%xm%o)
   call from_history("n2",sog%ym%n);  call from_history("d2",sog%ym%d);  call from_history("o2",sog%ym%o)
   call from_history("n6",nhorz,1)
   call from_aux("sou","n3",sou%w%n);call from_aux("sou","d3",sou%w%d);call from_aux("sou","o3",sou%w%o)
else
 call from_param("gat", gat)
 call from_history("n1",sog%xm%n);  call from_history("d1",sog%xm%d);  call from_history("o1",sog%xm%o)
 call from_history("n2",sog%ym%n);  call from_history("d2",sog%ym%d);  call from_history("o2",sog%ym%o)
 call from_history("n3",sog%shx%n); call from_history("d3",sog%shx%d); call from_history("o3",sog%shx%o)
 call from_history("n4",sog%shy%n); call from_history("d4",sog%shy%d); call from_history("o4",sog%shy%o)
 call from_history("n5",sog%z%n);   call from_history("d5",sog%z%d);   call from_history("o5",sog%z%o)  
 call from_history("n6",nhorz1,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","n1",v%z%n);	call from_aux("vel","d1",v%z%d);	call from_aux("vel","o1",v%z%o)

call from_param("lloyd",lloyd,.true.)
if (lloyd) then
 call from_param("nrefvel",nref) 
 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("nrefvel",nref) 
  call from_param("vsal",vsal)
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("#x source cigs .ne. #x receiver cigs")
 end if
 if (sog%ym%n .ne. recev%ym%n) then
    call seperr("#y source cigs .ne. #y receiver cigs")
 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. recev%shx%d) then
    call seperr("source x-offset interval .ne. receiver x-offset interval")
 end if
 if (sog%shy%d .ne. recev%shy%d) then
    call seperr("source y-offset interval .ne. receiver y-offset 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%d .ne. v%ym%d) then
    call seperr("source/receiv offset y-interval .ne. velocity trace y-interval")
 end if
 if (sog%xm%d .ne. recev%xm%d) then
    call seperr("source cig y-interval .ne. receiver cig y-interval")
 end if
 if (sog%ym%d .ne. recev%ym%d) then
    call seperr("source cig x-interval .ne. receiver cig x-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. image x-spacing")
 end if
 if (sog%ym%d .ne. sog%shy%d) then
    call seperr("ODCIG trace y-spacing .ne. image y-spacing")
 end if

pi=4*atan(1.)

 !!  get parameters

! memory blocking 
  call from_param("wbn",wbn,sou%w%n)
  call from_param("zbn",zbn,recev%z%n)

  call from_param("x_pad",x_pad,0)
  call from_param("y_pad",y_pad,0)
  call from_param("x_tap",x_tap,x_pad/4)
  call from_param("y_tap",y_tap,y_pad/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

  call from_param("encode",encode,.false.)
  call from_param("random",rand,.false.)
  call from_param("gold",gold,.false.)
  if (encode.and.rand.and.gold) 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("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
    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("nx_gap",nx_gap,1)
  call from_param("ny_gap",ny_gap,1)
  call from_param("node",node,9999)
  call from_param("maxang",maxang,80.)

 xmin=sog%xm%o-(x_pad+x_tap)*sog%xm%d+sog%shx%o
 xmax=((sog%xm%n-1)+(x_pad+x_tap))*sog%xm%d+sog%xm%o+(sog%shx%n-1)*sog%shx%d+sog%shx%o
 ymin=sog%ym%o-(y_pad+y_tap)*sog%ym%d+sog%shy%o
 ymax=((sog%ym%n-1)+(y_pad+y_tap))*sog%ym%d+sog%ym%o+(sog%shy%n-1)*sog%shy%d+sog%shy%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")
 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+0.5
 nzmin=(zmin-v%z%o)/sog%z%d+0.5

 write(0,*) "*************************"
 write(0,*) "* ",nz,nzmax,nzmin," *"
 write(0,*) "*************************"

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

!----------------------------------------------
!
! Wavefields
!
 call to_history("n1",(nx-2*(x_pad+x_tap)),"rec_out");   call to_history("d1",sog%shx%d,"rec_out");   call to_history("o1",xmin+(x_pad+x_tap)*sog%shx%d,"rec_out")  
 call to_history("n2",(ny-2*(y_pad+y_tap)),"rec_out");   call to_history("d2",sog%shy%d,"rec_out");   call to_history("o2",ymin+(y_pad+y_tap)*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",nreal*ngap,"rec_out");             call to_history("d4",1.,"rec_out");          call to_history("o4",1.,"rec_out")
 call to_history("n5",1,"rec_out");                      call to_history("d5",1.,"rec_out");          call to_history("o5",1.,"rec_out")  
 call to_history("esize",8,"rec_out")

 call to_history("n1",(nx-2*(x_pad+x_tap)),"sou_out");   call to_history("d1",sog%shx%d,"sou_out");   call to_history("o1",xmin+(x_pad+x_tap)*sog%shx%d,"sou_out")  
 call to_history("n2",(ny-2*(y_pad+y_tap)),"sou_out");   call to_history("d2",sog%shy%d,"sou_out");   call to_history("o2",ymin+(y_pad+y_tap)*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",nreal*ngap,"sou_out");             call to_history("d4",1.,"sou_out");          call to_history("o4",1.,"sou_out")
 call to_history("n5",1,"sou_out");	call to_history("d5",1.,"sou_out");	call to_history("o5",1.,"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",v%z%n,"velout");                   call to_history("d3",v%z%d,"velout");        call to_history("o3",v%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
! 
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 vectors
!
!kxm
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

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

write(0,*) "   PARM   "
write(0,*) "=========="
write(0,*) "PARAMETERS"
write(0,*) "=========="
write(0,*) //
write(0,*) "no.ref.velocities   = ",nref
write(0,*) "traces in taper(x,y)= ",x_tap,y_tap
write(0,*) "padded traces(x,y)  = ",x_pad,y_pad
if (encode) then
   if (rand) write(0,*) "phase encoding = ",encode," rand = ",rand
   if (gold) write(0,*) "phase encoding = ",encode," gold = ",gold
end if
write(0,*) " traces in gap(x,y) = ",nx_gap,ny_gap
write(0,*) " no.realizations    = ",nreal
write(0,*) " DATUMING:  wbott = ",wbott,"dwbott = ",dwbott,"nwbott = ",nwbott
write(0,*) " Lloyd's method vref = ",lloyd
write(0,*) "************************"
write(0,*) "*   GRID parameters    *"
write(0,*) "************************"
write(0,*) " nx = ",nx," ny = ",ny," nz = ",nz
write(0,*) " xmin = ",xmin," xmax = ",xmax
write(0,*) " ymin = ",ymin," ymax = ",ymax
write(0,*) " zmin = ",zmin," zmax = ",zmax
write(0,*) //
write(0,*) "============"
write(0,*) "HEADER WORDS"
write(0,*) "============"
write(0,*) //
write(0,*) "Rotated CIG input"
write(0,*) "================="
write(0,*) "no.offsets(x,y)  = ", sog%shx%n,sog%shy%n
write(0,*) "off.spacing(x,y) = ", sog%shx%d,sog%shy%d
write(0,*) "min.offset(x,y)  = ", sog%shx%o,sog%shy%o
write(0,*) "no.rotated CIG(x,y)= ", sog%xm%n,sog%ym%n
write(0,*) "CIG spacing(x,y) = ", sog%xm%d,sog%ym%d
write(0,*) "CIG origin (x,y) = ", sog%xm%o,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.of traces (x,y)  = ", v%xm%n, v%ym%n
write(0,*) "trace spacing(x,y)  = ", v%xm%d, v%ym%d
write(0,*) "trace origin(x,y)   = ", v%xm%o, 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,*) "FREQUENCIES"
write(0,*) "==========="
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,*) "no.wavenumbers(x,y) = ", kxm%n,kym%n
end if
end subroutine

end module
