Liren Large Software Subsidy 10/10.iso/l/l292/1.ddi/SPECFUNC.FOR

Text File | 1989-10-12 | 20.4 KB | 719 lines

      FUNCTION ArcTanh(x)
      REAL x
        ArcTanh = 0.5 * LOG((1.0 + x) / (1.0 - x))
      END !FUNCTION

      FUNCTION EvalAssocConf(a, b, z)
      INCLUDE 'STDHDR.FOR'
      REAL result, a, b, z, temp1, temp2,Calc1, Calc2, Calc3
      REAL PiCalc, Kum1, Kum2, Exp1, Subtotal, Gam1, Gam2
      INTEGER ErrorFlag
          Temp1 = 0
          Temp2 = 0
          PiCalc = pi / SIN(pi * b)
          Kum1 = KumrConf(a, b, z)
          Gam1 = Gamma(1.0 + a - b, ErrorFlag) * Gamma(b, ErrorFlag)
          Temp1 = Gamma(1.0 + a - b, ErrorFlag)
          Temp2 = Gamma(b, ErrorFlag)
          Gam1 = Gamma(1.0 + a - b, ErrorFlag) * Gamma(b, ErrorFlag)
          Exp1 = CalcPower(z, 1.0 - b)
          Kum2 = KumrConf(1.0 + a - b, 2.0 - b, z)
          Gam2 = Gamma(a, ErrorFlag) * Gamma(2.0 - b, ErrorFlag)
          Temp1 = Gamma(a, ErrorFlag)
          Temp2 = Gamma(2.0 - b, ErrorFlag)
          Calc1 = Kum1 / Gam1
          Calc2 = Exp1 * Kum2
          Calc3 = Calc2 / Gam2
          Subtotal = Calc1 - Calc3
          result = PiCalc * Subtotal
          EvalAssocConf = result
      END !FUNCTION


      FUNCTION AssocConf(a, b, z)
      INCLUDE 'STDHDR.FOR'
      REAL a, b, z, bu,bl
      result = 0.0
      IF ( ABS(z) .LT. 1E-10) THEN
        result = 0.0
      ELSE
        IF ( ABS(b - INT(b)) .LT. 0.0001 .OR.
     +       ABS(b - NINT(b)) .LT. 0.0001) THEN
            bu = b + 0.0005
            bl = b - 0.0005
          result  = EvalAssocConf(a, bl , z) +
     +       EvalAssocConf(a, bu, z) / 2.0
        ELSE
          result  = EvalAssocConf(a, b , z)

        END IF
      END IF
      AssocConf = result
      END !FUNCTION

      FUNCTION Bessel(order, x, ErrorFlag)
      INCLUDE 'STDHDR.FOR'
      REAL result,order, x, d, f, j, sum, NearZero, Tolerance
      INTEGER c, ErrorFlag
      PARAMETER (nearzero = 1E-20)
      PARAMETER (Tolerance = 1E-10)
      result = 0.0
      IF ( order .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        IF ( ((x .LT. 0) .AND. (ABS(x - NINT(x)) .LT. nearzero))) THEN
          ErrorFlag = 2
        ELSE
          IF ( ABS(x) .LT. nearzero) THEN
             IF ( ABS(order) .LT. nearzero) THEN
               result = 1
             ELSE
               result = 0
             END IF
          ELSE
            IF ( x .LE. 15.0) THEN
              f = 1.0 / Gamma(REAL(order) + 1.0, ErrorFlag)
              j = -((x ** 2) / 4.0)
              d = order + 1
              sum = f
              c = 1
              DO WHILE (ABS(f) .GE. Tolerance  )
                f = f * j / (REAL(c) * d)
                d = d + 1
                c = c + 1
                sum = sum + f
              END DO
              IF ( x .GT. 0) THEN
                result = sum * EXP(order * LOG(x / 2.0))
              ELSE
                result = (1.0 - 2.0 * (MOD(INT(order), 2))) *
     +                  sum * EXP(order * LOG(ABS(x / 2.0)))
              END IF
            ELSE
              result = SQRT(2.0 / (pi * x)) *
     +                   COS(x - pi / 4.0 - order * pi * 2.0)
            END IF
         END IF
        END IF
      END IF
      Bessel = result
      END !FUNCTION  Bessel

      FUNCTION Beta(z, w)
      REAL z, w, b
        b = LogGamma(z) + LogGamma(w) - LogGamma(z + w)
        Beta = EXP(b)
      END !FUNCTION

     !!!------------------------------------------------------------------
     !!! Calculated a polynomial
     !!!------------------------------------------------------------------
      FUNCTION CalcPoly(XIn, CoefVector, order)
      REAL CoefVector(*), xin
      INTEGER order, i

      YOut = CoefVector(order)
      DO i = order - 1, 0,  -1
        YOut = YOut * XIn + CoefVector(i)
      END DO
      CalcPoly = YOut
      END !FUNCTION

      FUNCTION CalcPower(x, a)
      REAL x, a
      CalcPower = EXP(a * LOG(x))
      END !FUNCTION


      FUNCTION  Combin (n, m)
      !!! Combinatorial bracket n-over-m, m an integer
      REAL n, nt, result
      INTEGER i,m,im
        result = 1
        nt = n
        DO I=1,M
          im = M-I+1
          result =  result * (nt/ REAL(im))
          nt = nt - 1.0
        END DO
        Combin = result
      END !FUNCTION   Combin

      FUNCTION Coshy (z)
      REAL z
       CoshY = (EXP(z) + EXP(-z)) / 2.0
      END !FUNCTION  Cosh

      SUBROUTINE DisplayErrorMessage(ErrorFlag)
      INTEGER ErrorFlag
      SELECT CASE (ErrorFlag)

       CASE (0)
          WRITE (*,*) 'No Errors Found.'
       CASE (1)
          WRITE (*,*) 'ERROR!!  Argument must be a nonzero number'
       CASE (2)
          WRITE (*,*) 'ERROR!!  Argument must not be a negative'
       CASE (3)
          WRITE (*,*)
     +     'ERROR!!  Order of polynomial must be greater than zero'
       CASE (4)
         WRITE (*,*) 'ERROR!!  Argument must be greater than -1'
       END SELECT
       END !SUB

       FUNCTION ErrFunc(x, ErrorFlag)
       INCLUDE 'STDHDR.FOR'
       REAL x, result, p, temp
       INTEGER ErrorFlag
       REAL PolyCoef(0 : 15), NearZero
       PARAMETER (nearzero = 1E-20)
       result = 0.0
       IF ( ABS(x) .LT. nearzero) THEN
         ErrorFlag = 1
       ELSE
         IF ( x .GT. 4.0) THEN
           result = 1.0
         ELSE
           IF ( x .GT. 1.5) THEN
             result = ErrFuncIter(x)
           ELSE
             p = x * x
             PolyCoef(0) = 0.0
             PolyCoef(0) = 0.0
             PolyCoef(1) = 0.6666666999999999
             PolyCoef(2) = 0.2666667
             PolyCoef(3) = 0.07619048
             PolyCoef(4) = 0.01693122
             PolyCoef(5) = 3.078403E-03
             PolyCoef(6) = 4.736005E-04
             PolyCoef(7) = 6.314673E-05
             PolyCoef(8) = 7.429027E-06
             PolyCoef(9) = 7.820028E-07
             PolyCoef(10) = 7.447646E-08
             PolyCoef(11) = 6.476214E-09
             temp = x + x * CalcPoly(p, PolyCoef, 11)
             result = 2.0 * temp * EXP(-p) / SQRT(pi)
           END IF
         END IF
       END IF
       ErrFunc = result
      END !FUNCTION ErrFunc


      FUNCTION ErrFuncComp(x, ErrorFlag)
      REAL x
      INTEGER ErrorFlag
      ErrFuncComp = 1 - ErrFunc(x, ErrorFlag)
      END !FUNCTION



      FUNCTION ErrFuncI(x, y)
      !!!   The imaginary part of complex error function, z = x + iy
      INCLUDE 'STDHDR.FOR'
      REAL x, y, jy,x2, xy2,sum, part, g
      INTEGER j
      x2 = 2.0 * x
      xy2 = x2 * y
      sum = 0.0
      IF ( ABS(x) .LT. 0.0000001) THEN
        part = y / pi
      ELSE
        part = EXP(-(x ** 2)) / (pi * x2) * SIN(xy2)
      END IF
      sum = sum + part
      part = 0
      DO j = 1, 20
        jy = j * y
        g = x2 * Cosh(jy) * SIN(xy2) + j * Sinh(jy) * COS(xy2)
        part = part * EXP(-((j ** 2) / 4.0)) / ((j**2) + (x2**2)) * g
      END DO
      sum = sum + part * EXP(-(x ** 2)) * 2 / pi
      ErrFuncI = sum
      END !FUNCTION   ErrFunci



      FUNCTION ErrFuncIter(x)
      INCLUDE 'STDHDR.FOR'
      REAL x, a, s1, t, result, Tolerance
      INTEGER i
      PARAMETER (Tolerance = 1E-10)
      result = 0.0
      a = x ** 2
      s1 = x
      t = x
      i = 1
      DO WHILE (t .GT. (Tolerance * s1))
        s2 = s1
        t = 2.0 * t * (a / (1.0 + 2.0 * i))
        s1 = t + s2
        i = i + 1
      END DO
      result = 2.0 * s1 * (EXP(-a)) / SQRT(pi)
      ErrFuncIter = result
      END !FUNCTION ErrFuncIter



      FUNCTION ErrFuncR(x, y)
      INCLUDE 'STDHDR.FOR'
      !!!    The real part of complex error function, z = x + iy
      REAL x, y, x2, xy2, sum, part, jy, f
      INTEGER ErrorFlag
      sum = 0.0
      x2 = 2.0 * x
      xy2 = x2 * y
      IF ( ABS(x) .GT. 0.0) THEN
        part = EXP(-(x ** 2))/x * (1.0 / (2.0*pi)) * (1.0- COS(xy2))
      ELSE
        part = 0.0
      END IF
      sum = ErrFunc(x, ErrorFlag)
      sum = sum + part
      part = 0.0
      DO j = 1, 10
        jy = j * y
        f = x2 - x2 * COS(xy2) * Cosh(jy) + j * Sinh(jy) * SIN(xy2)
        part = part + EXP(-(j ** 2/4.0))/(j ** 2 + 4 * x ** 2) * f
      END DO
      sum = sum + part * 2.0 / pi * EXP(-(x ** 2))
      ErrFuncR = sum
      END !FUNCTION  ErrFuncr

      FUNCTION EvalGamma(x)
      INCLUDE 'STDHDR.FOR'
      REAL result, x, y, f,NearZero
      PARAMETER (nearzero = 1E-20)

        result = 0.0
        y = x
        f = 1
        DO WHILE (y .GT. 1)
          y = y - 1
          f = f * y
        END DO
        IF ( y .EQ. 1) THEN
          result = f
        ELSE
          IF ( y .GT. 0.5) THEN
            result = (f * pi) / (SIN(pi * y) * StandGamma(1.0 - y))
          ELSE
            result = f * StandGamma(y)
          END IF
        END IF
        EvalGamma = result
      END !FUNCTION   EvalGamma

      FUNCTION Gamma(x, ErrorFlag)
      INCLUDE 'STDHDR.FOR'
      REAL result, x, NearZero
      INTEGER ErrorFlag
      PARAMETER (nearzero = 1E-20)
      result = 0.0
      IF (ABS(x) .LT. nearzero) THEN
        ErrorFlag = 1
      ELSE
        IF ((x .LT. 0) .AND. (ABS(x - NINT(x)) .LT. nearzero)) THEN
          ErrorFlag = 2
        ELSE
          IF ( x .EQ. 1) THEN
            result = 1
          ELSE
            IF ( x .LT. 0.0) THEN
              result = -(pi / (x * SIN(pi * x) * EvalGamma(-x)))
            ELSE
              result = EvalGamma(x)
            END IF
          END IF
        END IF
      END IF
      Gamma = result
      END !FUNCTION   Gamma

      FUNCTION GaussHyper(a, b, c, z)
      !!!    Gauss hypergeometric function
      REAL a,b,c,z,g,sum, Tolerance
      INTEGER j
      PARAMETER (Tolerance = 1E-10)
      g = 1
      j = 1
      sum = 0 .0
      DO WHILE (j .LE. 80 .AND. ABS(G) .GE. Tolerance)
        sum = sum + g
        g = g * a * b * z / (REAL(j) * c)
        j = j + 1
        a = a + 1
        b = b + 1
        c = c + 1
      END DO
      GaussHyper = sum
      END !FUNCTION   GaussHyper

      FUNCTION Hermite(n, x, ErrorFlag)
      !!! Hermite Polynomial of Order n
      REAL result, X
      IF ( n .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        SELECT CASE (n)
          CASE (0)
            result = 1.0
          CASE (1)
            result = 2.0 * x
          CASE (2)
            result = 4.0 * (x ** 2) - 2.0
          CASE (3)
            result = 8.0 * CalcPower(x, 3.0) - (12.0 * x)
          CASE (4)
            result = 16.0 * CalcPower(x, 4.0) -
     +           48.0 * (x ** 2) + 12.0
          CASE (5)
            result = 32.0 * CalcPower(x, 5.0) - 160.0 *
     +           CalcPower(x, 3.0) + 120.0 * x
        CASE DEFAULT
          result = x * EXP(n * LOG(2)) *
     +       AssocConf(0.5 - REAL(n) / 2.0, 3.0 / 2.0, (x ** 2))
        END SELECT
      END IF
      Hermite = result
      END !FUNCTION  Hermite

      !!!---------------------------------------
      !!!    I N C O M P L E T E  B E T A
      !!!---------------------------------------
      FUNCTION IBeta(a, b, x)
      REAL a,b,x,am,bm,az,qab,qap,qam,bz,m,em,tem,d,ap
      REAL bpp, aold, Tolerance
      PARAMETER (Tolerance = 1E-10)
      am = 1.0
      bm = 1.0
      az = 1.0
      aold = 0.0
      qab = a + b
      qap = a + 1.0
      qam = a - 1.0
      bz = 1.0 - qab * x / qap
      m = 1.0
      DO WHILE (ABS(az - aold) .GE. (Tolerance * ABS(az)))
        em = m
        tem = 2.0 * em
        d = em * (b - m) * x / ((qam + tem) * (a + tem))
        ap = az + d * am
        bp = bz + d * bm
        d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem))
        app = ap + d * az
        bpp = bp + d * bz
        aold = az
        am = ap / bpp
        bm = bp / bpp
        az = app / bpp
        bz = 1.0
        m = m + 1.0
      END DO
      IBeta = az
      END !FUNCTIONIBeta



      FUNCTION IncompleteBeta(a, b, x)
      REAL a,b,x,bt,result, NearZero
      PARAMETER (nearzero = 1E-20)
      result = 0.0
      IF (x .LT. 0.0 .OR. x .GT. 1.0) THEN
        WRITE (*,*) 'Improper Value passed IncompleteBeta'
      ELSE
        IF ((ABS(x) .LT. nearzero) .OR. ((1.0-x) .LT. nearzero)) THEN
          bt = 0.0
        ELSE
          bt = LogGamma(a + b) + a*LOG(x) +
     +         b * LOG(1.0 - x) - LogGamma(a) - LogGamma(b)
        END IF
        bt = EXP(bt)
        IF ( x .LT. (a + 1.0) / (a + b + 2.0)) THEN
          result = bt * (IBeta(a, b, x) / a)
        ELSE
          result = 1 - bt * IBeta(b, a, 1.0 - x) / b
        END IF
        IncompleteBeta = result
      END IF
      END !FUNCTION



      FUNCTION IncompleteGamma(a, x)
      REAL a, x, result, glog, s1, dif, ir, Tolerance
      PARAMETER (Tolerance = 1E-10)
      result = 0.0
      IF (a .LT. 0.0 .OR. x .LE. 0.0) THEN
        result = 0.0
      ELSE
        gLOG = LogGamma(a)
        ir = a
        s1 = 1.0 / a
        dif = s1
        DO WHILE (ABS(dif) .GT. ABS(s1) * Tolerance)
          ir = ir + 1.0
          dif = dif * x / ir
          s1 = s1 + dif
        END DO
        result = s1 * EXP(-x + a * LOG(x) - gLOG)
      END IF
      IncompleteGamma = result
      END !FUNCTION

      FUNCTION IncompleteGammaComp(a, x)
      REAL a,x
      IncompleteGammaComp = 1.0 - IncompleteGamma(a, x)
      END !FUNCTION

      FUNCTION Jacobi(n, a, b, x, ErrorFlag)
      REAL result,a,b,x
      INTEGER n, ErrorFlag
      !!! Jacobi polynomial of order n, superscripts (a,b)
      result = 0.0
      IF ( n .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        IF ((a .LE. -1.0) .OR. (b .LE. -1.0)) THEN
          ErrorFlag = 4
        ELSE
         result = Combin(REAL(n) + a, n) *
     +        GaussHyper(-REAL(n), n+a+b+1.0, a+1.0, (1.0-x)/2.0)
        END IF
      END IF
      Jacobi = result
      END !FUNCTION  Jacobi

      FUNCTION KumrConf(a, b, z)
      REAL a,b,z, sum, g, Tolerance
      INTEGER j
      PARAMETER (Tolerance = 1E-10)
      !!!  Kummer (confluent) hypergeometric function
      sum = 0.0
      g = 1.0
      j = 1
      DO WHILE ((j .LT. 80) .AND. (ABS(g) .LT. Tolerance))
        sum = sum + g
        g = g * a * z / (b * j)
        j = j + 1
        a = a + 1.0
        b = b + 1.0
      END DO

      KumrConf = sum
      END !FUNCTION    KumrConf

      FUNCTION Laguerre(n, a, x, ErrorFlag)
      !!! Laguerre function or order n, superscript a
      REAL result,a, x,num
      INTEGER n, ErrorFlag
      IF ( n .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        IF ( a .LE. -1) THEN
          ErrorFlag = 4
        ELSE
          IF ( (n/2) .EQ. 0) THEN
            num = 1
          ELSE
            num = -1
          END IF
          SELECT CASE (n)
            CASE (0)
              result = 1
            CASE (1)
              result = -x + 1
            CASE (2)
              result = 0.5 * ((x ** 2) - 4 * x + 2)
            CASE (3)
              result = 1.0 / 6.0 *
     +         (-CalcPower(x, 3.0) + 9.0 * (x ** 2) - 18.0 * x + 6)
            CASE (4)
              result = 1.0 / 24.0 * (CalcPower(x, 4.0) - 16 *
     +          CalcPower(x, 3.0) + 72.0 * (x ** 2) - 96.0*x + 24.0)
          CASE DEFAULT
              result = num * AssocConf(-REAL(n), a + 1, x) /
     +           Gamma(REAL(n) + 1.0, ErrorFlag)
          END SELECT
        END IF
      END IF
      Laguerre = result
      END !FUNCTION  Laguerre

      FUNCTION Legend(n, x, ErrorFlag)
      REAL x
      INTEGER result,n, ErrorFlag
      !!! Legendre Polynomial of Order n

      IF ( n .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        SELECT CASE (n)
          CASE (0)
             result = 1.0
          CASE (1)
             result = x
          CASE (2)
            result = 0.5 * (3.0 * x * x - 1.0)
          CASE (3)
             result = 0.5 * (5.0 * CalcPower(x, 3.0) - (3.0 * x))
         CASE (4)
             result = 1.0 / 8.0 * (35.0 * CalcPower(x, 4.0) -
     +                 30.0 * (4.0) + 3.0)
         CASE (5)
             result = 1.0 / 8.0 * (63.0 * CalcPower(x, 5.0) -
     +                 70.0 * CalcPower(x, 3.0) + 15.0 * x)
         CASE (6)
             result = 1.0 / 16.0 * (231.0 * CalcPower(x, 6.0) -
     +                 315.0 * CalcPower(x, 4.0) + 105.0 * 4.0 - 5.0)
         CASE DEFAULT
            result = GaussHyper(-REAL(n), REAL(n)+1.0,1.0,(1.0-x)/2.0)
          END SELECT
        END IF
        Legend = result
      END !FUNCTION   Legend



      FUNCTION LogGamma(x)
      REAL x,stp,xx,temp,ser,coef(0 : 15)
      INTEGER i
      coef(0) = 76.18009173
      coef(1) = -86.50532033
      coef(2) = 24.01409822
      coef(3) = -1.231739516
      coef(4) = 0.00120858003
      coef(5) = -5.36382E-06
      stp = 2.50662827565
      xx = x - 1.0
      temp = xx + 5.5
      temp = (xx + 0.5) * LOG(temp) - temp
      ser = 1.0
      DO i = 0, 5
        xx = xx + 1.0
        ser = ser + coef(i) / xx
      END DO
      LogGamma = temp + LOG(stp * ser)
      END !FUNCTION

      FUNCTION ModBesselI(n, x, ErrorFlag)
      !!! Modified Bessel function i of order n
      REAL result, n, x, NearZero
      INTEGER ErrorFlag
      PARAMETER (nearzero = 1E-20)
       result = 0.0
       IF ( ABS(x) .LT. nearzero) THEN
         IF ( ABS(n) .LT. nearzero) THEN
             result = 1.0
         ELSE
             result = 0.0
         END IF
       ELSE
         result = EXP(-x) * EXP(n * LOG(x / 2.0)) *
     +     KumrConf(n+0.5, 2.0*n+1.0, 2*x)/Gamma(n+1, ErrorFlag)
      END IF
      ModBesselI = result
      END !FUNCTION  ModBesselI

      FUNCTION ModBesselK(n, x)
      INCLUDE 'STDHDR.FOR'
      !!! Modified bessel function k of order n
      REAL result,n, x, NearZero
      PARAMETER (nearzero = 1E-20)

      IF ( ABS(x) .LT. nearzero) THEN
         IF ( ABS(n) .LT. nearzero) THEN
             result = 1.0
         ELSE
            result = 0.0
         END IF
      ELSE
        result = EXP(n * LOG(2 * x)) * EXP(-x) * SQRT(pi) *
     +      AssocConf(n + 0.5, 2.0 * n + 1.0, 2.0 * x)
      END IF
      ModBesselK = result
      END !FUNCTION ModBesselK

      FUNCTION Sech(z)
      REAL z
      Sech = 1.0 / Cosh(z)
      END !FUNCTION  Sech

      FUNCTION SinhY(z)
      REAL z
      SinhY = (EXP(z) - EXP(-z)) / 2.0
      END !FUNCTION  Sinh

      FUNCTION StandGamma(y)
      REAL BigArry(0 : 20), y, f, temp

      BigArry(0) = 0.5772156649015329
      BigArry(1) = -.6558780715202538
      BigArry(2) = -0.0420026350340952
      BigArry(3) = 0.1665386113822915
      BigArry(4) = -0.0421977345555443
      BigArry(5) = -9.621971527877001e-03
      BigArry(6) = 0.007218943246663
      BigArry(7) = -0.0011651675918591
      BigArry(8) = -0.0002152416741149
      BigArry(9) = 0.0001280502823882
      BigArry(10) = -0.0000201348547807
      BigArry(11) = -0.0000012504934821
      BigArry(12) = 0.000001133027232
      BigArry(13) = -0.0000002056338417
      BigArry(14) = 0.000000006116095
      BigArry(15) = 0.0000000050020075
      BigArry(16) = -0.0000000011812746
      BigArry(17) = 0.0000000001043427
      BigArry(18) = 7.782299999999999D-12
      BigArry(19) = -0.0000000000036968
      BigArry(20) = 0.000000000000051
      f = y
      temp = f
      DO i = 0, 20
        f = f * y
        temp = temp + f * BigArry(i)
      END DO
      StandGamma = 1.0 / temp
      END !FUNCTION  StandGamma

      FUNCTION Tangent(x)
      REAL x
      Tangent = SIN(x) / COS(x)
      END !FUNCTION   Tan

      FUNCTION Tcheb(n, x, ErrorFlag)
      !!! Tchebyshev polynomial of order n
      INTEGER n, ErrorFlag
      REAL x, result
      IF ( n .LT. 0) THEN
        ErrorFlag = 3
      ELSE
        SELECT CASE (n)
          CASE (0)
               result = 1.0
          CASE (1)
               result = x
          CASE (2)
               result = 2.0 * (x ** 2) - 1.0
          CASE (3 )
               result = 4.0 * CalcPower(x, 3.0) - 3 * x
          CASE (4)
               result = 8.0 * CalcPower(x, 4.0) - 8.0 * (x ** 2) + 1.0
          CASE (5)
               result = 16.0 * CalcPower(x, 5.0) - 20.0 *
     +          CalcPower(x, 3.0) + 5.0 * x
          CASE (6)
               result = 32.0 * CalcPower(x, 6.0) - 48 *
     +          CalcPower(x, 4.0) + 18 * (x ** 2) - 1.0
         CASE DEFAULT
           result = GaussHyper(-REAL(n), REAL(n), 0.5, (1.0 - x) / 2.0)
         END SELECT
      END IF
      Tcheb = result
      END !FUNCTION   Tcheb