Directory of image this file is from
This file as a plain text file

C
C     ..................................................................
C
C        SUBROUTINE CEL2
C
C        PURPOSE
C           COMPUTES THE GENERALIZED COMPLETE ELLIPTIC INTEGRAL OF
C           SECOND KIND.
C
C        USAGE
C           CALL CEL2(RES,AK,A,B,IER)
C
C        DESCRIPTION OF PARAMETERS
C           RES   - RESULT VALUE
C           AK    - MODULUS (INPUT)
C           A     - CONSTANT TERM IN NUMERATOR
C           B     - FACTOR OF QUADRATIC TERM IN NUMERATOR
C           IER   - RESULTANT ERROR CODE WHERE
C                   IER=0  NO ERROR
C                   IER=1  AK NOT IN RANGE -1 TO +1
C
C        REMARKS
C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1.7E38 IF B IS                0
C           POSITIVE, TO -1.7E38 IF B IS NEGATIVE.                              0
C           SPECIAL CASES ARE
C           K(K) OBTAINED WITH A = 1, B = 1
C           E(K) OBTAINED WITH A = 1, B = CK*CK WHERE CK IS
C           COMPLEMENTARY MODULUS.
C           B(K) OBTAINED WITH A = 1, B = 0
C           D(K) OBTAINED WITH A = 0, B = 1
C           WHERE K, E, B, D DEFINE SPECIAL CASES OF THE GENERALIZED
C           COMPLETE ELLIPTIC INTEGRAL OF SECOND KIND IN THE USUAL
C           NOTATION, AND THE ARGUMENT K OF THESE FUNCTIONS MEANS
C           THE MODULUS.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DEFINITION
C           RES=INTEGRAL((A+B*T*T)/(SQRT((1+T*T)*(1+(CK*T)**2))*(1+T*T))
C           SUMMED OVER T FROM 0 TO INFINITY).
C           EVALUATION
C           LANDENS TRANSFORMATION IS USED FOR CALCULATION.
C           REFERENCE
C           R.BULIRSCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS
C           AND ELLIPTIC FUNCTIONS', HANDBOOK SERIES SPECIAL FUNCTIONS,
C           NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.
C
C     ..................................................................
C
      SUBROUTINE CEL2(RES,AK,A,B,IER)
      IER=0
      ARI=2.
      GEO=(0.5-AK)+0.5
      GEO=GEO+GEO*AK
      RES=A
      A1=A+B
      B0=B+B
      IF(GEO)1,2,6
    1 IER=1
    2 IF(B)3,8,4
    3 RES=-1.7E38                                                               0
      RETURN
    4 RES=1.7E38                                                                0
      RETURN
    5 GEO=GEO*AARI
    6 GEO=SQRT(GEO)
      GEO=GEO+GEO
      AARI=ARI
      ARI=ARI+GEO
      B0=B0+RES*GEO
      RES=A1
      B0=B0+B0
      A1=B0/ARI+A1
      IF(GEO/AARI-0.9999)5,7,7
    7 RES=A1/ARI
      RES=RES+0.5707963E0*RES
    8 RETURN
      END