!!------------------------------------------------
!!.. This module computes optimized coefficients 
!!.. for FFD migration
!!------------------------------------------------
!!.. Author: Yaxun Tang
!!------------------------------------------------
module wem_ffdopt_mod
  implicit none
contains

  subroutine ffd_coeffts(alph, beta, nstep, cov)

  end subroutine

  !!.. locally optimized coefficients
  subroutine ffd_optloca()
  end subroutine

  !!.. globally optimized coefficients
  subroutine ffd_optglob()
  end subroutine

  !!.. only fit to 4th order Taylor expansion
  subroutine ffd_objloca_kernel(err, derr, alph, beta, p, theta)
    real :: err
    real :: alph, beta
    real :: p, theta

    real :: u, tgt, fit, err

    u = sin(theta)
    tgt = 1./2*u**2+1./8.*sigma(p,2)*u**4+1./16.*sigma(p,3)*u**6+5./128.*sigma(p,4)*u**8
    fit = alph*u*u/(1-beta*u*u)
    err = tgt - fit
    derr =  
  end subroutine


  real function sigma(p, n)
    real    :: p
    integer :: n
    integer :: i
    do i=0, 2*n-2
       sigma = sigma + p**i
    end do
  end function

  integer function factorial(n)
    integer :: n
    integer :: i
    n = 1 
    do i=1, n
       n = n * i
    end do 
  end function


end module
