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, kym, xm, ym, hx, hy, z, sh, kz, kshx,kshy
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,.true.)
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",img%hx%n); call from_aux("bkgimg","d2",img%hx%d); call from_aux("bkgimg","o2",img%hx%o)
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)  
img%hy%n=1
img%hy%d=img%ym%d
img%hy%o=0.
if (scan.or.resmigx) then
   call from_param("nrho",nrho,101);call from_param("orho",orho,0.8);call from_param("drho",drho,0.004)
   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",img%hx%n,"Dimage"); call to_history("d2",img%hx%d,"Dimage"); call to_history("o2",img%hx%o,"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",img%ym%n,"Dimage"); call to_history("d4",img%ym%d,"Dimage"); call to_history("o4",img%ym%o,"Dimage")  
      call to_history("n5",nrho,"Dimage"); call to_history("d5",drho,"Dimage"); call to_history("o5",orho,"Dimage")
   end if
   if (resmigx) then
     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")  
   end if
   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)

!!ky

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

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

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

end subroutine

end module
