module OW_ffd_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal

implicit none

contains

!---------------------------------------------------------
subroutine owffd(iw,inode,sign,ws,vel,tfld,gat) 

character*3					:: gat
integer						:: iw,ix,ikx,i,j
integer						:: inode,sign
integer						:: indref(1)
real						:: sqav(kxm%n)
complex						:: ws(kxm%n,1)
complex						:: tfld(:)
real						:: vel(kxm%n,1),vel1
!ffd variables
real                                                        :: coef1, coef2,coef3
real                                                        :: vw, vw1, vw0
complex                                                     :: url,urr,ror,rol
complex, dimension(kxm%n)                                   :: u,ra,rb,rc 
!ffd variables add
real                                                        :: sqb(kxm%n)
complex                                                     :: sqa(kxm%n),a(kxm%n),b(kxm%n),cc(kxm%n),r(kxm%n)
!ffd variables yx
    real                                                    :: alph(5), beta(5)
    real                                                    :: wov
    integer                                                 :: nstep

 if (gat.eq."sou") then
    tfld=ws(:,1)
    call owfftexes(tfld,inode)
    ws(:,1)=tfld/sqrt(1.*kxm%n)
 else
    tfld=ws(:,1)
    call owfftexer(tfld,inode)
    ws(:,1)=tfld/sqrt(1.*kxm%n)
 end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Phase-shift
!
 vel1=minval(1/vel(:,1)) !reference velocity
 sqav=(rec%w%all(iw) / vel1)**2 - kxm%all**2
 where(sqav >= 0) 
    ws(:,1)=ws(:,1)*exp(cmplx(0,sign*v%z%d*sqrt(sqav)))
 elsewhere
    ws(:,1)=0
 end where
 if (gat.eq."sou") then
    tfld=ws(:,1)
    call owifftexes(tfld,inode)
    ws(:,1)=tfld/sqrt(1.*kxm%n)
 else
    tfld=ws(:,1)
    call owifftexer(tfld,inode)
    ws(:,1)=tfld/sqrt(1.*kxm%n)
 end if

 ws(:,1) = ws(:,1) * exp(cmplx(0,sign*v%z%d*rec%w%all(iw))*(vel(:,1)-1./vel1))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   Fourier FD
!   from yaxun sep134
    call ifd_coeffts(alph, beta, nstep)
    call ffd_correct(ws(:,1), 1/vel(:,1), alph, beta, vel1, rec%w%all(iw), img%xm%d, 1., img%z%d, kxm%n, 1, nstep, sign)

end subroutine
!====================================================================================================
  !!--------------------------------------------------------------------------
  !!.. Fourier finite-difference correction
  !!--------------------------------------------------------------------------
  subroutine ffd_correct(wfld_slice, vmod_slice, alph, beta, vref, ws, dx, dy, dz, nx, ny, nstep, opsign)
    complex                          :: wfld_slice(kxm%n,1)
    real                             :: vmod_slice(kxm%n,1)
    real                             :: vref, ws
    real                             :: dx, dy, dz
    integer                          :: nx, ny, opsign
    real   , dimension(nx+ny)        :: rwrk1d
    complex, dimension(nx*5+ny*5)    :: cwrk1d

    real                             :: alph(5), beta(5)
    integer                          :: ix, iy, istep, nstep 
    integer                          :: lrx1, lry1
    integer                          :: lcx1, lcx2, lcx3, lcx4, lcx5
    integer                          :: lcy1, lcy2, lcy3, lcy4, lcy5

    !!.. asign work space
    lcx1 = 1           !!.. wfdx
    lcx2 = lcx1 + nx    !!.. ahx
    lcx3 = lcx2 + nx    !!.. bhx
    lcx4 = lcx3 + nx    !!.. chx
    lcx5 = lcx4 + nx    !!.. rhx
    lcy1 = lcx5 + nx    !!.. wfdy
    lcy2 = lcy1 + ny    !!.. ahy
    lcy3 = lcy2 + ny    !!.. bhy
    lcy4 = lcy3 + ny    !!.. chy
    lcy5 = lcy4 + ny    !!.. rhy

    lrx1 = 1
    lry1 = lrx1 + nx

    if (nx > 1) then
       do iy=1, ny
          cwrk1d(lcx1:lcx2-1) = wfld_slice(:,iy)
          rwrk1d(lrx1:lry1-1) = vmod_slice(:,iy)
          do istep=1, nstep
             call ffd_sysform(cwrk1d(lcx2), cwrk1d(lcx3), cwrk1d(lcx4), rwrk1d(lrx1), vref, ws, dx, dz, alph(istep), beta(istep), nx, opsign)
             call ifd_dfrterm(cwrk1d(lcx1), cwrk1d(lcx2), cwrk1d(lcx3), cwrk1d(lcx4), cwrk1d(lcx5), nx)
          end do
          wfld_slice(:,iy) = cwrk1d(lcx1:lcx2-1)
       end do
    end if
    if (ny > 1) then
       do ix=1, nx
          cwrk1d(lcy1:lcy2-1) = wfld_slice(ix,:)
          rwrk1d(lry1:nx+ny ) = vmod_slice(ix,:)
          do istep=1, nstep
             call ffd_sysform(cwrk1d(lcy2), cwrk1d(lcy3), cwrk1d(lcy4), rwrk1d(lry1), vref, ws, dy, dz, alph(istep), beta(istep), ny, opsign)
             call ifd_dfrterm(cwrk1d(lcy1), cwrk1d(lcy2), cwrk1d(lcy3), cwrk1d(lcy4), cwrk1d(lcy5), ny)
          end do
          wfld_slice(ix,:) = cwrk1d(lcy1:lcy2-1)
       end do
    end if

    return
  end subroutine


!====================================================================================================
  !!----------------------------------------------------------------------------------------------
  !!.. forming the tridiagonal matrix for FFD extrapolation
  !!----------------------------------------------------------------------------------------------
  subroutine ffd_sysform(ah, bh, ch, velh, vref, ws, dh, dz, alph, beta, nh, opsign)
    complex, dimension(nh) :: ah, bh, ch
    real   , dimension(nh) :: velh
    real                   :: vref, ws, dh, dz
    real                   :: alph, beta
    integer                :: nh, opsign

    integer                :: ih
    real                   :: wov, rov
    complex                :: aah, bbh

    ah = 0; bh = 0; ch = 0

    do ih=1, nh
       wov = ws/velh(ih)
       rov = vref/velh(ih)  
       aah = opsign * alph*cmplx(0.,1.)*dz/2./wov/dh/dh*(1-rov)
       bbh = beta/dh/dh/wov/wov*(rov*rov+rov+1.) + trick

       ah(ih) = bbh - aah
       ch(ih) = bbh - aah
       bh(ih) = 1 - 2.*(bbh - aah)
    end do

    return
  end subroutine

!====================================================================================================
  !!-------------------------------------------------------------------------------------------
  !!.. solving for diffraction term
  !!------------------------------------------------------------------------------------------- 
  subroutine ifd_dfrterm(wfdh, ah, bh, ch, rh, nh)
    complex, dimension(nh) :: wfdh
    complex, dimension(nh) :: ah, bh, ch, rh
    integer                :: nh

    integer                :: ih

    rh = 0.

    do ih=2, nh-1
       rh(ih) = conjg(ah(ih))*wfdh(ih-1) &
              + conjg(bh(ih))*wfdh(ih  ) &
              + conjg(ch(ih))*wfdh(ih+1)
    end do

    !bh(1)  = -ch( 1)
    !bh(nh) = -ah(nh)

    rh(1)  = conjg( bh(1))*wfdh(   1) &
           + conjg( ch(1))*wfdh(   2)

    rh(nh) = conjg(ah(nh))*wfdh(nh-1) &
           + conjg(bh(nh))*wfdh(nh  )

    call ctridag(ah, bh, ch, rh, wfdh, nh)

    return
  end subroutine
!====================================================================================================
  !!---------------------------------------------------------------------------
  !!.. complex tridiagonal solver
  !!---------------------------------------------------------------------------
  subroutine ctridag(a,b,c,r,u,n)
    integer            :: n
    complex            :: a(n), b(n), c(n), r(n), u(n)
    integer            :: j
    integer, parameter :: NMAX = 5000
    complex            :: bet, gam(NMAX)
    if (b(1) == 0.) pause "ctirdag: rewrite equations"
    bet  = b(1)
    u(1) = r(1)/bet
    do j=2, n
       gam(j) = c(j-1)/bet
       bet    = b(j)-a(j)*gam(j)
       if (bet == 0.) pause "ctridag failed"
       u(j) = (r(j)-a(j)*u(j-1))/bet
    end do 
    do j=n-1,1,-1
       u(j) = u(j) - gam(j+1)*u(j+1)
    end do
    return
  end subroutine
!============================================================================================
  !!-----------------------------------------------------------------------------------------
  !!.. set coefficients for implicit finite-difference extrapolator
  !!.. based on Zaitian Ma higher order splitting method
  !!.. Ma, Z., 1981, Finite-difference migration with higher-order approximation:
  !!.. Presented at the 1981 joint meeting of the China Geophysical Society and
  !!.. Soc. of Expl. Geophys., Beijing, China
  !!-----------------------------------------------------------------------------------------
  subroutine ifd_coeffts(alph, beta, nstep)
    real    :: alph(5), beta(5)
    integer :: nstep

    select case (dipop)
      case (15)
           alph(1) = 0.500000000; beta(1) = 0.000000000
           nstep = 1

      case (45)
           alph(1) = 0.500000000; beta(1) = 0.250000000
           nstep = 1

      case (65)
           alph(1) = 0.478242060; beta(1) = 0.376369527
           nstep = 1

      case (79)
           alph(1) = 0.457500000; beta(1) = 0.457500000
           nstep = 1

      case (80)
           alph(1) = 0.457289566; beta(1) = 0.222691983
           alph(2) = 0.040315157; beta(2) = 0.873981642
           nstep = 2

      case (87)
           alph(1) = 0.414236605; beta(1) = 0.150843924
           alph(2) = 0.081312882; beta(2) = 0.744418059
           alph(3) = 0.004210420; beta(3) = 0.972926132
           nstep = 3

      case (89)
           alph(1) = 0.367013245; beta(1) = 0.105756624
           alph(2) = 0.117592008; beta(2) = 0.614520676
           alph(3) = 0.014853510; beta(3) = 0.919432661
           alph(4) = 0.000523275; beta(4) = 0.994065088
           nstep = 4

      case (90)
           alph(1) = 0.318013812; beta(1) = 0.073588213
           alph(2) = 0.143798076; beta(2) = 0.483340757
           alph(3) = 0.033860918; beta(3) = 0.824918565
           alph(4) = 0.004172967; beta(4) = 0.964827992
           alph(5) = 0.000153427; beta(5) = 0.997370236
           nstep = 5

    end select
    return
  end subroutine
!======================================================================================
end module
