C
C     ..................................................................
C
C        SUBROUTINE DJELF
C
C        PURPOSE
C           COMPUTES THE THREE JACOBIAN ELLIPTIC FUNCTIONS SN, CN, DN.
C
C        USAGE
C           CALL DJELF(SN,CN,DN,X,SCK)
C
C        DESCRIPTION OF PARAMETERS
C           SN    - RESULT VALUE SN(X) IN DOUBLE PRECISION
C           CN    - RESULT VALUE CN(X) IN DOUBLE PRECISION
C           DN    - RESULT VALUE DN(X) IN DOUBLE PRECISION
C           X     - DOUBLE PRECISION ARGUMENT OF JACOBIAN ELLIPTIC
C                   FUNCTIONS
C           SCK   - SQUARE OF COMPLEMENTARY MODULUS IN DOUBLE PRECISION
C
C        REMARKS
C           NONE
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           DEFINITION
C           X=INTEGRAL(1/SQRT((1-T*T)*(1-(K*T)**2)), SUMMED OVER
C           T FROM 0 TO SN), WHERE K=SQRT(1-SCK).
C           SN*SN + CN*CN = 1
C           (K*SN)**2 + DN**2 = 1.
C           EVALUATION
C           CALCULATION IS DONE USING THE PROCESS OF THE ARITHMETIC
C           GEOMETRIC MEAN TOGETHER WITH GAUSS DESCENDING TRANSFORMATION
C           BEFORE INVERSION OF THE INTEGRAL TAKES PLACE.
C           REFERENCE
C           R. BULIRSCH, NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS AND
C                  ELLIPTIC FUNCTIOMS.
C                  HANDBOOK SERIES OF SPECIAL FUNCTIONS
C                  NUMERISCHE MATHEMATIK VOL. 7, 1965, PP. 78-90.
C
C     ..................................................................
C
      SUBROUTINE DJELF(SN,CN,DN,X,SCK)
C
      DIMENSION ARI(12),GEO(12)
      DOUBLE PRECISION SN,CN,DN,X,SCK,ARI,GEO,CM,Y,A,B,C,D
C
C        TEST MODULUS
C
      CM=SCK
      Y=X
      IF(SCK)3,1,4
    1 D=DEXP(X)
      A=1.D0/D
      B=A+D
      CN=2.D0/B
      DN=CN
      A=(D-A)/2.D0
      SN=A*CN
C        DEGENERATE CASE SCK=0 GIVES RESULTS
C           CN X = DN X = 1/COSH X
C           SN X = TANH X
    2 RETURN
C
C        JACOBIS MODULUS TRANSFORMATION
C
    3 D=1.D0-SCK
      CM=-SCK/D
      D=DSQRT(D)
      Y=D*X
    4 A=1.D0
      DN=1.D0
      DO 6 I=1,12
      L=I
      ARI(I)=A
      CM=DSQRT(CM)
      GEO(I)=CM
      C=(A+CM)*.5D0
      IF(DABS(A-CM)-1.D-9*A)7,7,5
    5 CM=A*CM
    6 A=C
C
C        START BACKWARD RECURSION
C
    7 Y=C*Y
      SN=DSIN(Y)
      CN=DCOS(Y)
      IF(SN)8,13,8
    8 A=CN/SN
      C=A*C
      DO 9 I=1,L
      K=L-I+1
      B=ARI(K)
      A=C*A
      C=DN*C
      DN=(GEO(K)+A)/(B+A)
    9 A=C/B
      A=1.D0/DSQRT(C*C+1.D0)
      IF(SN)10,11,11
   10 SN=-A
      GOTO 12
   11 SN=A
   12 CN=C*SN
   13 IF(SCK)14,2,2
   14 A=DN
      DN=CN
      CN=A
      SN=SN/D
      RETURN
      END
