module rwekernel_
  use sep
  use fft
  use util_
  use tridiag

  implicit none 
  real, parameter, private :: sixth = 0.122996

  integer, private :: nz,nx,nw,nk
  real,    private :: dz,dx,dw,dk
  real,    private :: oz,ox,ow,ok
  real,    private :: pi

  complex, pointer, private :: av(:),bv(:),cv(:)
  real,    pointer, private :: ko(:,:),b2a(:,:),s2a(:,:),d2a(:,:),c2a(:,:)
  complex, pointer, private :: kda(:,:),kca(:,:)
  complex, pointer, private :: nu(:,:),mu(:,:)

  complex, pointer, private :: dxt(:,:)
  complex, pointer, private :: wfl(:)
contains
  !----------------------------------------------------------------
  subroutine rwemig_init(aw,ax,az)
    type(myaxis), intent(in) :: aw,ax,az

    pi = acos(-1.)

    nw=aw%n
    dw=aw%d * 2*pi
    ow=aw%o * 2*pi

    nx=ax%n
    dx=ax%d
    ox=ax%o

    nz=az%n
    dz=az%d
    oz=az%o

    nk = nx
    ok =   -pi/dx
    dk = 2.*pi/dx/nx

    allocate(av(ax%n),bv(ax%n),cv(ax%n))
    allocate(ko(az%n,ax%n),b2a(az%n,ax%n),s2a(az%n,ax%n),d2a(az%n,ax%n),c2a(az%n,ax%n))
    allocate(kda(az%n,ax%n),kca(az%n,ax%n))
    allocate(nu(az%n,ax%n),mu(az%n,ax%n))

    allocate(dxt(az%n,ax%n))
    allocate(wfl(     ax%n))
  end subroutine rwemig_init
  !----------------------------------------------------------------
  real function getw(iw)
    integer, intent(in) :: iw
    real :: w
    w=ow + dw*(iw-1)
    getw=w
  end function getw

  real function getk(ik)
    integer, intent(in) :: ik
    real :: k
    k=ok + dk*(ik-1)
    getk=k
  end function getk
  !----------------------------------------------------------------
  subroutine rwemig_f15(img,rcp)
    complex, pointer     :: img(:,:)
    real,    pointer     :: rcp(:,:,:)
    real                 :: w
    integer              :: iz,ix,iw
    complex              :: km,kn
    complex              :: el,er
    complex, dimension(size(img,2)) :: cd

    real, pointer :: a(:,:),b(:,:),c(:,:),d(:,:),s(:,:)

    call ctris_init(nx) 
    km = (0.,1.)*(-dz)/(4*dx)
    kn = (0.,1.)*(-dz)/(2*dx**2)

    a => rcp(:,:,1)
    b => rcp(:,:,2)
    c => rcp(:,:,3)
    d => rcp(:,:,4)
    s => rcp(:,:,5)

    b2a=b/(2*a)
    d2a=d/(2*a)
    c2a=c/(2*a)
    s2a=(2*s)**2/a

    do iw=1,nw
       w=getw(iw)
       write(0,*) "rwe-f15",iw,nw,minval(cabs(img)),maxval(cabs(img))

       call sep_read(dxt)
       write(0,*) minval(cabs(dxt)),maxval(cabs(dxt))

       ko = sqrt( s2a * w**2 - b2a**2)

       mu = km/ko *   c2a
       nu =-kn/ko * ((c2a/ko)**2/2 - d2a)

       wfl=(0.,0.)
       do iz= 1, nz 
          wfl(:) = wfl(:) + dxt(iz,:) !! feed-in data

          av =   -  (nu(iz,:)-mu(iz,:)) +   sixth
          bv = 1.+2.*nu(iz,:)           - 2*sixth
          cv =   -  (nu(iz,:)+mu(iz,:)) +   sixth

          cd( 1) = 0.
          do ix=1+1,nx-1
             cd(ix) = (      (nu(iz,ix)-mu(iz,ix))+  sixth)*wfl(ix-1) &
             +        ( 1.-2.*nu(iz,ix)           -2*sixth)*wfl(ix  ) &
             +        (      (nu(iz,ix)+mu(iz,ix))+  sixth)*wfl(ix+1)
          end do
          cd(nx) = 0.

          el = -(av( 1)+bv( 1))/cv( 1)
          er = -(cv(nx)+bv(nx))/av(nx)
          call tris(el, av, bv, cv, er, cd, wfl)

          wfl = wfl*cexp( cmplx(-b2a(iz,:),ko(iz,:))*(-dz))
          img(iz,:) = img(iz,:)+wfl

       end do

       call sseek('out',0,0)
       call sep_write(img)
    end do
    call ctris_close()
  end subroutine rwemig_f15
  !----------------------------------------------------------------
  subroutine rwemig_ssf(img,rcp)
    complex, pointer     :: img(:,:)
    real,    pointer     :: rcp(:,:,:)
    real                 :: w,a0,b0,c0,d0,s0
    integer              :: iz,ix,iw,ikx

    real                 :: kz, kx
    complex              :: arg
    complex              :: ikz,cc
    real, pointer :: a(:,:),b(:,:),c(:,:),d(:,:),s(:,:)
    real                 :: dkzds, dkz
    real                 :: ko0
    complex              :: dkzs,dkzb,dkza

    a => rcp(:,:,1)
    b => rcp(:,:,2)
    c => rcp(:,:,3)
    d => rcp(:,:,4)
    s => rcp(:,:,5)

    do iw=1,nw
       w=getw(iw)
       write(0,*) "rwe-ssf",iw,nw,minval(cabs(img)),maxval(cabs(img))

       call sep_read(dxt)
       write(0,*) minval(cabs(dxt)),maxval(cabs(dxt))

       wfl=(0.,0.)
       do iz= 1, nz 
          a0=(minval(a(iz,:))+maxval(a(iz,:)))/2.
          b0=(minval(b(iz,:))+maxval(b(iz,:)))/2.
          c0=(minval(c(iz,:))+maxval(c(iz,:)))/2.
          d0=(minval(d(iz,:))+maxval(d(iz,:)))/2.
          s0=(minval(s(iz,:))+maxval(s(iz,:)))/2.

          wfl(:) = wfl(:) + dxt(iz,:) !! feed-in data

          !----- FOURIER -----
          call fth(.false.,.false.,wfl)
          do ikx  = 2, nx
             kx=getk(ikx)

             arg = (2*w*s0)**2/a0 - (b0/(2*a0))**2 &
             +     cmplx(0.,c0/a0)*kx      &
             -     cmplx(d0/a0,0.)*kx**2

             ikz = cmplx(-b0/(2*a0),sqrt(arg))
             cc = cexp( ikz*(-dz))

             wfl(ikx) = wfl(ikx) * cc
          end do
          call fth( .true.,.false.,wfl)
          !----- FOURIER -----

          ko0 = sqrt( 4*a0*(2*w*s0)**2 - b0**2 )
          if(   ko0 <epsilon(ko0)) ko0 = ko0 + epsilon(ko0)
          if(abs(a0)<epsilon( a0)) a0  = a0  + epsilon( a0)

          dkzs = (0.,1.)*      2*2*w *(2*w*s0)      /          ko0
          dkzb =-(0.,1.)*  b0                       / (2*a0   *ko0) - 1 /(2*a0   )
          dkza = (0.,1.)*( b0**2-2*a0*(2*w*s0)**2 ) / (2*a0**2*ko0) + b0/(2*a0**2)

          wfl = wfl * cexp(  (   &
          dkzs * (s(iz,:)-s0)  + &
          dkzb * (b(iz,:)-b0)  + &
          dkza * (a(iz,:)-a0)) * (-dz))

          img(iz,:) = img(iz,:)+wfl
       end do

       call sseek('out',0,0)
       call sep_write(img)
    end do
    call ctris_close()
  end subroutine rwemig_ssf
  !----------------------------------------------------------------  
end module rwekernel_
