!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  C2PACK/ISPACK-F90
!
!履歴  2001/10/07  竹広真一
!      2001/12/26  竹広真一  関数, 変数命名法変更
!
module c2pack
  implicit none

  integer            :: im=32, jm=8     ! 格子点の設定(X,Y)
  integer            :: km=10, lm=5     ! 切断波数の設定(X,Y)
  double precision   :: xl=2.0, yl=1.0   ! 領域の大きさ

  integer,dimension(5)                       :: itj
  double precision,dimension(:),allocatable  :: tj
  integer,dimension(5)                       :: iti
  double precision,dimension(:),allocatable  :: ti

  double precision, dimension(:),   allocatable  :: wg, ws, wgj
  double precision, dimension(:,:), allocatable  :: gg_work,es_work,ec_work
  double precision, parameter  ::  pi=3.1415926535897932385D0

  private
  public c2initial
  public gg_es, gg_ec, es_gg, ec_gg
  public es_lapla_es, es_laplainv_es, es_dx_es, ec_dy_es, es_jacobian_es_es
  public ec_lapla_ec, ec_dx_ec, es_dy_ec, ec_jacobian_es_ec

  save im, jm, km, lm, itj, tj, iti, ti, xl, yl

  contains
  !--------------- 初期化 -----------------
    subroutine c2initial(i,j,k,l,xlength,ylength)

      integer,intent(in) :: i, j           ! 格子点の設定(X,Y)
      integer,intent(in) :: k, l           ! 切断波数の設定(X,Y)

      double precision,intent(in)  :: xlength, ylength   ! 領域の大きさ

      im = i       ; jm = j
      km = k       ; lm = l
      xl = xlength ; yl = ylength

      allocate(tj(jm*6),ti(im*2))
      allocate(wg((jm+1)*im))
      allocate(ws((2*km+1)*(lm+1)),wgj((jm+1)*im*3))
      allocate(gg_work(0:jm,0:im-1))
      allocate(es_work(-km:km,lm),ec_work(-km:km,0:lm))

      call c2init(jm,im,itj,tj,iti,ti)
    end subroutine c2initial

  !--------------- 基本変換 -----------------
    function gg_es(es) ! スペクトル SIN(Y) -> 台形格子
      double precision, dimension(0:jm,0:im-1)              :: gg_es
      double precision, dimension(-km:km,lm), intent(in)    :: es

      call c2s2ga(lm,km,jm,im,es,gg_es,wg,itj,tj,iti,ti,1)
    end function gg_es

    function gg_ec(ec)  ! スペクトル COS(Y) -> 台形格子
      double precision, dimension(0:jm,0:im-1)              :: gg_ec
      double precision, dimension(-km:km,0:lm), intent(in)  :: ec

      call c2s2ga(lm,km,jm,im,ec,gg_ec,wg,itj,tj,iti,ti,2)
    end function gg_ec

    function es_gg(gg)  ! 台形格子 -> スペクトル SIN(Y)
      double precision, dimension(-km:km,lm)                :: es_gg
      double precision, dimension(0:jm,0:im-1), intent(in)  :: gg

      gg_work = gg
      call c2g2sa(lm,km,jm,im,gg_work,es_gg,wg,itj,tj,iti,ti,1)
    end function es_gg

    function ec_gg(gg)  ! 台形格子 -> スペクトル COS(Y)
      double precision, dimension(-km:km,0:lm)              :: ec_gg
      double precision, dimension(0:jm,0:im-1), intent(in)  :: gg

      gg_work = gg
      call c2g2sa(lm,km,jm,im,gg_work,ec_gg,wg,itj,tj,iti,ti,2)
    end function ec_gg

  !--------------- 微分計算 -----------------
    function es_lapla_es(es)   ! スペクトル SINY に作用する \lapla 演算子
      double precision, dimension(-km:km,lm)                :: es_lapla_es
      double precision, dimension(-km:km,lm), intent(in)    :: es
      integer k,l

      do l=1,lm
         do k=-km,km
            es_lapla_es(k,l) = -((2*pi*k/xl)**2+(pi*l/yl)**2)*es(k,l)
         enddo
      enddo
    end function es_lapla_es

    function ec_lapla_ec(ec)   ! スペクトル COSY に作用する \lapla 演算子
      double precision, dimension(-km:km,0:lm)                :: ec_lapla_ec
      double precision, dimension(-km:km,0:lm), intent(in)    :: ec
      integer k,l

      do l=0,lm
         do k=-km,km
            ec_lapla_ec(k,l) = -((2*pi*k/xl)**2+(pi*l/yl)**2)*ec(k,l)
         enddo
      enddo
    end function ec_lapla_ec

    function es_laplainv_es(es)   ! スペクトル SINY に作用する逆 \lapla 演算子
      double precision, dimension(-km:km,lm)                :: es_laplainv_es
      double precision, dimension(-km:km,lm), intent(in)    :: es
      integer k,l

      do l=1,lm
         do k=-km,km
            es_laplainv_es(k,l) = -es(k,l)/((2*pi*k/xl)**2+(pi*l/yl)**2)
         enddo
      enddo
    end function es_laplainv_es

    function es_dx_es(es)   ! スペクトル SINY に作用する x 微分演算子
      double precision, dimension(-km:km,lm)                :: es_dx_es
      double precision, dimension(-km:km,lm), intent(in)    :: es
      integer k,l

      do l=1,lm
         do k=-km,km
            es_dx_es(k,l)  =  (-2*pi*k/xl)*es(-k,l)
         enddo
      enddo
    end function es_dx_es

    function ec_dx_ec(ec)   ! スペクトル COS(Y) に作用する x 微分演算子
      double precision, dimension(-km:km,0:lm)                :: ec_dx_ec
      double precision, dimension(-km:km,0:lm), intent(in)    :: ec
      integer k,l

      do l=0,lm
         do k=-km,km
            ec_dx_ec(k,l)  =  (-2*pi*k/xl)*ec(-k,l)
         enddo
      enddo
    end function ec_dx_ec

    function ec_dy_es(es)   ! スペクトル SINY に作用する y 微分演算子
      double precision, dimension(-km:km,0:lm)              :: ec_dy_es
      double precision, dimension(-km:km,lm), intent(in)    :: es
      integer k,l

      do k=-km,km
         ec_dy_es(k,0)  =  0.0
      enddo
      do l=1,lm
         do k=-km,km
            ec_dy_es(k,l)  =  (pi*l/yl)*es(k,l)
         enddo
      enddo
    end function ec_dy_es

    function es_dy_ec(ec)   ! スペクトル COSY に作用する y 微分演算子
      double precision, dimension(-km:km,lm)                 :: es_dy_ec
      double precision, dimension(-km:km,0:lm), intent(in)   :: ec
      integer k,l

      do l=1,lm
         do k=-km,km
            es_dy_ec(k,l)  =  -(pi*l/yl)*ec(k,l)
         enddo
      enddo
    end function es_dy_ec

    function es_jacobian_es_es(es_a,es_b) !スペクトル SINY に作用するヤコビアン
      double precision, dimension(-km:km,lm)                :: es_jacobian_es_es
      double precision, dimension(-km:km,lm), intent(in)    :: es_a,es_b

      integer k,l

      call c2ajcb(lm,km,jm,im,es_a,es_b,es_work,ws,wgj,itj,tj,iti,ti)

      do l=1,lm
         do k=-km,km
            es_jacobian_es_es(k,l) = (2*pi/xl)*(pi/yl)*es_work(k,l)
         enddo
      enddo
    end function es_jacobian_es_es

    function ec_jacobian_es_ec(es,ec)  ! スペクトル COS(Y) に作用するヤコビアン
      double precision, dimension(-km:km,0:lm)              :: ec_jacobian_es_ec
      double precision, dimension(-km:km,lm), intent(in)    :: es
      double precision, dimension(-km:km,0:lm), intent(in)  :: ec
      integer k,l

      call c2ajcc(lm,km,jm,im,es,ec,ec_work,ws,wgj,itj,tj,iti,ti)

      do l=0,lm
         do k=-km,km
            ec_jacobian_es_ec(k,l) = (2*pi/xl)*(pi/yl)*ec_work(k,l)
         enddo
      enddo
    end function ec_jacobian_es_ec

  end module c2pack
