module OW_parms_areal

use sep
use OW_type

implicit none

logical,public, save					:: encode,souwfld,every,encz,encw,combine,enclin
type (axis) ,public, save					:: kxm, h, w, x, xm, z, sh
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,tapertr,padtraces,dref,isp,rndcte
integer, public, save					:: imgxn,nphi,nx,ngap,seed1,seed2,funct
real ,public, save						:: xmin,xmax,ymin,ymax,zmin,zmax,pi,wbott
real, public, save						:: vsal,dphi,ophi,hmax,zmx,t0,eps

contains

subroutine owparam(migrate)

logical							:: migrate
real							:: t1
integer							:: i,m

if (not(migrate)) then

 call from_aux("sou_cig","n1",sog%z%n);		call from_aux("sou_cig","d1",sog%z%d);	call from_aux("sou_cig","o1",sog%z%o)
 call from_aux("sou_cig","n2",sog%sh%n);	call from_aux("sou_cig","d2",sog%sh%d);	call from_aux("sou_cig","o2",sog%sh%o)
 call from_aux("sou_cig","n3",sog%xm%n);	call from_aux("sou_cig","d3",sog%xm%d);	call from_aux("sou_cig","o3",sog%xm%o)  

 call from_aux("rec_cig","n1",recev%z%n);	call from_aux("rec_cig","d1",recev%z%d);	call from_aux("rec_cig","o1",recev%z%o)
 call from_aux("rec_cig","n2",recev%sh%n);	call from_aux("rec_cig","d2",recev%sh%d);call from_aux("rec_cig","o2",recev%sh%o)
 call from_aux("rec_cig","n3",recev%xm%n);	call from_aux("rec_cig","d3",recev%xm%d);call from_aux("rec_cig","o3",recev%xm%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("vel","n1",v%z%n);	call from_aux("vel","d1",v%z%d);	call from_aux("vel","o1",v%z%o)
 call from_aux("vel","n2",v%xm%n);	call from_aux("vel","d2",v%xm%d);	call from_aux("vel","o2",v%xm%o)

 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 cigs .ne. # receiver cigs")
 end if
 if (sog%sh%n .ne. recev%sh%n) then
    call seperr("# source offsets. .ne. # 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%sh%d .ne. recev%sh%d) then
    call seperr("source offset interval .ne. receiver offset interval")
 end if
 if (sog%sh%d .ne. v%xm%d) then
    call seperr("source/receiv offset interval .ne. velocity trace interval")
 end if
 if (sog%xm%d .ne. recev%xm%d) then
    call seperr("source cig interval .ne. receiver cig 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 origin .ne. receiver cig origin")
 end if
 if (sog%sh%o .ne. recev%sh%o) then
    call seperr("source offset origin .ne. receiver offset 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%sh%d) then
    call seperr("ODCIG trace spacing .ne. image spacing")
 end if

pi=4*atan(1.)

 !!  get parameters
  call from_param("nrefvel",nref) 
  call from_param("vsal",vsal)
!  call from_param("traces_in_taper",tapertr,0)
  call from_param("padtraces",padtraces,0)
!  if (padtraces.eq.0) then
!     write(0,*) "Must pad traces around the CIG to account for the propagation"
!     call seperr("It depends on the focusing, complexity of the geology ...")
!  end if
  call from_param("encode",encode)
  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)
     else
        call from_param("encz",encz,.false.)
        if (encz) then
           call from_param("rndcte",rndcte,1371)
           call from_param("funct",funct,0)
           call from_param("eps",eps,0.05)
        else
           call from_param("rndcte",rndcte,1371)
           call from_param("encw",encw,.false.)
           if (not(encw)) call seperr('How to encode the areal shots ?')
        end if
     end if
  end if
  call from_param("ngap",ngap,sog%sh%n)
  call from_param("combine",combine,.true.)
  call from_param("every",every,.false.)

 
 xmin=sog%xm%o-padtraces*sog%xm%d+sog%sh%o
 xmax=((sog%xm%n-1)+padtraces)*sog%xm%d+sog%xm%o+(sog%sh%n)*sog%sh%d+sog%sh%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")

 zmin=sog%z%o
 zmax=(sog%z%n-1)*sog%z%d+sog%z%o
 if (zmin < v%xm%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=int((xmax-xmin)/sog%xm%d)+1.5

  write(0,*) "xmin ",xmin
  write(0,*) "xmax ",xmax
  write(0,*) "nx ",nx
  write(0,*) "zmin ",zmin
  write(0,*) "zmax ",zmax
  write(0,*) "nrefvel ",nref
  write(0,*) "vsal ",vsal
  write(0,*) "traces_in_taper ",tapertr
  write(0,*) "padtraces ",padtraces
  write(0,*) "encode ",encode
  if (encode.and.encw) then; write(0,*) "encw  ",encw
  else if (encode.and.encz) then; write(0,*) "encz  ",encz
  end if
  write(0,*) "t0 ",t0
  write(0,*) "ngap ",ngap

!----------------------------------------------
!
! Wavefield at z=0
!

 call to_history("n1",sou%w%n,"rec_out");	call to_history("d1",sou%w%d,"rec_out");	call to_history("o1",sou%w%o,"rec_out")
 call to_history("n2",1,"rec_out");	call to_history("d2",1,"rec_out");	call to_history("o2",0.,"rec_out")
 call to_history("n3",nx,"rec_out");	call to_history("d3",sog%sh%d,"rec_out");	call to_history("o3",xmin,"rec_out")  
 if (combine) then
   call to_history("n4",ngap,"rec_out");	call to_history("d4",1.,"rec_out");	call to_history("o4",1.,"rec_out")
 else
   call to_history("n4",sog%xm%n,"rec_out");	call to_history("d4",1.,"rec_out");	call to_history("o4",1.,"rec_out")
 end if
 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",sou%w%n,"sou_out");	call to_history("d1",sou%w%d,"sou_out");	call to_history("o1",sou%w%o,"sou_out")
 call to_history("n2",1,"sou_out");	call to_history("d2",1,"sou_out");	call to_history("o2",0.,"sou_out")
 call to_history("n3",nx,"sou_out");	call to_history("d3",sog%sh%d,"sou_out");	call to_history("o3",xmin,"sou_out")  
 if (combine) then
   call to_history("n4",ngap,"sou_out");	call to_history("d4",1.,"sou_out");	call to_history("o4",1.,"sou_out")
 else
   call to_history("n4",sog%xm%n,"sou_out");	call to_history("d4",1.,"sou_out");	call to_history("o4",1.,"sou_out")
 end if
 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("wr")) then
 call to_history("n1",rec%w%n,"wr") ; call to_history("d1",rec%w%d,"wr") ; call to_history("o1",rec%w%o,"wr")
 call to_history("n2",img%z%n,"wr"); call to_history("d2",img%z%d,"wr"); call to_history("o2",img%z%o,"wr")
 call to_history("n3",rec%h%n,"wr"); call to_history("d3",rec%h%d,"wr"); call to_history("o3",rec%h%o,"wr")  
 call to_history("n4",1,"wr"); call to_history("d4",1.,"wr"); call to_history("o4",0.,"wr")
 call to_history("n5",1,"wr"); call to_history("d5",1.,"wr"); call to_history("o5",0.,"wr")  
 call to_history("esize",8,"wr")
end if

if(exist_file("velout")) then
 call to_history("n1",sog%z%n,"cigset") ; call to_history("d1",sog%z%d,"cigset") ; call to_history("o1",sog%z%o,"cigset")
 call to_history("n2",nx,"cigset"); call to_history("d2",sog%sh%d,"cigset"); call to_history("o2",xmin,"cigset")
 call to_history("n3",ngap,"cigset"); call to_history("d3",1.,"cigset"); call to_history("o3",1.,"cigset")  
 call to_history("n4",1,"cigset"); call to_history("d4",1.,"cigset"); call to_history("o4",0.,"cigset")
 call to_history("n5",1,"cigset"); call to_history("d5",1.,"cigset"); call to_history("o5",0.,"cigset")  
 call to_history("esize",4,"cigset")
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
!
!kxm%n = 2*padtraces+sog%sh%n
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

write(0,*) "   PARM   "
write(0,*) "=========="
write(0,*) "PARAMETERS"
write(0,*) "=========="
write(0,*) //
write(0,*) "no.ref.velocities   = ",nref
write(0,*) "max.limit velocity  = ",vsal
write(0,*) "traces in taper     = ",tapertr
write(0,*) "padded traces       = ",padtraces
write(0,*) "phase encoding      = ",encode
write(0,*) "no.traces in gap    = ",ngap
write(0,*) //
write(0,*) "============"
write(0,*) "HEADER WORDS"
write(0,*) "============"
write(0,*) //
write(0,*) "Rotated CIG input"
write(0,*) "================="
write(0,*) "no.offsets      = ", sog%sh%n
write(0,*) "offset spacing  = ", sog%sh%d
write(0,*) "min.offset      = ", sog%sh%o
write(0,*) "no.rotated CIGs = ", sog%xm%n
write(0,*) "CIG spacing     = ", sog%xm%d
write(0,*) "min.CIG x-coord = ", sog%xm%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,*) "AREAL SPACE"
write(0,*) "==========="
write(0,*) "no.of traces    = ", nx
write(0,*) "trace spacing   = ", sog%sh%d
write(0,*) "min.trace-x     = ", xmin
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  = ", kxm%n
end if

end subroutine

end module
