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

C
C     ..................................................................
C
C        SUBROUTINE DCEL1
C
C        PURPOSE
C           CALCULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C
C        USAGE
C           CALL DCEL1(RES,AK,IER)
C
C        DESCRIPTION OF PARAMETERS
C           RES   - RESULT VALUE IN DOUBLE PRECISION
C           AK    - MODULUS (INPUT) IN DOUBLE PRECISION
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           THE RESULT IS SET TO 1.E75 IF ABS(AK) GE 1
C           FOR MODULUS AK AND COMPLEMENTARY MODULUS CK,
C           EQUATION AK*AK+CK*CK=1.D0 IS USED.
C           AK MUST BE IN THE RANGE -1 TO +1
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DEFINITION
C           CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CK*T)**2)), SUMMED
C           OVER T FROM 0 TO INFINITY).
C           EQUIVALENT ARE THE DEFINITIONS
C           CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED
C           OVER T FROM 0 TO PI/2),
C           CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER T
C           FROM 0 TO PI/2), WHERE K=SQRT(1.-CK*CK).
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 DCEL1(RES,AK,IER)
      DOUBLE PRECISION RES,AK,GEO,ARI,AARI
      IER=0
      ARI=2.D0
      GEO=(0.5D0-AK)+0.5D0
      GEO=GEO+GEO*AK
      RES=0.5D0
      IF(GEO)1,2,4
    1 IER=1
    2 RES=1.7D38                                                                0
      RETURN
    3 GEO=GEO*AARI
    4 GEO=DSQRT(GEO)
      GEO=GEO+GEO
      AARI=ARI
      ARI=ARI+GEO
      RES=RES+RES
      IF(GEO/AARI-0.999999995D0)3,5,5
    5 RES=RES/ARI*6.2831853071795865D0
      RETURN
      END