module fdmigtt_
  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 fdmigtt_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 fdmigtt_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 migtt(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,*) "migrc15",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 migtt
  !----------------------------------------------------------------
end module fdmigtt_
