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

C
C     ..................................................................
C
C        SUBROUTINE ELI1
C
C        PURPOSE
C           COMPUTES THE ELLIPTIC INTEGRAL OF FIRST KIND
C
C        USAGE
C           CALL ELI1(RES,X,CK)
C
C        DESCRIPTION OF PARAMETERS
C           RES   - RESULT VALUE
C           X     - UPPER INTEGRATION BOUND (ARGUMENT OF ELLIPTIC
C                   INTEGRAL OF FIRST KIND)
C           CK    - COMPLEMENTARY MODULUS
C
C        REMARKS
C           MODULUS K = SQRT(1.-CK*CK).
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DEFINITION
C           RES=INTEGRAL(1/SQRT((1+T*T)*(1+(CK*T)**2)), SUMMED
C           OVER T FROM 0 TO X).
C           EQUIVALENT ARE THE DEFINITIONS
C           RES=INTEGRAL(1/(COS(T)*SQRT(1+(CK*TAN(T))**2)), SUMMED
C           OVER T FROM 0 TO ATAN(X)),
C           RES=INTEGRAL(1/SQRT(1-(K*SIN(T))**2), SUMMED OVER
C           T FROM 0 TO ATAN(X)).
C           EVALUATION
C           LANDENS TRANSFORMATION IS USED FOR CALCULATION.
C           REFERENCE
C           R. BULIRSCH, NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS AND
C                  ELLIPTIC FUNCTIONS.
C                  HANDBOOK SERIES OF SPECIAL FUNCTIONS
C                  NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.
C
C     ..................................................................
C
      SUBROUTINE ELI1(RES,X,CK)
C
      IF(X)2,1,2
    1 RES=0.
      RETURN
    2 IF(CK)4,3,4
    3 RES=ALOG(ABS(X)+SQRT(1.+X*X))
      GOTO 13
    4 ANGLE=ABS(1./X)
      GEO=ABS(CK)
      ARI=1.
      PIM=0.
    5 SQGEO=ARI*GEO
      AARI=ARI
      ARI=GEO+ARI
      ANGLE=-SQGEO/ANGLE+ANGLE
      SQGEO=SQRT(SQGEO)
      IF(ANGLE)7,6,7
C     REPLACE 0 BY SMALL VALUE
    6 ANGLE=SQGEO*1.E-8
    7 TEST=AARI*1.E-4
      IF(ABS(AARI-GEO)-TEST)10,10,8
    8 GEO=SQGEO+SQGEO
      PIM=PIM+PIM
      IF(ANGLE)9,5,5
    9 PIM=PIM+3.1415927
      GOTO 5
   10 IF(ANGLE)11,12,12
   11 PIM=PIM+3.1415927
   12 RES=(ATAN(ARI/ANGLE)+PIM)/ARI
   13 IF(X)14,15,15
   14 RES=-RES
   15 RETURN
      END