module fdmigcc_
  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 :: ws(:)

  complex, pointer, private :: nu(:)
contains
  !----------------------------------------------------------------
  subroutine fdmigcc_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(ws(ax%n))

    allocate(nu(ax%n))
  end subroutine fdmigcc_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
  !----------------------------------------------------------------
  !! split-step
  subroutine migccss(dat,img,slo)
    complex, pointer     :: img(:,:)
    real,    pointer     :: slo(:,:)
    complex, pointer     :: dat(:,:)
    real                 :: w,so
    integer              :: iz,iw,ikx
    complex, pointer     :: frq(:)

    real                 :: kz, kx, arg
    complex              :: ikz,cc
    real                 :: dkzds

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

       frq => dat(:,iw)
       do iz= 1, nz
          so = (minval(slo(iz,:))+maxval(slo(iz,:)))/2.

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

             arg = (2*w*so)**2 - kx**2

!!$             if(arg<0.) cycle
!!$             kz = sqrt(arg)
!!$             cc = cexp( cmplx(0.,kz)*(-dz))

             if(arg<0.) then
                ikz=cmplx(sqrt(-arg),0.)
             else
                ikz=cmplx(0.,sqrt(+arg))
             end if
             cc = cexp( ikz*(-dz))

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

          dkzds = 2*w
          frq = frq*cexp( cmplx( 0.,dkzds*(slo(iz,:)-so))*(-dz))
          img(iz,:) = img(iz,:)+frq
       end do
    end do

  end subroutine migccss
  !----------------------------------------------------------------
  subroutine migcc15(dat,img,slo)
    complex, pointer     :: img(:,:)
    real,    pointer     :: slo(:,:)
    complex, pointer     :: dat(:,:)
    real                 :: w
    integer              :: iz,ix,iw
    complex, pointer     :: frq(:)
    complex              :: kn
    complex              :: el,er
    complex, dimension(size(img,2)) :: cd



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

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

       frq => dat(:,iw)
       do iz= 1, nz 
          ws=2*w*slo(iz,:)
          nu = kn/ws

          av =   -   nu +   sixth
          bv = 1.+2.*nu - 2*sixth
          cv =   -   nu +   sixth

          cd( 1) = 0.
          do ix=1+1,nx-1
             cd(ix) = (       nu(ix)+  sixth)*frq(ix-1) &
             +        ( 1.-2.*nu(ix)-2*sixth)*frq(ix  ) &
             +        (       nu(ix)+  sixth)*frq(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, frq)

          frq = frq*cexp( cmplx(0.,ws)*(-dz))
          img(iz,:) = img(iz,:)+frq

       end do
    end do
    call ctris_close()
  end subroutine migcc15
  !----------------------------------------------------------------
  subroutine migcc45(dat,img,slo)
    complex, pointer     :: img(:,:)
    real,    pointer     :: slo(:,:)
    complex, pointer     :: dat(:,:)
    real                 :: w
    integer              :: iz,ix,iw
    complex, pointer     :: frq(:)
    complex              :: kn,kb,bb
    complex              :: el,er
    complex, dimension(size(img,2)) :: cd

    call ctris_init(nx) 

    kn = -(0.,1.)*dz/(4*dx**2)
    kb =           1/(4*dx**2)

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

       frq => dat(:,iw)
       do iz= 1, nz 
          ws=2*w*slo(iz,:)
          nu = kn/ws

          av =   -   nu +    kb/ws**2 +   sixth
          bv = 1.+2.*nu - 2.*kb/ws**2 - 2*sixth
          cv =   -   nu +    kb/ws**2 +   sixth

          cd( 1) = 0.
          do ix=1+1,nx-1
             bb = kb/(ws(ix)**2)
             cd(ix) = (       nu(ix)+   bb+  sixth)*frq(ix-1) &
             +        ( 1.-2.*nu(ix)-2.*bb-2*sixth)*frq(ix  ) &
             +        (       nu(ix)+   bb+  sixth)*frq(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, frq)

          frq = frq*cexp( cmplx(0.,ws)*(-dz))
          img(iz,:) = img(iz,:)+frq
       end do
    end do
    call ctris_close()
  end subroutine migcc45
  !----------------------------------------------------------------
end module fdmigcc_













