! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!-----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief 1 æ¬¡å…ƒå‘¨æœŸå¢ƒç•Œ: å®Ÿãƒ•ãƒ¼ãƒªã‚¨å¤‰æ›
!>
!> @details
!> 1 æ¬¡å…ƒå‘¨æœŸå¢ƒç•Œæ¡ä»¶ã®ä¸‹ã§ã®æµä½“é‹å‹•ã‚’
!> å®Ÿãƒ•ãƒ¼ãƒªã‚¨å¤‰æ›ã«ã‚ˆã‚‹ã‚¹ãƒšã‚¯ãƒˆãƒ«æ³•ã§æ•°å€¤è¨ˆç®—ã™ã‚‹ãŸã‚ã®
!> Fortran90 é–¢æ•°ã‚’æä¾›ã™ã‚‹.
!> 2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã® 1 æ¬¡å…ƒã«é–¢ã—ã¦åŒæ™‚ã«ã‚¹ãƒšã‚¯ãƒˆãƒ«è¨ˆç®—ã‚’å®Ÿè¡Œã™ã‚‹ãŸã‚ã®
!> é–¢æ•°ã‚‚æä¾›ã—ã¦ãŠã‚Š, 2, 3 æ¬¡å…ƒé ˜åŸŸã§ã®è¨ˆç®—ã®ãƒ™ãƒ¼ã‚¹ã‚‚æä¾›ã™ã‚‹.
!>
!> ã“ã®ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã¯å†…éƒ¨ã§ ISPACK3/FXPACK ã®ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ã‚’å‘¼ã‚“ã§ã„ã‚‹.
!> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®æ ¼ç´æ–¹æ³•ã«ã¤ã„ã¦ã¯ ISPACK3/FXPACK
!> ã¨ç•°ãªã£ã¦ã„ã‚‹ã®ã§ä»¥ä¸‹ã®ã‚³ãƒ¡ãƒ³ãƒˆã«æ³¨æ„ã•ã‚ŒãŸã„.
!>
!> ## é–¢æ•°ãƒ»å¤‰æ•°ã®åå‰ã¨åž‹ã«ã¤ã„ã¦
!>
!> ### å‘½åæ³•
!>
!> * é–¢æ•°åã®å…ˆé  (`e_`, `g_`, `ae_`, `ag_`) ã¯, è¿”ã™å€¤ã®å½¢ã‚’ç¤ºã—ã¦ã„ã‚‹.
!>   * `e_`  :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿,
!>   * `g_`  :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿,
!>   * `ae_` :: 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ãŒè¤‡æ•°ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿,
!>   * `ag_` :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ãŒè¤‡æ•°ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿.
!>
!> * é–¢æ•°åã®é–“ã®æ–‡å—åˆ—(`Dx`)ã¯, ãã®é–¢æ•°ã®ä½œç”¨ã‚’è¡¨ã—ã¦ã„ã‚‹.
!>
!> * é–¢æ•°åã®æœ€å¾Œ (`_e`, `_ae`, `_g`, `_ag`) ã¯,
!>   å…¥åŠ›å¤‰æ•°ã®å½¢ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ãŠã‚ˆã³æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã§ã‚ã‚‹ã“ã¨ã‚’ç¤ºã—ã¦ã„ã‚‹.
!>   * `_e`  :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
!>   * `_g`  :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
!>   * `_ae` :: 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ãŒè¤‡æ•°ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿
!>   * `_ag` :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ãŒè¤‡æ•°ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿
!>
!> ### å„ãƒ‡ãƒ¼ã‚¿ã®ç¨®é¡žã®èª¬æ˜Ž
!>
!> * ä»¥ä¸‹ `DP  = kind(1.0D)` ã§ã‚ã‚‹
!>
!> * `g` : 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿.
!>   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ `real(DP), dimension(0:im-1)`.
!>   * `im` ã¯ X åº§æ¨™ã®æ ¼åç‚¹æ•°ã§ã‚ã‚Š, ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ ae_initial() ã«ã¦
!>     ã‚ã‚‰ã‹ã˜ã‚è¨å®šã—ã¦ãŠã.
!>
!> * `e` : ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿.
!>   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ `real(DP), dimension(-km:km)`.
!>   * `km` ã¯ X æ–¹å‘ã®æœ€å¤§æ³¢æ•°ã§ã‚ã‚Š, ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ ae_initial() ã«ã¦
!>     ã‚ã‚‰ã‹ã˜ã‚è¨å®šã—ã¦ãŠã.
!>   * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®æ ¼ç´ã®ã•ã‚Œæ–¹ã¯
!>     `0:km` ã¾ã§ãŒ \f$\cos(kx)\f$ ã®ä¿‚æ•°,
!>     `-km:-1` ãŒ  \f$sin(kx)\f$ ã®ä¿‚æ•°ã¨ãªã£ã¦ã„ã‚‹.
!>
!> * `ag` : 1 æ¬¡å…ƒ(X)æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿.
!>   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ `real(DP), dimension(:,0:im-1)`.
!>     ç¬¬ 2 æ¬¡å…ƒãŒ X æ–¹å‘ã‚’è¡¨ã™.
!>
!> * `ae` : 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿.
!>   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ `real(DP), dimension(:,-km:km)`.
!>     ç¬¬ 2 æ¬¡å…ƒãŒã‚¹ãƒšã‚¯ãƒˆãƒ«ã‚’è¡¨ã™.
!>
!> * `g_` ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
!>
!> * `e_` ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
!>
!> * `ag_` ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã
!>   2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
!>
!> * `ae_` ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã
!>   2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
!>
!> * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã«å¯¾ã™ã‚‹å¾®åˆ†ç‰ã®ä½œç”¨ã¨ã¯, å¯¾å¿œã™ã‚‹æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã«
!>   å¾®åˆ†ãªã©ã‚’ä½œç”¨ã•ã›ãŸãƒ‡ãƒ¼ã‚¿ã‚’ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã—ãŸã‚‚ã®ã®ã“ã¨ã§ã‚ã‚‹.
!>
module ae_module
  use dc_message, only: MessageNotify
  use dc_types, only: DP
  use spml_utils, only: pi
  implicit none

  private
  public ae_Initial
  public ag_ae
  public ae_ag
  public g_e
  public e_g
  public ae_Dx_ae
  public e_Dx_e
  public a_Int_ag
  public Int_g
  public a_Avr_ag
  public Avr_g
  public Interpolate_e
  public a_Interpolate_ae
  public g_X
  public g_X_Weight

  !> æ ¼åç‚¹åº§æ¨™(X)ã‚’æ ¼ç´ã—ãŸ1æ¬¡å…ƒé…åˆ—.
  real(DP), allocatable    :: g_x(:)
  !> é‡ã¿åº§æ¨™ã‚’æ ¼ç´ã—ãŸ 1 æ¬¡å…ƒé…åˆ—. Xæ–¹å‘ã®æ ¼åç‚¹ã®é–“éš”ãŒæ ¼ç´ã—ã¦ã‚ã‚‹.
  real(DP), allocatable    :: g_x_weight(:)

  integer(8)               :: im = 32    ! æ ¼åç‚¹ã®æ•°
  integer(8)               :: km = 10    ! åˆ‡æ–æ³¢æ•°
  real(DP)                 :: xl = 1.0D0 ! é ˜åŸŸã®å¤§ãã•
  integer(8), allocatable  :: it(:)
  real(DP),   allocatable  :: t(:)

  save im, km, it, t, xl, g_X, g_X_Weight

contains
  !-----------------------------------------------------------------------
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã®æ ¼åç‚¹æ•°, æ³¢æ•°, é ˜åŸŸã®å¤§ãã•ã‚’è¨å®šã™ã‚‹.
  !>
  !> ä»–ã®é–¢æ•°ã‚„å¤‰æ•°ã‚’å‘¼ã¶å‰ã«, æœ€åˆã«ã“ã®ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ã‚’å‘¼ã‚“ã§
  !> åˆæœŸè¨å®šã‚’ã—ãªã‘ã‚Œã°ãªã‚‰ãªã„.
  !>
  subroutine ae_Initial(i,k,xmin,xmax)
    !> æ ¼åç‚¹ã®æ•°
    integer, intent(in) :: i
    !> åˆ‡æ–æ³¢æ•°
    integer, intent(in) :: k
    !> Xåº§æ¨™ä¸‹ç«¯
    real(DP),intent(in) :: xmin
    !> Xåº§æ¨™ä¸Šç«¯
    real(DP),intent(in) :: xmax
    integer(8)          :: ii

    im = i
    km = k
    xl = xmax - xmin

    if ( im <= 0 .or. km <= 0 ) then
       call MessageNotify('E','ae_Initial', &
            'Number of grid points and waves should be positive')
    elseif ( mod(im,2_8) /= 0 ) then
       call MessageNotify('E','ae_Initial', &
            'Number of grid points should be even')
    elseif ( km > im/2 ) then
       call MessageNotify('E','ae_Initial', &
            'KM shoud be less equal to IM/2')
    endif

    allocate(it(im/2),t(im*3/2))

    call fxrini(im,it,t)

    allocate(g_x(0:im-1))
    do ii=0,im-1
       g_X(ii) = xmin + (xl/im)*ii
    enddo

    allocate(g_x_weight(0:im-1))
    g_X_Weight = xl/im

    call MessageNotify(&
         'M','ae_initial','ae_module (2019/09/17) is initialized')

  end subroutine ae_Initial

  !-----------------------------------------------------------------------
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã¸é€†å¤‰æ›ã™ã‚‹(2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ç”¨)
  !>
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«é€†å¤‰æ›ã®å®šç¾©ã¯ä»¥ä¸‹ã®ã¨ãŠã‚Š:
  !> * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿: \f$e_k \, (k=-km,...,km)\f$
  !> * æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿:     \f$g_j \, (j=-0,...,im-1)\f$
  !> \f{align*}{
  !>  g_j = e_0
  !>      + 2\sum_{k=1}^{km}e_k\cos(2\pi jk/im) - e_{-k}\sin(2\pi jk/im))
  !>                                               \quad (j=0,1,...,im-1).
  !> \f}
  function ag_ae(ae)
    !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,-km:), intent(in)  :: ae
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(size(ae,1),0:im-1)   :: ag_ae

    integer(8) :: nm, k
    nm=size(ae,1)

    if ( size(ae,2) < 2*km+1 ) then
       call MessageNotify('E','ag_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call MessageNotify('W','ag_ae', &
            'The Fourier dimension of input data too large.')
    endif

    ag_ae       = 0.0D0
    ag_ae(:,0)  = ae(:,0)
    do k=1,min(km,im/2-1)
       ag_ae(:,2*k)   = ae(:,k)
       ag_ae(:,2*k+1) = ae(:,-k)
    enddo
    if ( km == im/2 ) then
      ag_ae(:,1) = ae(:,km)
    end if

    call fxrtba(nm,im,ag_ae,it,t)

  end function ag_ae

  !-----------------------------------------------------------------------
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã¸é€†å¤‰æ›ã™ã‚‹(1 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ç”¨)
  !>
  !>  * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿: \f$e_k \, (k=-km,...,km)\f$
  !>  * æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿:  \f$g_j \, (j=-0,...,im-1)\f$
  !> \f{align*}{
  !>  g_j = e_0
  !>      + 2\sum_{k=1}^{km}(e_k\cos(2\pi jk/im) - e_{-k}\sin(2\pi jk/im))
  !>                                                  \quad (j=0,1,...,im-1).
  !> \f}
  function g_e(e)
    !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(-km:km), intent(in) :: e
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(0:im-1)             :: g_e

    real(DP), dimension(1,size(e))  :: ae_work
    real(DP), dimension(1,0:im-1)   :: ag_work

    ae_work(1,:) = e
    ag_work = ag_ae(ae_work)
    g_e = ag_work(1,:)

  end function g_e

  !-----------------------------------------------------------------------
  !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã¸æ£å¤‰æ›ã™ã‚‹(2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ç”¨)
  !>
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«æ£å¤‰æ›ã®å®šç¾©ã¯ä»¥ä¸‹ã®ã¨ãŠã‚Š:
  !> * æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿     \f$g_j \, (j=-0,...,im-1)\f$
  !> * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ \f$e_k \, (k=-km,...,km)\f$
  !> \f{align*}{
  !>   e_0 & = (1/im)\sum_{j=0}^{im-1} g_j \\
  !>   e_k & = (1/im)\sum_{j=0}^{im-1} g_j \cos(2\pi jk/im)
  !>          \quad (k=1,2,...,km) \\
  !>   e_{-k} & = - (1/im)\sum_{j=0}^{im-1} g_j \sin(2\pi jk/im)
  !>          \quad (k=1,2,...,km)
  !> \f}
  function ae_ag(ag)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,:), intent(in)   :: ag
    !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(size(ag,1),-km:km) :: ae_ag

    real(DP), dimension(size(ag,1),0:im-1) :: ag_work
    integer(8) :: nm, k

    nm = size(ag,1)
    if ( size(ag,2) < im ) then
       call MessageNotify('E','ae_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','ae_ag', &
            'The Grid points of input data too large.')
    endif
    ag_work = ag

    call fxrtfa(nm,im,ag_work,it,t)

    do k=1,min(km,im/2-1)
       ae_ag(:, k) = ag_work(:,2*k)
       ae_ag(:,-k) = ag_work(:,2*k+1)
    enddo
    ae_ag(:,0) = ag_work(:,0)
    if ( km == im/2 ) then
      ae_ag(:,km) = ag_work(:,1)
    end if

  end function ae_ag

  !-----------------------------------------------------------------------
  !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã¸æ£å¤‰æ›ã™ã‚‹(1 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ç”¨)
  !>
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«æ£å¤‰æ›ã®å®šç¾©ã¯ä»¥ä¸‹ã®ã¨ãŠã‚Š:
  !> * æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿:     \f$g_j \, (j=-0,...,im-1)\f$
  !> * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿: \f$e_k \, (k=-km,...,km)\f$
  !> \f{align*}{
  !>  e_0    & = (1/im)\sum_{j=0}^{im-1} g_j \\
  !>  e_k    & = (1/im)\sum_{j=0}^{im-1} g_j \cos(2\pi jk/im)
  !>             \quad (k=1,2,...,km) \\
  !>  e_{-k} & = - (1/im)\sum_{j=0}^{im-1} g_j \sin(2\pi jk/im)
  !>              \quad (k=1,2,...,km)
  !> \f}
  function e_g(g)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(0:im-1), intent(in)  :: g
    !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(-km:km)              :: e_g

    real(DP), dimension(1,size(g))        :: ag_work
    real(DP), dimension(1,-km:km)         :: ae_work

    ag_work(1,:) = g
    ae_work = ae_ag(ag_work)
    e_g = ae_work(1,:)

  end function e_g

  !-----------------------------------------------------------------------
  !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã« X å¾®åˆ†ã‚’ä½œç”¨ã™ã‚‹(2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿).
  !>
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã® X å¾®åˆ†ã¨ã¯, å¯¾å¿œã™ã‚‹æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã« X å¾®åˆ†ã‚’
  !> ä½œç”¨ã•ã›ãŸãƒ‡ãƒ¼ã‚¿ã®ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã®ã“ã¨ã§ã‚ã‚‹.
  function ae_Dx_ae(ae)
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,-km:), intent(in)  :: ae
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã® X å¾®åˆ†
    real(DP), dimension(size(ae,1),-km:km)   :: ae_dx_ae

    integer(8) ::  k

    if ( size(ae,2) < 2*km+1 ) then
       call MessageNotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call MessageNotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too large.')
    endif

    do k=-km,km
       ae_Dx_ae(:,k) = -(2*pi*k/xl)*ae(:,-k)
    enddo
  end function ae_dx_ae

  !-----------------------------------------------------------------------
  !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã« X å¾®åˆ†ã‚’ä½œç”¨ã™ã‚‹(1 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿).
  !>
  !> ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã® X å¾®åˆ†ã¨ã¯, å¯¾å¿œã™ã‚‹æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã« X å¾®åˆ†ã‚’
  !> ä½œç”¨ã•ã›ãŸãƒ‡ãƒ¼ã‚¿ã®ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã®ã“ã¨ã§ã‚ã‚‹.
  function e_Dx_e(e)
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(-km:km), intent(in)     :: e
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã® X å¾®åˆ†
    real(DP), dimension(-km:km)                 :: e_Dx_e

    real(DP), dimension(1,-km:km)               :: ae_work

    ae_work(1,:) = e
    ae_work = ae_Dx_ae(ae_work)
    e_Dx_e = ae_work(1,:)

  end function e_Dx_e

  !-----------------------------------------------------------------------
  !> 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ãŒä¸¦ã‚“ã  2 æ¬¡å…ƒé…åˆ—ã®ç©åˆ†
  function a_Int_ag(ag)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,0:), intent(in)     :: ag
    !> ç©åˆ†çµæžœ
    real(DP), dimension(size(ag,1))           :: a_Int_ag

    integer(8) :: i

    if ( size(ag,2) < im ) then
       call MessageNotify('E','ae_Int_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','ae_Int_ag', &
            'The Grid points of input data too large.')
    endif

    a_Int_ag = 0.0D0
    do i=0,im-1
       a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i)
    enddo
  end function a_Int_ag

  !-----------------------------------------------------------------------
  !> 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ç©åˆ†
  function Int_g(g)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(0:im-1), intent(in)   :: g
    !> ç©åˆ†çµæžœ
    real(DP)                                  :: Int_g

    Int_g = sum(g*g_X_Weight)

  end function Int_g

  !-----------------------------------------------------------------------
  !> 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ãŒä¸¦ã‚“ã  2 æ¬¡å…ƒé…åˆ—ã®å¹³å‡
  function a_Avr_ag(ag)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,0:), intent(in)     :: ag
    !> å¹³å‡çµæžœ
    real(DP), dimension(size(ag,1))           :: a_Avr_ag

    a_Avr_ag = a_Int_ag(ag)/sum(g_X_Weight)

  end function a_Avr_ag

  !-----------------------------------------------------------------------
  !> 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®å¹³å‡
  function Avr_g(g)
    !> æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(0:im-1), intent(in)   :: g
    !> å¹³å‡å€¤
    real(DP)                                  :: Avr_g

    Avr_g = Int_g(g)/sum(g_X_Weight)

  end function Avr_g

  !-----------------------------------------------------------------------
  !> è£œé–“è¨ˆç®—(1æ¬¡å…ƒ)
  function Interpolate_e(e_data,xval)
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(-km:km), intent(IN) :: e_data
    !> è£œé–“ã™ã‚‹ç‚¹ã®åº§æ¨™
    real(DP), intent(in)                    :: xval
    !> è£œé–“ã—ãŸçµæžœã®å€¤
    real(DP)                                :: Interpolate_e

    integer(8) :: k

    Interpolate_e = e_data(0)
    do k=1,km
       Interpolate_e = Interpolate_e                             &
            &  + 2.0D0 * (   e_data(k)  * cos(2*PI*k/xl*xval)    &
            &              - e_data(-k) * sin(2*PI*k/xl*xval) )
    end do

  end function Interpolate_e

  !-----------------------------------------------------------------------
  !> è£œé–“è¨ˆç®—(2æ¬¡å…ƒ)
  function a_Interpolate_ae(ae_data,xval)
    !> å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
    real(DP), dimension(:,-km:), intent(IN)   :: ae_data
    !> è£œé–“ã™ã‚‹ç‚¹ã®åº§æ¨™
    real(DP), intent(in)                      :: xval
    !> è£œé–“ã—ãŸçµæžœã®å€¤
    real(DP), dimension(size(ae_data,1))      :: a_Interpolate_ae

    integer(8) :: k

    a_Interpolate_ae = ae_data(:,0)

    do k=1,km
       a_Interpolate_ae = a_Interpolate_ae                          &
            &  + 2.0D0 * (   ae_data(:,k)  * cos(2*PI*k/xl*xval)    &
            &              - ae_data(:,-k) * sin(2*PI*k/xl*xval) )
    end do

  end function a_Interpolate_ae

end module ae_module