module OW_parms_areal

use sep
use OW_type

implicit none

character*6,public, save                                    :: extrap
logical,public, save				:: adj,dot,marine
logical,public, save				:: light
logical,public, save				:: areal
type (axis) ,public, save				:: kxm, h, w, xd, xm, z, sh
type (dat_grid)	 ,public, save			:: rec, sou
type (mod_grid)	 ,public, save			:: img
type (velo) ,public, save				:: v
integer	 ,public, save				:: nref,tapertr,padtraces,nz,nzmin,nzmax,nsecdim
integer	 ,public, save				:: isp,nwbott,node
real ,public, save					:: ymin,ymax,zmin,zmax,pi,dwbott
real ,public, save					:: osecdim,dsecdim
real, public, save					:: vsal,maxang
real ,public, save					:: trick
integer, public, save				:: nah,dipop,iwb
real, public, save					:: dah,oah
! 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						:: wbott,dd1,dd2,od1,od2,dv1,dv2,ov1,ov2
real						:: deltat, d_sampt,dzref,ozref
integer						:: iw,nt, i, m,nd1,nd2,nv1,nv2,nzref,jh
logical						:: inversion,migration

call from_param("dot",dot,.false.)
call from_param("marine",marine,.false.)
call from_param("adj",adj,.true.)
call from_param("extrapolator",extrap,"ffd")
if (extrap(1:3).eq."ffd") then
   call from_param("dipop",dipop,45)
   call from_param("trick",trick,0.14)
end if
call from_param("areal",areal,.false.)
call from_param("zmin",zmin,0.)
call from_param("zmax",zmax,(v%z%n-1)*v%z%d+v%z%o)
call from_param("node",node,9999)
call from_param("traces_in_taper",tapertr,0)
call from_param("padtraces",padtraces,tapertr)
call from_param("maxang",maxang,89.)

 call from_aux("sou","n1",sou%h%n);	call from_aux("sou","d1",sou%h%d);	call from_aux("sou","o1",sou%h%o)
 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("sou","n4",sou%xd%n);	call from_aux("sou","d4",sou%xd%d);	call from_aux("sou","o4",sou%xd%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","n3",v%z%n);	call from_aux("vel","d3",v%z%d);	call from_aux("vel","o3",v%z%o)

if (adj) then

   if (not(marine).or.areal) then
      img%xm%n=rec%h%n+2*padtraces+2*tapertr
      img%xm%d=rec%h%d
      img%xm%o=rec%h%o-(padtraces+tapertr)*img%xm%d
   else
      img%xm%n=rec%h%n+2*padtraces+2*tapertr+(rec%xd%n-1)*rec%xd%d/rec%h%d
      img%xm%d=rec%h%d
      img%xm%o=rec%h%o-(padtraces+tapertr)*img%xm%d+rec%xd%o
   end if

   nz=int((zmax-zmin)/v%z%d)+1.5
   nzmax=int((zmax-v%z%o)/v%z%d)+1.5
   nzmin=int((zmin-v%z%o)/v%z%d)+1.5
   img%z%o=zmin
   img%z%n=nz
   img%z%o=zmin
   img%z%d=v%z%d

   call from_param("nh",img%sh%n,1)
   call from_param("jh",jh,1)
   call from_param("dh",img%sh%d,jh*img%xm%d)
   call from_param("oh",img%sh%o,-int(img%sh%n/2)*img%sh%d)
   nsecdim=img%sh%n
   dsecdim=img%sh%d
   osecdim=img%sh%o
   call auxclose("bkgimg")
   call auxclose("rec")
   call auxinout("bkgimg")
   call auxinout("rec")
   call to_history("n1",img%xm%n,"bkgimg"); call to_history("d1",img%xm%d,"bkgimg"); call to_history("o1",img%xm%o,"bkgimg")  
   call to_history("n2",img%sh%n,"bkgimg") ; call to_history("d2",img%sh%d,"bkgimg") ; call to_history("o2",img%sh%o,"bkgimg")
   call to_history("n3",nz,"bkgimg") ; call to_history("d3",img%z%d,"bkgimg") ; call to_history("o3",zmin,"bkgimg")
   call to_history("n4",1,"bkgimg"); call to_history("d4",1,"bkgimg"); call to_history("o4",0,"bkgimg")
   call to_history("n5",1,"bkgimg"); call to_history("d5",1,"bkgimg"); call to_history("o5",0,"bkgimg")  
   call to_history("esize",4,"bkgimg")
   call from_aux("rec","n1",rec%h%n);	call from_aux("rec","d1",rec%h%d);call from_aux("rec","o1",rec%h%o)
   call from_aux("rec","n3",rec%w%n);	call from_aux("rec","d3",rec%w%d);call from_aux("rec","o3",rec%w%o)
   call from_aux("rec","n4",rec%xd%n);call from_aux("rec","d4",rec%xd%d);call from_aux("rec","o4",rec%xd%o)  
   if (sou%w%n .ne. rec%w%n) then
      call seperr("# source freq. .ne. # receiver freq.")
   end if
   if (sou%xd%n.ne.1) then
      if (sou%xd%n .ne. rec%xd%n) then
         call seperr("# source planes. .ne. # receiver planes")
      end if
      if (sou%h%n .ne. rec%h%n) then
         call seperr("# source traces. .ne. # receiver traces")
      end if
      if (sou%h%d .ne. rec%h%d) then
         call seperr("source trace interval .ne. areal shot trace interval")
      end if
   end if
   if (sou%w%d .ne. rec%w%d) then
      call seperr("source freq.interval .ne. receiver freq.interval")
   end if
   if (sou%w%o .ne. rec%w%o) then
      call seperr("source freq.origin .ne. receiver freq.origin")
   end if
else
   call auxclose("bkgimg")
   call auxclose("rec")
   call auxinout("bkgimg")
   call auxinout("rec")
   call from_aux("bkgimg","n1",img%xm%n); call from_aux("bkgimg","d1",img%xm%d); call from_aux("bkgimg","o1",img%xm%o)  
   call from_aux("bkgimg","n2",img%sh%n) ; call from_aux("bkgimg","d2",img%sh%d) ; call from_aux("bkgimg","o2",img%sh%o)
   call from_aux("bkgimg","n3",nz) ; call from_aux("bkgimg","d3",img%z%d) ; call from_aux("bkgimg","o3",zmin)
   img%z%n=nz
   zmax=int((img%z%n-1)*img%z%d+img%z%o)
   nzmax=int((zmax-v%z%o)/v%z%d)+1.5
   nzmin=int((zmin-v%z%o)/v%z%d)+1.5
   call from_param("nrec",rec%h%n);	call from_param("drec",rec%h%d);	call from_param("orec",rec%h%o)
   call from_param("nsht",rec%xd%n);	call from_param("dsht",rec%xd%d);	call from_param("osht",rec%xd%o)  
   rec%w%n=sou%w%n; rec%w%d=sou%w%d; rec%w%o=sou%w%o

   call to_history("n1",rec%h%n,"rec");	call to_history("d1",rec%h%d,"rec");
   if (marine) then
      call to_history("o1",rec%h%o,"rec")
   else   
      call to_history("o1",rec%h%o+rec%xd%o,"rec")
   end if
   call to_history("n2",1,"rec");	call to_history("d2",1.,"rec");call to_history("o2",.0,"rec")
   call to_history("n3",rec%w%n,"rec");	call to_history("d3",rec%w%d,"rec");call to_history("o3",rec%w%o,"rec")
   call to_history("n4",rec%xd%n,"rec");call to_history("d4",rec%xd%d,"rec");call to_history("o4",rec%xd%o,"rec")  
   call to_history("n5",1,"rec");	call to_history("d5",1.,"rec");call to_history("o5",.0,"rec")
   call to_history("esize",8,"rec")
   nsecdim=img%sh%n
   dsecdim=img%sh%d
   osecdim=img%sh%o
   if (exist_file("bws")) then
     call to_history("n1",rec%h%n,"bws");	call to_history("d1",rec%h%d,"bws");
     if (marine) then
        call to_history("o1",rec%h%o,"bws")
     else   
        call to_history("o1",rec%h%o+rec%xd%o,"bws")
     end if
     call to_history("n2",img%z%n,"bws");	call to_history("d2",img%z%d,"bws");call to_history("o2",img%z%o,"bws")
     call to_history("n3",rec%w%n,"bws");	call to_history("d3",rec%w%d,"bws");call to_history("o3",rec%w%o,"bws")
     call to_history("n4",rec%xd%n,"bws");call to_history("d4",rec%xd%d,"rec");call to_history("o4",rec%xd%o,"rec")  
     call to_history("n5",1,"bws");	call to_history("d5",1.,"bws");call to_history("o5",.0,"bws")
     call to_history("esize",8,"bws")
   end if  
   if (rec%h%o+rec%xd%o-(padtraces+tapertr)*rec%h%d < v%xm%o) call seperr("minx_data < minx_velocity")
!   if (rec%h%o+rec%xd%o-(padtraces+tapertr)*rec%h%d + < v%xm%o) call seperr("minx_data < minx_velocity")
end if
if (extrap(1:4).eq."pspi") then
   call from_param("nrefvel",nref,1) 
   call from_param("lloyd",lloyd,.false.)
   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
end if

! to be changed for real 3D
  img%ym%n=1
  img%ym%d=1.
  img%ym%o=0.

  if (img%xm%o < v%xm%o ) call seperr("Min.x-image space < Min.x-velocity")
  if ((img%xm%n-1)*img%xm%d+img%xm%o > v%xm%o+(v%xm%n-1)*v%xm%d ) call seperr("Max.x-image space > Max.x-velocity")
 
  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")

  wbott=int(wbott/v%z%d)*v%z%d
  nwbott=int(wbott/v%z%d)+1
  dwbott=wbott/2.

pi=4*atan(1.)
isp=abs(floor(rec%h%o/rec%h%d)) + 1

kxm%n = rec%h%n+2*(tapertr+padtraces)
kxm%d = 1./(kxm%n*rec%h%d)
kxm%o = 1./(2*rec%h%d)-kxm%d
allocate(kxm%all(kxm%n))
do i=1,kxm%n
   kxm%all(i) = mod(kxm%o+(i-1)*kxm%d,1./rec%h%d)-kxm%o
end do
kxm%all = 2*pi * kxm%all


!--------------------------------------------------
! Frequency vector
! no zero frequency component
!
rec%w%o   = 2*pi * rec%w%o
rec%w%d   = 2*pi * rec%w%d
allocate(rec%w%all(rec%w%n))
do iwb = 1,rec%w%n
   rec%w%all(iwb) = rec%w%o + (iwb - 1) * rec%w%d
end do

!----------------------------------------------------
write(0,*) "   PARM   "
write(0,*) "=========="
write(0,*) "PARAMETERS"
write(0,*) "=========="
write(0,*) //
write(0,*) "Marine              = ",marine
write(0,*) "no.ref.velocities   = ",nref
write(0,*) "min.x-coord         = ",img%xm%o
write(0,*) "max.x-coord         = ",(img%xm%n-1)*img%xm%d+img%xm%o
write(0,*) "min.y-coord         = ",ymin
write(0,*) "max.y-coord         = ",ymax
write(0,*) "min.depth           = ",zmin
write(0,*) "max.depth           = ",zmax
write(0,*) "traces in taper     = ",tapertr
write(0,*) "padded traces       = ",padtraces
write(0,*) //
write(0,*) "============"
write(0,*) "HEADER WORDS"
write(0,*) "============"
write(0,*) //
write(0,*) "SHOTS"
write(0,*) "====="
write(0,*) "number of recv  = ", rec%h%n
write(0,*) "recv spacing    = ", rec%h%d
write(0,*) "min.recv offset = ", rec%h%o
write(0,*) "number of freq  = ", rec%w%n
write(0,*) "freq spacing    = ", rec%w%d/(2*pi)
write(0,*) "minimum freq    = ", rec%w%o/(2*pi)
write(0,*) "number of shots = ", rec%xd%n
write(0,*) "shot spacing    = ", rec%xd%d
write(0,*) "min.shot-x      = ", rec%xd%o
write(0,*) //
write(0,*) "SOURCE FUNCTION"
write(0,*) "==============="
write(0,*) "no.of traces    = ", sou%h%n
write(0,*) "trace spacing   = ", sou%h%d
write(0,*) "min.trace-x     = ", sou%h%o
write(0,*) "number of freq  = ", sou%w%n
write(0,*) "freq spacing    = ", sou%w%d
write(0,*) "minimum freq    = ", sou%w%o
write(0,*) "no.of s.fnct    = ", sou%xd%n
write(0,*) "s.fnct spacing  = ", sou%xd%d
write(0,*) "min.s.fnct-x    = ", sou%xd%o
write(0,*) //
write(0,*) "VELOCITY MODEL"
write(0,*) "=============="
write(0,*) "no.of traces    = ", v%xm%n
write(0,*) "trace spacing   = ", v%xm%d
write(0,*) "min.trace-x     = ", v%xm%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,*) "IMAGE SPACE"
write(0,*) "==========="
write(0,*) "no.of traces    = ", img%xm%n
write(0,*) "trace spacing   = ", img%xm%d
write(0,*) "min.trace-x     = ", img%xm%o
write(0,*) "no.of offsets    = ", img%sh%n
write(0,*) "offset spacing   = ", img%sh%d
write(0,*) "offset origin     = ", img%sh%o
write(0,*) "no.depth samples= ", img%z%n
write(0,*) "depth interval  = ", v%z%d
write(0,*) "minimum depth   = ", img%z%o

write(0,*) //
write(0,*) "WAVENUMBER"
write(0,*) "=========="
write(0,*) "no.wavenumbers  = ", kxm%n
write(0,*) "waveno.spacing  = ", kxm%d
write(0,*) "waveno.origin   = ", kxm%o
write(0,*) //

write(0,*) "*****************************************************"
end subroutine

end module
