module off2ang_fd_mod
  use adj_mod
  use wei_myfft_mod

  implicit none

  complex, dimension(:,:), allocatable, private :: cmod, cdat
  real   , dimension(:)  , allocatable, private :: kz, khx
  integer                             , private :: nkz, nkhx, nax
  real                                , private :: axmin, axmax, dkhx, okhx, oax, dax, dz
  real   , parameter                  , private :: PI = 3.14159265359
contains
  subroutine off2ang_fd_init(dz_in, dhx_in, &
                             nz_in, nhx_in, axmin_in, axmax_in)

    integer  :: nz_in, nhx_in
    real     :: dz_in, dhx_in, axmin_in, axmax_in

    real     :: dhx
    integer  :: ikz, ikhx, nz, nhx

    nz = nz_in
    nhx = nhx_in

    nkz  = nz_in
    nkhx = nhx_in
    nax  = nhx_in
  
    dz  = dz_in
    dhx = dhx_in

    axmin = axmin_in * PI / 180.
    axmax = axmax_in * PI / 180.

    allocate(kz(nkz))
    allocate(khx(nkhx))

    if (nax == 1) then
        dax = 0
    else
        dax = (axmax - axmin) / (nax - 1)
    end if
    oax = axmin

    !! kz
    do ikz=1, nkz/2 + 1
       kz(ikz) = ikz - 1 
    end do
    do ikz=nkz/2+2, nkz
       kz(ikz) = ikz - nkz - 1
    end do 
    kz = kz / (nz * dz)
    kz = 2 * PI * kz 

    !! khx
    do ikhx = 1, nkhx
       khx(ikhx) = - (nkhx - 1) / 2 + ikhx - 1
    end do
    khx = khx / (nhx * dhx)
    khx = 2 * PI * khx
    okhx = khx(1)
    dkhx = khx(2) - khx(1)

!    do ikhx=1, nkhx/2 + 1
!       khx(ikhx) = ikhx - 1
!    end do
!    do ikhx=nkhx/2+2, nkhx
!       khx(ikhx) = ikhx - nkhx - 1
!    end do
!    khx = khx / (nhx * dhx)
!    khx = 2 * PI * khx
!    okhx = 0.
!    dkhx = khx(2) - khx(1)

    write(0,*)"nkz, nkhx, nax", nkz, nkhx, nax
    write(0,*)"dkhx, dax", dkhx, dax
    write(0,*)"okhx, oax", okhx, oax

  end subroutine
!============================================================================
  integer function off2ang_fd_oper(adj, add, modl, data) result(stat)
    logical    :: adj, add
    real, dimension(nkz*nax) :: modl
    real, dimension(nkz*nkhx) :: data
    allocate(cmod(nkz, nax), cdat(nkz, nkhx))
    call adjnull(adj,add,modl,data)
    call off2ang_fd_oper2(adj, modl, data)
    stat = 0 
  end function
!============================================================================
  subroutine off2ang_fd_oper2(adj, modl, data)
    logical  :: adj
    real, dimension(nkz,nax)  :: modl
    real, dimension(nkz,nkhx) :: data
    integer                  :: ikhx, iax, ikz
    real                     :: ax, khx_map

    if (adj) then
       cmod = 0.
       modl = 0.
       cdat = cmplx(data, 0.)
       call FFT2DX(cdat,  1)
    else
    end if


    do ikz=1, nkz
       do iax=1, nax
          ax      = oax + (iax-1)*dax
          khx_map = - tan(ax) * kz(ikz)

          if (abs(khx_map) > 1./abs(tan(ax)*dz)) cycle
          ikhx    = 1.5 + (khx_map - okhx)/dkhx

          if (ikhx > (nkhx-1)/2) then
             ikhx = ikhx - (nkhx-1)/2 
          else
             ikhx = ikhx + nkhx/2 - 1
          end if

          if (ikhx < 1 .or. ikhx > nkhx) cycle
          cmod(ikz, iax) = cmod(ikz, iax) + cdat(ikz, ikhx)

       end do 
    end do

    if (adj) then
       do iax=1, nax
          call FFT1DX(cmod(:, iax), -1)
       end do
       modl = modl + real(cmod)
    else
    end if
  end subroutine
!==============================================================================
  subroutine shuffle(array)
    complex, dimension(:) :: array
    complex, allocatable, dimension(:) :: tmp
    integer               :: i, n
    n = size(array)
    allocate(tmp(n))
    tmp = array
    do i=1, n/2 + 1
       array(i) = tmp((n-1)/2+i)
    end do
    do i=n/2+2, n
       array(i) = tmp(i-(n-1)/2-1)
    end do
  end subroutine
end module