* PACKAGE UCHEB   !"  チェビシェフ(Chebyshev)多項式近似
*
*"  [HIS] 92/09/18(takepiro) Numerical Recipes より
*"        92/11/09(takepiro) ベクトル化チューニング
*"        93/01/14(takepiro) 和のとりかたを修正
*
**********************************************************************
      SUBROUTINE CHBC2G  !" チェビシェフ係数→グリッド(停留点)
     M              (  WDATA , 
     I                 CDATA , 
     C                 QXSI  , DRFAC ,
     F                 HGRAD , HFUNC , 
     D                 NMDIM , KMAX  , 
     W                 CHEB  , WORK       )
*
*   [PARAM] 
      INTEGER    NMDIM
      INTEGER    KMAX
*
*   [MODIFY]       
      REAL       WDATA ( NMDIM, KMAX )    !" チェビシェフグリッドデータ
*
*   [INPUT]
      REAL       CDATA ( NMDIM, KMAX )    !" チェビシェフ係数
*
*   [CONST]
      REAL       QXSI  ( KMAX )           !" チェビシェフポイント x_k
      REAL       DRFAC                    !" 微分のファクター
*
*   [FUNC]
      CHARACTER  HGRAD*4                  !" 微分フラグ
      CHARACTER  HFUNC*4                  !" 符号フラグ
*
*   [WORK] 
      REAL       CHEB  ( NMDIM, KMAX )    !" 作業用チェビシェフ係数
      REAL       WORK  ( NMDIM, KMAX )
*
*   [INTERNAL WORK] 
      INTEGER    NMAX
        PARAMETER ( NMAX = 10000 )
      REAL       DJ(NMAX), DJ1(NMAX), DJ2(NMAX)
      INTEGER    NM, J, K 
*
      IF ( NMAX .LT. NMDIM )THEN
         CALL MSGDMP( 'E','CHBC2G','WORK AREA IS NOT ENOUGH.' )
      ENDIF
*
*"   < 1. 微分の計算 >
*
      IF  ( HGRAD(1:1) .EQ. 'D' ) THEN
         CALL CHBDER
     O        ( CHEB  , 
     I          CDATA , DRFAC , 
     M          NMDIM , KMAX    )
      ELSE
         CALL COPY ( CHEB , CDATA , NMDIM*KMAX )
      ENDIF
*
*"   < 2. Σc_jnm Tj(x_k) の計算 >
*
*" Clenshaw's reccurence
*
      DO 1000 K = 1, KMAX
         CALL RESET( DJ1, NMAX )
         CALL RESET( DJ,  NMAX )
         DO 1100 J = KMAX-1 , 2, -1
            DO 1100 NM = 1, NMDIM
               DJ2(NM) = DJ1(NM)
               DJ1(NM) = DJ(NM)
               DJ(NM) = 2.0 * QXSI(K) * DJ1(NM) - DJ2(NM) + CHEB( NM,J )
 1100    CONTINUE
*
         DO 1000 NM = 1, NMDIM
            WORK( NM,K ) = QXSI(K) * DJ(NM) - DJ1(NM)
     &                    + 0.5 * CHEB( NM,1 )
     &                    + 0.5 * CHEB( NM,KMAX ) * (-1) ** (K-1)
 1000 CONTINUE
*
*"    < 3. オプション >
*
*"       < 3.1  足し込み >
*
      IF  ( HFUNC(1:1) .EQ. 'A' ) THEN
         DO 4100 K = 1, KMAX
            DO 4100 NM = 1, NMDIM
               WDATA( NM, K ) = WDATA( NM, K ) + WORK( NM, K )
 4100    CONTINUE
*
*"       < 3.2  差し引き >
*
      ELSE IF  ( HFUNC(1:1) .EQ. 'S' ) THEN
         DO 4200 K = 1, KMAX
            DO 4200 NM = 1, NMDIM
               WDATA( NM, K ) = WDATA( NM, K ) - WORK( NM, K )
 4200    CONTINUE
*
*"       < 3.3. マイナス >
*
      ELSE IF  ( HFUNC(1:1) .EQ. 'N' ) THEN
         DO 4300 K = 1, KMAX
            DO 4300 NM = 1, NMDIM
               WDATA( NM, K ) =  - WORK( NM, K )
 4300    CONTINUE
*
*"       < 3.4  通常 >
*
      ELSE 
         DO 4400 K = 1, KMAX
            DO 4400 NM = 1, NMDIM
               WDATA( NM, K ) =  WORK( NM, K )
 4400    CONTINUE
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE CHBG2C     !" グリッド(停留点)→チェビシェフ係数
     M              (  CDATA , 
     I                 WDATA , 
     C                 QCHEB , DRFAC ,
     F                 HGRAD , HFUNC , 
     D                 NMDIM , KMAX  , 
     W                 CHEB  , WORK       )
*
*   [PARAM] 
      INTEGER    NMDIM
      INTEGER    KMAX
*
*   [MODIFY]       
      REAL       CDATA ( NMDIM, KMAX )    !" チェビシェフグリッドデータ
*
*   [INPUT]
      REAL       WDATA ( NMDIM, KMAX )    !" チェビシェフ係数
*
*   [CONST]
      REAL       QCHEB ( KMAX , KMAX )    !" チェビシェフポイント x_k
      REAL       DRFAC                    !" 微分のファクター
*
*   [FUNC]
      CHARACTER  HGRAD*4                  !" 微分フラグ
      CHARACTER  HFUNC*4                  !" 符号フラグ
*
*   [WORK] 
      REAL       CHEB  ( NMDIM, KMAX )    !" 作業用チェビシェフ係数
      REAL       WORK  ( NMDIM, KMAX )
*
*   {INTERNAL WORK]
      INTEGER  NM, J, K, KDIM
*
      KDIM = KMAX - 1
*
*" チェビシェフ係数の計算
*
      DO 2000 J = 1, KMAX
         DO 2000 NM = 1, NMDIM
            WORK( NM,J ) = 0.5 * WDATA( NM,1 )
     &                   + 0.5 * WDATA( NM,KMAX ) * QCHEB( J,KMAX )
 2000 CONTINUE 
*
      DO 2100 K = 2, KMAX-1
         DO 2100 J = 1, KMAX
            DO 2100 NM = 1, NMDIM
               WORK( NM,J ) = WORK( NM,J ) 
     &                      + WDATA( NM,K ) * QCHEB( J,K )
 2100 CONTINUE
*
      DO 2200 J = 1, KMAX
         DO 2200 NM = 1, NMDIM
            CHEB( NM,J ) = ( 2.0 / KDIM ) * WORK( NM,J )
 2200 CONTINUE
*
      IF  ( HGRAD(1:1) .EQ. 'D' ) THEN
         CALL CHBDER
     O          ( CDATA, 
     I            CHEB  , DRFAC, 
     M            NMDIM , KMAX   )
      ELSE
         CALL COPY ( CDATA , CHEB , NMDIM*KMAX )
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE CHBDER    !" 微分のチェビシェフ係数を計算
     O              ( CDDATA, 
     I                CDATA , DRFAC , 
     M                NMDIM , KMAX    )
*
*   [PARAM] 
      INTEGER    NMDIM
      INTEGER    KMAX
*
*   [OUTPUT]       
      REAL       CDDATA( NMDIM, KMAX )    !" 微分のチェビシェフ係数
*
*   [INPUT]
      REAL       CDATA ( NMDIM, KMAX )    !" チェビシェフ係数
      REAL       DRFAC                    !" 微分のファクター
*
      INTEGER    NM, J
*
      DO 1000 NM = 1, NMDIM
         CDDATA( NM, KMAX )   = 0.        !" 最後の項の和のとりかたに注意
         CDDATA( NM, KMAX-1 ) = 2 *( KMAX-1 ) * CDATA( NM, KMAX ) / 2
 1000 CONTINUE 
*
      IF( KMAX .GE. 3 )THEN
         DO 1100 J = KMAX-2, 1, -1
            DO 1100 NM = 1, NMDIM
               CDDATA( NM,J ) = 
     &             CDDATA( NM,J+2 ) + 2 * J * CDATA( NM,J+1 )
 1100    CONTINUE
      ENDIF
*
      DO 2000 J = 1, KMAX
         DO 2000 NM = 1, NMDIM
            CDDATA( NM,J ) = DRFAC * CDDATA ( NM,J ) 
 2000 CONTINUE 
*
      RETURN
      END
