mfcoef.f90

数值计算和数值分析在Fortran下的特殊函数库,是数值计算的必备

MODULE fcoef_func
 
! From the book "Computation of Special Functions"
!      by Shanjie Zhang and Jianming Jin
!   Copyright 1996 by John Wiley & Sons, Inc.
! The authors state:
!   "However, we give permission to the reader who purchases this book
!    to incorporate any of these programs into his or her programs
!    provided that the copyright is acknowledged."
 
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
 
CONTAINS


SUBROUTINE fcoef(kd, m, q, a, fc)

!       =====================================================
!       Purpose: Compute expansion coefficients for Mathieu
!                functions and modified Mathieu functions
!       Input :  m  --- Order of Mathieu functions
!                q  --- Parameter of Mathieu functions
!                KD --- Case code
!                       KD=1 for cem(x,q)  ( m = 0,2,4,...)
!                       KD=2 for cem(x,q)  ( m = 1,3,5,...)
!                       KD=3 for sem(x,q)  ( m = 1,3,5,...)
!                       KD=4 for sem(x,q)  ( m = 2,4,6,...)
!                A  --- Characteristic value of Mathieu
!                       functions for given m and q
!       Output:  FC(k) --- Expansion coefficients of Mathieu
!                       functions ( k= 1,2,...,KM )
!                       FC(1),FC(2),FC(3),... correspond to
!                       A0,A2,A4,... for KD=1 case, A1,A3,
!                       A5,... for KD=2 case, B1,B3,B5,...
!                       for KD=3 case and B2,B4,B6,... for
!                       KD=4 case
!       =====================================================

INTEGER, INTENT(IN)        :: kd
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN)      :: q
REAL (dp), INTENT(IN)      :: a
REAL (dp), INTENT(OUT)     :: fc(251)

REAL (dp)  :: f, f1, f2, f3, qm, s, s0, sp, ss, u, v
INTEGER    :: i, j, k, kb, km

fc = 0.0_dp
IF (q <= 1.0_dp) THEN
  qm = 7.5 + 56.1 * SQRT(q) - 134.7 * q + 90.7 * SQRT(q) * q
ELSE
  qm = 17.0 + 3.1 * SQRT(q) - .126 * q + .0037 * SQRT(q) * q
END IF
km = INT(qm+0.5*m)
IF (q == 0.0_dp) THEN
  IF (kd == 1) THEN
    fc((m+2)/2) = 1.0_dp
    IF (m == 0) fc(1) = 1.0_dp / SQRT(2.0_dp)
  ELSE IF (kd == 4) THEN
    fc(m/2) = 1.0_dp
  ELSE
    fc((m+1)/2) = 1.0_dp
  END IF
  RETURN
END IF
kb = 0
s = 0.0_dp
f = 1.0D-100
u = 0.0_dp
fc(km) = 0.0_dp
IF (kd == 1) THEN
  DO  k = km, 3, -1
    v = u
    u = f
    f = (a - 4*k*k) * u / q - v
    IF (ABS(f) < ABS(fc(k+1))) THEN
      kb = k
      fc(1) = 1.0D-100
      sp = 0.0_dp
      f3 = fc(k+1)
      fc(2) = a / q * fc(1)
      fc(3) = (a-4.0_dp) * fc(2) / q - 2.0_dp * fc(1)
      u = fc(2)
      f1 = fc(3)
      DO  i = 3, kb
        v = u
        u = f1
        f1 = (a - 4*(i-1)**2) * u / q - v
        fc(i+1) = f1
        IF (i == kb) f2 = f1
        IF (i /= kb) sp = sp + f1 * f1
      END DO
      sp = sp + 2.0_dp * fc(1) ** 2 + fc(2) ** 2 + fc(3) ** 2
      ss = s + sp * (f3/f2) ** 2
      s0 = SQRT(1.0_dp/ss)
      DO  j = 1, km
        IF (j <= kb+1) THEN
          fc(j) = s0 * fc(j) * f3 / f2
        ELSE
          fc(j) = s0 * fc(j)
        END IF
      END DO
      GO TO 160
    ELSE
      fc(k) = f
      s = s + f * f
    END IF
  END DO
  fc(2) = q * fc(3) / (a - 4.0_dp - 2.0_dp*q*q/a)
  fc(1) = q / a * fc(2)
  s = s + 2.0_dp * fc(1) ** 2 + fc(2) ** 2
  s0 = SQRT(1.0_dp/s)
  DO  k = 1, km
    fc(k) = s0 * fc(k)
  END DO
ELSE IF (kd == 2 .OR. kd == 3) THEN
  DO  k = km, 3, -1
    v = u
    u = f
    f = (a - (2*k-1)**2) * u / q - v
    IF (ABS(f) >= ABS(fc(k))) THEN
      fc(k-1) = f
      s = s + f * f
    ELSE
      kb = k
      f3 = fc(k)
      GO TO 80
    END IF
  END DO
  fc(1) = q / (a - 1.0_dp- (-1)**kd*q) * fc(2)
  s = s + fc(1) * fc(1)
  s0 = SQRT(1.0_dp/s)
  DO  k = 1, km
    fc(k) = s0 * fc(k)
  END DO
  GO TO 160

  80 fc(1) = 1.0D-100
  fc(2) = (a - 1.0_dp - (-1)**kd*q) / q * fc(1)
  sp = 0.0_dp
  u = fc(1)
  f1 = fc(2)
  DO  i = 2, kb - 1
    v = u
    u = f1
    f1 = (a - (2*i-1)**2) * u / q - v
    IF (i /= kb-1) THEN
      fc(i+1) = f1
      sp = sp + f1 * f1
    ELSE
      f2 = f1
    END IF
  END DO
  sp = sp + fc(1) ** 2 + fc(2) ** 2
  ss = s + sp * (f3/f2) ** 2
  s0 = 1.0_dp / SQRT(ss)
  DO  j = 1, km
    IF (j < kb) fc(j) = s0 * fc(j) * f3 / f2
    IF (j >= kb) fc(j) = s0 * fc(j)
  END DO
ELSE IF (kd == 4) THEN
  DO  k = km, 3, -1
    v = u
    u = f
    f = (a - 4*k*k) * u / q - v
    IF (ABS(f) >= ABS(fc(k))) THEN
      fc(k-1) = f
      s = s + f * f
    ELSE
      kb = k
      f3 = fc(k)
      GO TO 130
    END IF
  END DO
  fc(1) = q / (a - 4.0_dp) * fc(2)
  s = s + fc(1) * fc(1)
  s0 = SQRT(1.0_dp/s)
  DO  k = 1, km
    fc(k) = s0 * fc(k)
  END DO
  GO TO 160
  130   fc(1) = 1.0D-100
  fc(2) = (a-4.0_dp) / q * fc(1)
  sp = 0.0_dp
  u = fc(1)
  f1 = fc(2)
  DO  i = 2, kb - 1
    v = u
    u = f1
    f1 = (a - 4*i*i) * u / q - v
    IF (i /= kb-1) THEN
      fc(i+1) = f1
      sp = sp + f1 * f1
    ELSE
      f2 = f1
    END IF
  END DO
  sp = sp + fc(1) ** 2 + fc(2) ** 2
  ss = s + sp * (f3/f2) ** 2
  s0 = 1.0_dp / SQRT(ss)
  DO  j = 1, km
    IF (j < kb) fc(j) = s0 * fc(j) * f3 / f2
    IF (j >= kb) fc(j) = s0 * fc(j)
  END DO
END IF
160 IF (fc(1) < 0.0_dp) THEN
  DO  j = 1, km
    fc(j) = -fc(j)
  END DO
END IF
RETURN
END SUBROUTINE fcoef



SUBROUTINE cva2(kd, m, q, a)

!       ======================================================
!       Purpose: Calculate a specific characteristic value of
!                Mathieu functions
!       Input :  m  --- Order of Mathieu functions
!                q  --- Parameter of Mathieu functions
!                KD --- Case code
!                       KD=1 for cem(x,q)  ( m = 0,2,4,...)
!                       KD=2 for cem(x,q)  ( m = 1,3,5,...)
!                       KD=3 for sem(x,q)  ( m = 1,3,5,...)
!                       KD=4 for sem(x,q)  ( m = 2,4,6,...)
!       Output:  A  --- Characteristic value
!       Routines called:
!             (1) REFINE for finding accurate characteristic
!                 value using an iteration method
!             (2) CV0 for finding initial characteristic
!                 values using polynomial approximation
!             (3) CVQM for computing initial characteristic
!                 values for q