module dimwgh3_
  !! image perturbation weight

  use utilities
  implicit none

  type(paxis) :: laz,lahx,lahy,lamx,lamy

  real,    private :: pi
  real,    private :: tana, tang
  integer, private :: lntap

  real,    dimension(:),         pointer :: kmx,kmy,khx,khy, kz
  real,    dimension(:,:),       pointer :: ksx,ksy,krx,kry
  real,    dimension(:,:,:,:),   pointer :: kss,krr
  real,    dimension(:,:,:,:),   pointer ::  uu,wgh
  real,    dimension(:,:,:,:),   pointer :: tms,tmr
contains
  !----------------------------------------------------------------
  subroutine dimwgh3_init(az,ahx,ahy,amx,amy,amax,gmax,ntap)
    type(paxis)         ::az,ahx,ahy,amx,amy
    real,    intent(in) :: amax,gmax
    integer, intent(in) :: ntap
    integer :: iz,ihx,ihy,imx,imy
    integer :: er

    laz=az
    lahx=ahx
    lahy=ahy
    lamx=amx
    lamy=amy

    lntap = ntap

    pi = acos(-1.)

    tana = abs(tan(pi*amax/180))
    tang = abs(tan(pi*gmax/180))

    allocate(kmx(amx%n),stat=er)
    if(er/=0) call seperr('cant allocate kmx')
    call klinemap(kmx,amx%n,amx%d)
    allocate(kmy(amy%n),stat=er)
    if(er/=0) call seperr('cant allocate kmy')
    call klinemap(kmy,amy%n,amy%d)
    allocate(khx(ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate khx')
    call klinemap(khx,ahx%n,ahx%d)

    allocate(khy(ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate khy')
    call klinemap(khy,ahy%n,ahy%d)

    allocate(kz(az%n),stat=er)
    if(er/=0) call seperr('cant allocate kz')
    call klinemap( kz, az%n, az%d)

    !! ksx(hx,mx), krx(hx,mx)
    allocate(ksx(ahx%n, amx%n),stat=er)
    if(er/=0) call seperr('cant allocate ksx')
    allocate(krx(ahx%n, amx%n),stat=er)
    if(er/=0) call seperr('cant allocate krx')

    !! ksy(hy,my), kry(hy,my)
    allocate(ksy(ahy%n, amy%n),stat=er)
    if(er/=0) call seperr('cant allocate ksy')
    allocate(kry(ahy%n, amy%n),stat=er)
    if(er/=0) call seperr('cant allocate kry')

    !! ks(hx,hy,mx,my)
    allocate(kss(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate kss')

    !! kr(hx,hy,mx,my)
    allocate(krr(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate krr')

    !!  uu(hx,hy,mx,my)
    allocate( uu(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate  uu')

    !! wgh(hx,hy,mx,my)
    allocate(wgh(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate wgh')

    !! tms(hx,hy,mx,my)
    allocate(tms(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate tms')
    !! tmr(hx,hy,mx,my)
    allocate(tmr(ahx%n,ahy%n,amx%n,amy%n),stat=er)
    if(er/=0) call seperr('cant allocate tmr')
    !----------------------------------------------------------------    

    !! ksx = kmx - khx
    !! krx = kmx + khx
    do ihx=1,ahx%n
       do imx=1,amx%n
          ksx(ihx,imx) = kmx(imx)-khx(ihx)
          krx(ihx,imx) = kmx(imx)+khx(ihx)
       end do
    end do

    !! ksy = kmy - khy
    !! kry = kmy + khy
    do ihy=1,ahy%n
       do imy=1,amy%n
          ksy(ihy,imy) = kmy(imy)-khy(ihy)
          kry(ihy,imy) = kmy(imy)+khy(ihy)
       end do
    end do

    !! kss = sqrt(ksx^2 + ksy^2) = |ks|
    !! krr = sqrt(krx^2 + kry^2) = |kr|
    do ihy=1,ahy%n
       do ihx=1,ahx%n
          do imy=1,amy%n
             do imx=1,amx%n
                kss(ihx,ihy,imx,imy) = &
                ksx(ihx,imx)**2      + &
                ksy(ihy,imy)**2

                krr(ihx,ihy,imx,imy) = &
                krx(ihx,imx)**2      + &
                kry(ihy,imy)**2
             end do
          end do
       end do
    end do
    kss = sqrt(kss)
    krr = sqrt(krr)

  end subroutine dimwgh3_init
  !----------------------------------------------------------------  
  subroutine dimwgh3(cdata)
    complex, dimension(:,:,:,:,:), pointer :: cdata
    integer :: iz
    real    :: kz2

    cdata(1,:,:,:,:)=0.0
    do iz=2,laz%n
       kz2=kz(iz)**2

       if(iz<=lntap) then
          cdata(iz,:,:,:,:) = &
          cdata(iz,:,:,:,:) * (1.*       iz /lntap)**3
          write(0,*) "*",iz,  (1.*       iz /lntap)**3
       end if

       if(iz>laz%n-lntap+1) then
          cdata(iz,:,:,:,:) = &
          cdata(iz,:,:,:,:) * (1.*(laz%n-iz)/lntap)**3
          write(0,*) "*",iz,  (1.*(laz%n-iz)/lntap)**3
       end if

       !! tan^2 gamma
       tms = ((krr-kss)/2)**2 / kz2
       where(tms>tang) 
          tms=tang
       end where
       !! tan^2 alpha
       tmr = ((krr+kss)/2)**2 / kz2
       where(tmr>tana) 
          tmr=tana
       end where

       uu = kz2 + kz2*(tms*tmr) + ((krr-kss)/2)**2 + ((krr+kss)/2)**2

       tms=uu - kss**2
       where(tms>1*epsilon(tmr))
          tms=uu/sqrt(tms)
       elsewhere
          tms=0.
       end where

       tmr=uu - krr**2
       where(tmr>1*epsilon(tmr))
          tmr=uu/sqrt(tmr)
       elsewhere 
          tmr=0.
       end where
       
       wgh=0.5*(tms + tmr)

       cdata(iz,:,:,:,:) = &
       cdata(iz,:,:,:,:) * wgh 

    end do
  end subroutine dimwgh3
  !----------------------------------------------------------------
  subroutine klinemap(k,n,d)
    real, dimension(:), pointer :: k
    integer, intent(in) :: n
    real,    intent(in) :: d
    integer             :: i,m

    if(n/=1) then
       m = n/2

       do i=1,m
          k(i  ) =    i-1
          k(i+m) = -m+i-1
       end do

    else
       k=0.
    end if

    k = k / (n*d)  !! 1/length
    k = 2 * pi * k !! radians
  end subroutine klinemap
  !----------------------------------------------------------------
end module dimwgh3_


