module OW_parms_areal

use sep
use OW_type

implicit none

logical,public, save				:: resmig,resmigtau,resmig_lin,scan,resmigx
type (axis) ,public, save				:: kxm, xm, ym, z, sh, kz, ksh
type (mod_grid)	 ,public, save			:: img
integer	 ,public, save				:: nsecdim,node,nrho
real ,public, save					:: s0,rho,osecdim,dsecdim,pi,drho,orho

contains

subroutine owparam()

integer						:: m,i

call from_param("scan",scan,.false.)
call from_param("resmigx",resmigx,.false.)
call from_param("node",node,9999)
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",nsecdim) ; call from_aux("bkgimg","d2",dsecdim) ; call from_aux("bkgimg","o2",osecdim)
call from_aux("bkgimg","n3",img%z%n) ; call from_aux("bkgimg","d3",img%z%d) ; call from_aux("bkgimg","o3",img%z%o)
call from_aux("bkgimg","n4",img%ym%n); call from_aux("bkgimg","d4",img%ym%d); call from_aux("bkgimg","o4",img%ym%o)  
if (scan.or.resmigx) then
   call from_param("nrho",nrho,101);call from_param("orho",orho,0.9);call from_param("drho",drho,0.002)
   call to_history("n1",img%xm%n,"Dimage"); call to_history("d1",img%xm%d,"Dimage"); call to_history("o1",img%xm%o,"Dimage")  
   call to_history("n2",nsecdim,"Dimage") ; call to_history("d2",dsecdim,"Dimage") ; call to_history("o2",osecdim,"Dimage")
   call to_history("n3",img%z%n,"Dimage") ; call to_history("d3",img%z%d,"Dimage") ; call to_history("o3",img%z%o,"Dimage")
   if (scan) then
      call to_history("n4",nrho,"Dimage"); call to_history("d4",drho,"Dimage"); call to_history("o4",orho,"Dimage")
      call to_history("n5",img%ym%n,"Dimage"); call to_history("d5",img%ym%d,"Dimage"); call to_history("o5",img%ym%o,"Dimage")  
   end if
   if (resmigx) call to_history("n4",1,"Dimage"); call to_history("d4",drho,"Dimage"); call to_history("o4",orho,"Dimage")
   call to_history("n5",1,"Dimage"); call to_history("d5",1,"Dimage"); call to_history("o5",0,"Dimage")  
   call to_history("esize",4,"Dimage")
else
   call from_param("resmig",resmig,.false.)
   call from_param("resmig_lin",resmig_lin,.false.)
   call from_param("resmigtau",resmigtau,.false.)
   if (not(resmig.or.resmig_lin.or.resmigtau.or.resmigx)) call seperr("How to compute residual migration ?")
   call from_param("s0",s0,0.)
   call from_param("rho",rho)
   call to_history("n1",img%xm%n,"Dimage"); call to_history("d1",img%xm%d,"Dimage"); call to_history("o1",img%xm%o,"Dimage")  
   call to_history("n2",nsecdim,"Dimage") ; call to_history("d2",dsecdim,"Dimage") ; call to_history("o2",osecdim,"Dimage")
   call to_history("n3",img%z%n,"Dimage") ; call to_history("d3",img%z%d,"Dimage") ; call to_history("o3",img%z%o,"Dimage")
   call to_history("n4",1,"Dimage"); call to_history("d4",1,"Dimage"); call to_history("o4",0,"Dimage")
   call to_history("n5",1,"Dimage"); call to_history("d5",1,"Dimage"); call to_history("o5",0,"Dimage")  
   call to_history("esize",4,"Dimage")
end if

pi=4*atan(1.)

!--------------------------------------------------
! Wavenumber vectors
!
! BUG DANGER: img%z%n ==> nz !!
!!!kz
   kz%n = img%z%n
   allocate(kz%all(kz%n))
   kz%all=0.
   if(kz%n/=1) then
     m = kz%n/2
     do i=1,m+1
        kz%all(i)=i-1
     end do
     do i=m+2,kz%n
        kz%all(i)=i-kz%n-1
     end do
   else
     kz%all=0. !! make sure k=0 if n=1
   end if
   kz%all = kz%all / (kz%n*img%z%d)      !! 1/length
   kz%all = 2 * PI * kz%all     !! radians
   kz%d = kz%all(2)-kz%all(1)
!!kx

   kxm%n = img%xm%n
   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*img%xm%d)      !! 1/length
   kxm%all = 2 * PI * kxm%all     !! radians
   kxm%d = kxm%all(2)-kxm%all(1)

!   kxm%n = img%xm%n
!   kxm%d = 1./(kxm%n*img%xm%d)
!   kxm%o = 1./(2*img%xm%d)-kxm%d
!   allocate(kxm%all(kxm%n))
!   do i=1,kxm%n
!      kxm%all(i) = mod(kxm%o+(i-1)*kxm%d,1./img%xm%d)-kxm%o
!   end do
!   kxm%all = 2*pi * kxm%all

!!! ksh
   ksh%n = nsecdim
   allocate(ksh%all(nsecdim))
   ksh%all=0.
   if(ksh%n/=1) then
     m = ksh%n/2
     do i=1,m+1
        ksh%all(i)=i-1
     end do
     do i=m+2,ksh%n
        ksh%all(i)=i-ksh%n-1
     end do
   else
     ksh%all=0. !! make sure k=0 if n=1
   end if
   ksh%all = ksh%all / (ksh%n*dsecdim)      !! 1/length
   ksh%all = 2 * PI * ksh%all     !! radians
   if (ksh%n/=1) then
      ksh%d = ksh%all(2)-ksh%all(1)
   else
      ksh%d=1/(nsecdim*dsecdim)
   end if

end subroutine

end module
