C C This file contain of algorithms for the calculation of elliptic C integrals and the Jacobian elliptic functions by papers R.Bulirsch: C Numerical Calculation of Elliptic Integrals and Elliptic Functions, C Numerische Mathematik 7,78-90 (1965) and R.Bulirsch: Numerical C Calculation of Elliptic Integrals and Elliptic Functions III, C Numerische Mathematik 13,305-315 (1969) C C F.Hroch, May 1997 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the elliptic integral of the first kind. FUNCTION EL1(X,KC) DOUBLE PRECISION EL1,X,KC C Number of significant digits. INTEGER D PARAMETER( D = 16 ) Constants associated to D CA = 10**(-D/2) CB = 10**(-(D + 2)) DOUBLE PRECISION CA,CB PARAMETER( CA = 1D-8, CB = 1D-18 ) C PI DOUBLE PRECISION PI PARAMETER( PI = 3.14159 26535 89793 23846 ) INTEGER L DOUBLE PRECISION E,G,M,Y,KC1,ASINH IF( X .EQ. 0.0 )THEN EL1 = 0.0 ELSE IF( KC .NE. 0.0 )THEN Y = ABS(1.0/X) KC1 = ABS(KC) M = 1.0 L = 0 00001 CONTINUE E = M*KC1 G = M M = KC1 + M Y = -E/Y + Y IF( Y .EQ. 0.0 ) Y = SQRT(E)*CB IF( ABS(G - KC1) .GT. CA*G )THEN KC1 = 2.0*SQRT(E) L = 2*L IF( Y .LT. 0.0 ) L = L + 1 GOTO 00001 ENDIF IF( Y .LT. 0.0 ) L = L + 1 E = (ATAN(M/Y) + PI*L)/M EL1 = SIGN(E,X) ELSE EL1 = ASINH(X) ENDIF ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the complete elliptic integral of the first kind. FUNCTION CEL1(KC) DOUBLE PRECISION CEL1,KC C Number of significant digits. INTEGER D PARAMETER( D = 16 ) Constant associated to D CA = 10**(-D/2) DOUBLE PRECISION CA PARAMETER( CA = 1D-8 ) C PI DOUBLE PRECISION PI PARAMETER( PI = 3.14159 26535 89793 23846 ) DOUBLE PRECISION H,M,KC1 KC1 = KC IF( KC1 .NE. 0.0 )THEN KC1 = ABS(KC1) M = 1.0 00001 CONTINUE H = M M = KC1 + M IF( ABS(H - KC1) .GT. CA*H )THEN KC1 = SQRT(H*KC1) M = 0.5*M GOTO 00001 ENDIF CEL1 = PI/M ELSE CEL1 = 1.0/CA ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the general elliptic integral of the second kind. C Yields for a*b >= 0 about D valid digits, which are significant C if |a| >= |b|. The c*z part of value "el2" gives the value C Jacobian Zeta function for z = k**2. FUNCTION EL2(X,KC,A,B) DOUBLE PRECISION EL2,X,KC,A,B C Number of significant digits. INTEGER DD PARAMETER( DD = 16 ) Constants associated to DD CA = 10**(-DD/2) CB = 10**(-(DD + 2)) DOUBLE PRECISION CA,CB PARAMETER( CA = 1D-8, CB = 1D-18 ) C PI DOUBLE PRECISION PI PARAMETER( PI = 3.14159 26535 89793 23846 ) INTEGER L DOUBLE PRECISION C,D,E,F,G,I,M,P,Y,Z,KC1,A1,B1 IF( X .EQ. 0.0 )THEN EL2 = 0.0 ELSE IF( KC .NE. 0.0 )THEN C = X**2 D = C + 1.0 P = SQRT((1.0 + C*KC**2)/D) D = X/D C = D/(2.0*P) Z = A - B I = A A1 = (A + B)/2.0 Y = ABS(1.0/X) F = 0.0 L = 0 M = 1.0 KC1 = ABS(KC) B1 = B 00001 CONTINUE B1 = I*KC1 + B1 E = M*KC1 G = E/P D = F*G + D F = C I = A1 P = G + P C = (D/P + C)/2.0 G = M M = KC1 + M A1 = (B1/M + A1)/2.0 Y = -E/Y + Y IF( Y .EQ. 0.0 ) Y = SQRT(E)*CB IF( ABS(G - KC1) .GT. CA*G )THEN KC1 = 2.0*SQRT(E) L = 2*L IF( Y .LT. 0.0 ) L = L + 1 GOTO 00001 ENDIF IF( Y .LT. 0.0 ) L = L + 1 E = (ATAN(M/Y) + PI*L)*A1/M IF( X .LT. 0.0 ) E = -E EL2 = E + C*Z ELSE C EL2 = -1 EL2 = 0.5*(X*(A + B) + 0.5*(A - B)*SIN(2.0*X)) ENDIF ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the general complete elliptic integral of the second kind. C Yields about D significant digit for a*b >= 0. FUNCTION CEL2(KC,A,B) DOUBLE PRECISION CEL2,KC,A,B C Number of significant digits. INTEGER D PARAMETER( D = 16 ) Constant associated to D CA = 10**(-D/2) DOUBLE PRECISION CA PARAMETER( CA = 1D-8 ) C PI/4 DOUBLE PRECISION PIFO PARAMETER( PIFO = 0.78539 81633 97448 30962 ) DOUBLE PRECISION C,M,M0,KC1,A1,B1 IF( KC .NE. 0.0 )THEN M = 1.0 C = A A1 = A + B KC1 = ABS(KC) B1 = B 00001 CONTINUE B1 = 2.0*(C*KC1 + B1) C = A1 M0 = M M = M + KC1 A1 = B1/M + A1 IF( ABS(M0 - KC1) .GT. CA*M0 )THEN KC1 = 2.0*SQRT(KC1*M0) GOTO 00001 ENDIF CEL2 = PIFO*A1/M ELSE CEL2 = -1 ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the values of the tree Jacobian elliptic functions. C NO ANY TESTS ! SUBROUTINE SNCNDN(XX,MCC,SN,CN,DN) DOUBLE PRECISION XX,MCC,SN,CN,DN C Number of significant digits. C INTEGER D C PARAMETER( D = 16 ) Constant associated to D CA = 10**(-D/2) DOUBLE PRECISION CA PARAMETER( CA = 1D-8 ) INTEGER I,L LOGICAL BO DOUBLE PRECISION A,B,C,D,M(0:12),N(0:12),X,MC X = XX MC = MCC IF( MC .NE. 0.0 )THEN BO = MC .LT. 0.0 IF( BO )THEN D = 1.0 - MC MC = - MC/D D = SQRT(D) X = X*D ENDIF DN = 1.0 A = 1.0 DO I = 0,12 L = I M(I) = A MC = SQRT(MC) N(I) = MC C = 0.5*(A + MC) IF( ABS(A - MC) .LE. CA*A ) GOTO 001 MC = A*MC A = C ENDDO 001 X = X*C SN = SIN(X) CN = COS(X) IF( SN .EQ. 0.0 ) GOTO 002 A = CN/SN C = C*A DO I = L,0,-1 B = M(I) A = A*C C = C*DN DN = (N(I) + A)/(A + B) A = C/B ENDDO A = 1.0/SQRT(1.0 + C**2) SN = SIGN(A,SN) CN = C*SN 002 IF( BO )THEN A = DN DN = CN CN = A SN = SN/D ENDIF ELSE D = EXP(X) A = 1.0/D B = A + D CN = 2.0/B DN = CN IF( ABS(X) .LT. 0.3 )THEN D = X**3 SN = CN*(D*(0.33333333333333333333 + + 3.96825396825D-4*D*X) + SIN(X)) ELSE SN = (D - A)/B ENDIF ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the incomplete elliptic integral of the third kind. C FUNCTION EL3(XX,KCC,PP) DOUBLE PRECISION EL3,XX,KCC,PP DOUBLE PRECISION LN2,PI,CA,CB PARAMETER( LN2 = 0.6931471805599453 ) PARAMETER( PI = 3.1415926535897932 ) C D = 16, CA = 10**(-D/2), CB = 10**(-(D + 2)), ND = D - 2 PARAMETER( CA = 1D-8, CB = 1D-18 ) INTEGER ND PARAMETER( ND = 14 ) INTEGER K,N,M,L LOGICAL BO,BK DOUBLE PRECISION AM,AP,C,D,DE,E,F,FA,G,H,HH,P1,PM,PZ,Q,R,S,T DOUBLE PRECISION U,V,W,Y,YE,Z,ZD,X,KC,P DOUBLE PRECISION RA(2:ND),RB(2:ND),RR(2:ND) DOUBLE PRECISION ASINH SAVE RA,RB,ZD DATA ZD/-1.0/ IF( XX .EQ. 0.0 )THEN EL3 = 0.0 RETURN ENDIF EL3 = -1 X = XX KC = KCC P = PP HH = X**2 F = P*HH IF( KC .EQ. 0.0 )THEN S = CA/(1.0 + ABS(X)) ELSE S = KC ENDIF T = S**2 PM = 0.5*T E = HH*T Z = ABS(F) R = ABS(P) H = 1.0 + HH IF( E.LT.0.1 .AND. Z.LT.0.1 .AND. T.LT.1.0 .AND. R.LT.1.0 )THEN IF( ZD .LT. 0.0 )THEN DO K = 2,ND RB(K) = 0.5/K RA(K) = 1.0 - RB(K) ENDDO ZD = 0.5/(ND + 1.0) ENDIF S = P + PM DO K = 2,ND RR(K) = S PM = PM*T*RA(K) S = S*P + PM ENDDO S = S*ZD U = S BO = .FALSE. DO K = ND,2,-1 U = U + (RR(K) - U)*RB(K) BO = .NOT. BO IF( BO )THEN V = -U ELSE V = U ENDIF S = S*HH + V ENDDO IF( BO ) S = -S U = 0.5*(U + 1.0) EL3 = (U - S*H)*SQRT(H)*X + U*ASINH(X) ELSE W = 1.0 + F IF( W .EQ. 0.0 ) RETURN IF( P .EQ. 0.0 )THEN P1 = CB/HH ELSE P1 = P ENDIF S = ABS(S) Y = ABS(X) G = P1 - 1.0 IF( G .EQ. 0.0 ) G = CB F = P1 - T IF( F .EQ. 0.0 ) F = CB*T AM = 1.0 - T AP = 1.0 + E R = P1*H FA = G/F/P1 BO = FA .GT. 0.0 FA = ABS(FA) PZ = ABS(G*F) DE = SQRT(PZ) Q = SQRT(ABS(P1)) IF( PM .GT. 0.5 ) PM = 0.5 PM = P1 - PM IF( PM .GE. 0.0 )THEN U = SQRT(R*AP) V = Y*DE IF( G .LT. 0.0 ) V = - V D = 1.0/Q C = 1.0 ELSE U = SQRT(H*AP*PZ) YE = Y*Q V = AM*YE Q = - DE/G D = - AM/DE C = 0.0 PZ = AP - R ENDIF IF( BO )THEN R = V/U Z = 1.0 K = 1 IF( PM .LT. 0.0 )THEN H = Y*SQRT(H/AP/FA) H = 1.0/H - H Z = H - 2.0*R R = 2.0 + R*H IF( R .EQ. 0.0 ) R = CB IF( Z .EQ. 0.0 ) Z = H*CB R = R/Z Z = R W = PZ ENDIF U = U/W V = V/W ELSE T = U + ABS(V) BK = .TRUE. IF( P1 .LT. 0.0 )THEN DE = V/PZ YE = U*YE YE = 2.0*YE U = T/PZ V = (- F - G*E)/T T = PZ*ABS(W) Z = (HH*R*F - G*AP + YE)/T YE = YE/T ELSE DE = V/W YE = 0.0 U = (E + P1)/T V = T/W Z = 1.0 ENDIF IF( S .GT. 1.0 )THEN H = U U = V V = H ENDIF ENDIF Y = 1.0/Y E = S N = 1 M = 0 L = 0 T = 1.0 001 Y = Y - E/Y IF( Y .EQ. 0.0 ) Y = CB*SQRT(E) F = C C = D/Q + C G = E/Q D = 2.0*(F*G + D) Q = G + Q G = T T = S + T N = 2*N M = 2*M IF( BO )THEN IF( Z .LT. 0.0 ) M = M + K K = R/ABS(R) H = E/(U**2 + V**2) U = U*(1.0 + H) V = V*(1.0 - H) ELSE R = U/V H = Z*R Z = Z*H HH = E/V IF( BK )THEN DE = DE/U YE = YE*(H + 1.0/H) + DE*(1.0 + R) DE = DE*(U - HH) BK = ABS(YE) .LT. 1.0 ELSE CALL CRACK(Z,Z,K) M = M + K ENDIF ENDIF IF( ABS(G - S) .GT. CA*G )THEN IF( BO )THEN G = 0.5*(1.0/R - R) HH = U + V*G H = G*U - V IF( HH .EQ. 0.0 ) HH = U*CB IF( H .EQ. 0.0 ) H = V*CB Z = R*H R = HH/H ELSE U = U + E/U V = V + HH ENDIF S = 2.0*SQRT(E) E = S*T L = 2*L IF( Y .LT. 0.0 ) L = L + 1 GOTO 001 ENDIF IF( Y .LT. 0.0 ) L = L + 1 E = ATAN(T/Y) + PI*L E = E*(C*T + D)/(T*(T + Q)) IF( BO )THEN H = V/(T + U) Z = 1.0 - R*H H = R + H IF( Z .EQ. 0.0 ) Z = CB IF( Z .LT. 0.0 ) M = M + H/ABS(H) S = ATAN(H/Z) + M*PI ELSE IF( BK )THEN S = ASINH(YE) ELSE S = LOG(ABS(Z)) + M*LN2 ENDIF S = 0.5*S ENDIF E = (E + SQRT(FA)*S)/N EL3 = SIGN(E,X) ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Evaluates the general complete elliptic integral. FUNCTION CEL(KCC,PP,AA,BB) DOUBLE PRECISION CEL,KCC,PP,AA,BB C Number of significant digits. C INTEGER D C PARAMETER( D = 16 ) Constant associated to D CA = 10**(-D/2) DOUBLE PRECISION CA PARAMETER( CA = 1D-8 ) DOUBLE PRECISION KC,P,A,B,E,F,G,M,Q KC = KCC P = PP A = AA B = BB IF( KC .NE. 0.0 )THEN KC = ABS(KC) E = KC M = 1.0 IF( P .GT. 0.0 )THEN P = SQRT(P) B = B/P ELSE F = KC**2 Q = 1.0 - F G = 1.0 - P F = F - P Q = Q*(B - A*P) P = SQRT(F/G) A = (A - B)/G B = - Q/(P*G**2) + A*P ENDIF 001 F = A A = B/P + A G = E/P B = 2.0*(F*G + B) P = P + G G = M M = M + KC IF( ABS(G - KC) .GT. G*CA )THEN KC = 2.0*SQRT(E) E = KC*M GOTO 001 ENDIF CEL = 1.570796326794897*(A*M + B)/(M*(M + P)) ELSE CEL = 0.0 ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION ASINH(X) DOUBLE PRECISION ASINH,X C Arcsinh(x) = sign(x)*ln(|x| + sqrt(1.0 + x**2)) ASINH = SIGN(LOG(ABS(X) + SQRT(X**2 + 1.0)),X) END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC