!--
!----------------------------------------------------------------------
!     Copyright (c) 2009-2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!è¡¨é¡Œ  ae_module_fftj
!
!      spml/ae_module_fftj ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã¯ 1 æ¬¡å…ƒå‘¨æœŸå¢ƒç•Œæ¡ä»¶ã®ä¸‹ã§ã®æµä½“é‹å‹•ã‚’
!      å®Ÿãƒ•ãƒ¼ãƒªã‚¨å¤‰æ›ã«ã‚ˆã‚‹ã‚¹ãƒšã‚¯ãƒˆãƒ«æ³•ã§æ•°å€¤è¨ˆç®—ã™ã‚‹ãŸã‚ã® Fortran90 é–¢æ•°
!      ã‚’æä¾›ã™ã‚‹.
!
!      2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã® 1 æ¬¡å…ƒã«é–¢ã—ã¦åŒæ™‚ã«ã‚¹ãƒšã‚¯ãƒˆãƒ«è¨ˆç®—ã‚’å®Ÿè¡Œã™ã‚‹ãŸã‚ã®
!      é–¢æ•°ã‚‚æä¾›ã—ã¦ãŠã‚Š, 2, 3 æ¬¡å…ƒé ˜åŸŸã§ã®è¨ˆç®—ã®ãƒ™ãƒ¼ã‚¹ã‚‚æä¾›ã™ã‚‹.
!
!      ã“ã®ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã¯å†…éƒ¨ã§ fftj ã® Fortran77 ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ã‚’
!      å‘¼ã‚“ã§ã„ã‚‹.
!
!      æ£å¤‰æ›ã®å®šç¾©ãŠã‚ˆã³ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®æ ¼ç´æ–¹æ³•ã«ã¤ã„ã¦ã¯ fftj ã¨
!      ç•°ãªã£ã¦ã„ã‚‹ã®ã§ä»¥ä¸‹ã®ã‚³ãƒ¡ãƒ³ãƒˆã«æ³¨æ„ã•ã‚ŒãŸã„.
!
!
!å±¥æ´  2009/09/09  ç«¹åºƒçœŸä¸€  fftj ã® Fortran 90 ç‰ˆ. ae_module.f90 ã‚ˆã‚Šæ”¹é€
!      2011/12/03  ç«¹åºƒçœŸä¸€  è£œé–“é–¢æ•°ã‚’è¿½åŠ
!      2013/05/22  ä½ã€…æœ¨æ´‹å¹³ use dc_types DP instead of real(8)
!      2016/08/28  ç«¹åºƒçœŸä¸€  OpenMP å¯¾å¿œ
!
!++
module ae_module_fftj
  !
  != ae_module_fftj
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Copyright :: 2009-2012 SPMODEL Development Group
  ! License   :: MIT/X11. Please see COPYRIGHT[link:../COPYRIGHT]
  !
  !== æ¦‚è¦
  !
  !
  ! spml/ae_module_fftj ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã¯ 1 æ¬¡å…ƒå‘¨æœŸå¢ƒç•Œæ¡ä»¶ã®ä¸‹ã§ã®æµä½“é‹å‹•ã‚’
  ! å®Ÿãƒ•ãƒ¼ãƒªã‚¨å¤‰æ›ã«ã‚ˆã‚‹ã‚¹ãƒšã‚¯ãƒˆãƒ«æ³•ã§æ•°å€¤è¨ˆç®—ã™ã‚‹ãŸã‚ã® Fortran90 é–¢æ•°
  ! ã‚’æä¾›ã™ã‚‹.
  !
  ! 2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã® 1 æ¬¡å…ƒã«é–¢ã—ã¦åŒæ™‚ã«ã‚¹ãƒšã‚¯ãƒˆãƒ«è¨ˆç®—ã‚’å®Ÿè¡Œã™ã‚‹ãŸã‚ã®
  ! é–¢æ•°ã‚‚æä¾›ã—ã¦ãŠã‚Š, 2, 3 æ¬¡å…ƒé ˜åŸŸã§ã®è¨ˆç®—ã®ãƒ™ãƒ¼ã‚¹ã‚‚æä¾›ã™ã‚‹.
  !
  ! ã“ã®ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã¯å†…éƒ¨ã§ fftj ã® Fortran77 ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ã‚’
  ! å‘¼ã‚“ã§ã„ã‚‹.
  !
  ! æ£å¤‰æ›ã®å®šç¾©ãŠã‚ˆã³ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®æ ¼ç´æ–¹æ³•ã«ã¤ã„ã¦ã¯ fftj ã¨
  ! ç•°ãªã£ã¦ã„ã‚‹ã®ã§ä»¥ä¸‹ã®ã‚³ãƒ¡ãƒ³ãƒˆã«æ³¨æ„ã•ã‚ŒãŸã„.
  !
  !
  !== é–¢æ•°ãƒ»å¤‰æ•°ã®åå‰ã¨åž‹ã«ã¤ã„ã¦
  !
  !=== å‘½åæ³•
  !
  ! * é–¢æ•°åã®å…ˆé  (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 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿
  !
  !=== å„ãƒ‡ãƒ¼ã‚¿ã®ç¨®é¡žã®èª¬æ˜Ž
  !
  ! * g : 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿.
  !   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ real(8), dimension(0:im-1).
  !   * im ã¯ X åº§æ¨™ã®æ ¼åç‚¹æ•°ã§ã‚ã‚Š, ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ ae_Initial ã«ã¦
  !     ã‚ã‚‰ã‹ã˜ã‚è¨å®šã—ã¦ãŠã.
  !
  ! * e : ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿.
  !   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ real(8), dimension(-km:km).
  !   * km ã¯ X æ–¹å‘ã®æœ€å¤§æ³¢æ•°ã§ã‚ã‚Š, ã‚µãƒ–ãƒ«ãƒ¼ãƒãƒ³ ae_Initial ã«ã¦
  !     ã‚ã‚‰ã‹ã˜ã‚è¨å®šã—ã¦ãŠã.
  !   * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®æ ¼ç´ã®ã•ã‚Œæ–¹ã¯ 0:km ã¾ã§ãŒ cos(kx) ã®ä¿‚æ•°,
  !     -km:-1 ãŒ sin(kx) ã®ä¿‚æ•°ã¨ãªã£ã¦ã„ã‚‹.
  !
  ! * ag : 1 æ¬¡å…ƒ(X)æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿.
  !   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ real(8), dimension(:,0:im-1).
  !     ç¬¬ 2 æ¬¡å…ƒãŒ X æ–¹å‘ã‚’è¡¨ã™.
  !
  ! * ae : 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã  2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿.
  !   * å¤‰æ•°ã®ç¨®é¡žã¨æ¬¡å…ƒã¯ real(8), dimension(:,-km:km).
  !     ç¬¬ 2 æ¬¡å…ƒãŒã‚¹ãƒšã‚¯ãƒˆãƒ«ã‚’è¡¨ã™.
  !
  ! * g_ ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
  !
  ! * e_ ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
  !
  ! * ag_ ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã
  !   2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
  !
  ! * ae_ ã§å§‹ã¾ã‚‹é–¢æ•°ãŒè¿”ã™å€¤ã¯ 1 æ¬¡å…ƒã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã
  !   2 æ¬¡å…ƒãƒ‡ãƒ¼ã‚¿ã«åŒã˜.
  !
  ! * ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã«å¯¾ã™ã‚‹å¾®åˆ†ç‰ã®ä½œç”¨ã¨ã¯, å¯¾å¿œã™ã‚‹æ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã«
  !   å¾®åˆ†ãªã©ã‚’ä½œç”¨ã•ã›ãŸãƒ‡ãƒ¼ã‚¿ã‚’ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã—ãŸã‚‚ã®ã®ã“ã¨ã§ã‚ã‚‹.
  !
  !== å¤‰æ•°ãƒ»æ‰‹ç¶šãç¾¤ã®è¦ç´„
  !
  !==== åˆæœŸåŒ–
  !
  ! ae_Initial       :: ã‚¹ãƒšã‚¯ãƒˆãƒ«å¤‰æ›ã®æ ¼åç‚¹æ•°, æ³¢æ•°, é ˜åŸŸã®å¤§ãã•ã®è¨å®š
  !
  !==== åº§æ¨™å¤‰æ•°
  !
  ! g_X              :: æ ¼åç‚¹åº§æ¨™(X)ã‚’æ ¼ç´ã—ãŸ 1 æ¬¡å…ƒé…åˆ—
  ! g_X_Weight       :: é‡ã¿åº§æ¨™ã‚’æ ¼ç´ã—ãŸ 1 æ¬¡å…ƒé…åˆ—
  !
  !==== åŸºæœ¬å¤‰æ›
  !
  ! g_e, ag_ae       :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰æ ¼åãƒ‡ãƒ¼ã‚¿ã¸ã®å¤‰æ›
  ! e_g, ae_ag       :: æ ¼åãƒ‡ãƒ¼ã‚¿ã‹ã‚‰ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã¸ã®å¤‰æ›
  !
  !==== å¾®åˆ†
  !
  ! e_Dx_e, ae_Dx_ae :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã« X å¾®åˆ†ã‚’ä½œç”¨ã•ã›ã‚‹
  !
  !==== ç©åˆ†ãƒ»å¹³å‡
  !
  ! a_Int_ag, a_Avr_ag :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ä¸¦ã‚“ã  2 æ¬¡å…ƒé…åˆ—ã®ç©åˆ†ãŠã‚ˆã³å¹³å‡
  ! Int_g, Avr_g       :: 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®ç©åˆ†ãŠã‚ˆã³å¹³å‡
  !
  !==== è£œé–“
  !
  ! Interpolate_e      :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰ã®ç‰¹å®šã®åº§æ¨™ã®å€¤ã®è£œé–“(1 æ¬¡å…ƒ)
  ! a_Interpolate_ae   :: ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿ã‹ã‚‰ã®ç‰¹å®šã®åº§æ¨™ã®å€¤ã®è£œé–“(2 æ¬¡å…ƒ)
  !
  !
  use dc_message
  use dc_types, only: DP
  implicit none

  private
  public ae_Initial                       ! åˆæœŸåŒ–ãƒ«ãƒ¼ãƒãƒ³
  public ag_ae, ae_ag, g_e, e_g           ! åŸºæœ¬å¤‰æ›
  public ae_Dx_ae, e_Dx_e                 ! å¾®åˆ†
  public a_Int_ag, Int_g, a_Avr_ag, Avr_g ! ç©åˆ†ãƒ»å¹³å‡
  public Interpolate_e, a_Interpolate_ae  ! è£œé–“
  public g_X, g_X_Weight                  ! åº§æ¨™å¤‰æ•°

  integer, parameter :: index2max=11      ! Max(im) = 2048 = 2^{11}
  integer            :: im=32             ! æ ¼åç‚¹ã®æ•°
  integer            :: km=10             ! åˆ‡æ–æ³¢æ•°
  real(DP)           :: xl=1.0            ! é ˜åŸŸã®å¤§ãã•

  real(DP), allocatable :: xtb(:)
  integer               :: ipb(2)
  real(DP), allocatable :: xtf(:)
  integer               :: ipf(2)

  real(DP), parameter                 :: pi=3.1415926535897932385D0

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

  save im, km, xl, xtb, ipb, xtf, ipf, g_X, g_X_Weight

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

      integer :: ii, index2

      im = i
      km = k
      xl = xmax-xmin

      index2 = 1
      do while ( im /= 2**index2 .AND. index2 <= index2max )
         index2=index2+1
      enddo

      if ( index2 > index2max ) then
         call MessageNotify('E','ae_Initial', &
              'Number of grid points should be 2^N and <= 2048')
      else if ( km <= 0 ) then
         call MessageNotify('E','ae_Initial', &
              'Number of waves should be positive')
      elseif ( km >= im/2 ) then
         call MessageNotify('E','ae_Initial', &
              'KM shoud be less than IM/2')
      endif

      allocate(xtb(im*3),xtf(im*3))
      call fjrini(im,2,1,ipb,xtb)
      call fjrini(im,2,2,ipf,xtf)

      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_fftj (2016/08/28) is initialized')

    end subroutine ae_Initial

  !--------------- åŸºæœ¬å¤‰æ› -----------------

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

      real(DP), dimension(size(ae,1),0:im-1)   :: ag_in  ! å…¥åŠ›ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿
      real(DP), dimension(im)                  :: g_work
      integer :: m, k

      if ( im > 2048 ) then
         call MessageNotify('E','ag_ae', &
              'The dimension of output data exceeds 2048.')
      else if ( mod(im,2) /= 0 ) then
         call MessageNotify('E','ag_ae', &
              'The dimension of output data should be 2^N.')
      else 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_in = 0.0D0
      ag_in(:,0)=ae(:,0)
      ag_in(:,1)=0.0D0
      do k=1,km
         ag_in(:,2*k)   = ae(:,k)
         ag_in(:,2*k+1) = ae(:,-k)
      enddo
      ag_in(:,2*km+2:im-1)=0.0D0

!$omp parallel do default(shared) private(g_work)
      do m=1,size(ae,1)
         call fjrrun(ag_in(m,:),ag_ae(m,:),g_work,xtb,ipb)
      enddo
!$omp end parallel do

    end function ag_ae

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

      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

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

      real(DP), dimension(size(ag,1),0:im-1)   :: ag_out ! ã‚¹ãƒšã‚¯ãƒˆãƒ«ãƒ‡ãƒ¼ã‚¿

      real(DP), dimension(im)                  :: g_work
      integer :: m, k

      if ( im > 2048 ) then
         call MessageNotify('E','ag_ae', &
              'The dimension of output data exceeds 2048.')
      else if ( mod(im,2) /= 0 ) then
         call MessageNotify('E','ag_ae', &
              'The dimension of output data should be 2^N.')
      else 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

!$omp parallel do default(shared) private(g_work)
      do m=1,size(ag,1)
         call fjrrun(ag(m,:),ag_out(m,:),g_work,xtf,ipf)
      enddo
!$omp end parallel do

      do k=1,km
         ae_ag(:,k)  = ag_out(:,2*k)
         ae_ag(:,-k) = ag_out(:,2*k+1)
      enddo
      ae_ag(:,0) = ag_out(:,0)

      ae_ag = ae_ag/im

    end function ae_ag

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

      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

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

      integer 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

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

      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

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

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

      Int_g = sum(g*g_X_Weight)

    end function Int_g

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

      a_Avr_ag = a_Int_ag(ag)/sum(g_X_Weight)

    end function a_Avr_ag

    function Avr_g(g)
      !
      ! 1 æ¬¡å…ƒæ ¼åç‚¹ãƒ‡ãƒ¼ã‚¿ã®å¹³å‡
      !
      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

  !--------------- è£œé–“è¨ˆç®— -----------------
    function Interpolate_e(e_data,xval)
      real(DP), dimension(-km:km), intent(IN) :: e_data
      !(in) å…¥åŠ›ãƒ•ãƒ¼ãƒªã‚¨ãƒ‡ãƒ¼ã‚¿

      real(DP), intent(in)                    :: xval
      ! è£œé–“ã™ã‚‹ç‚¹ã®åº§æ¨™

      real(DP)                                :: Interpolate_e
      ! è£œé–“ã—ãŸçµæžœã®å€¤

      integer :: k
      ! DO æ–‡å¤‰æ•°

      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) )
      enddo

    end function Interpolate_e

    function a_Interpolate_ae(ae_data,xval)
      real(DP), dimension(:,-km:), intent(IN)   :: ae_data
      !(in) å…¥åŠ›ãƒ•ãƒ¼ãƒªã‚¨ãƒ‡ãƒ¼ã‚¿

      real(DP), intent(in)                      :: xval
      ! è£œé–“ã™ã‚‹ç‚¹ã®åº§æ¨™

      real(DP), dimension(size(ae_data,1))      :: a_Interpolate_ae
      ! è£œé–“ã—ãŸçµæžœã®å€¤

      integer :: k
      ! DO æ–‡å¤‰æ•°

      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) )
      enddo

    end function a_Interpolate_ae

end module ae_module_fftj