* * $Id: elin1.F,v 1.1.1.1 1996/04/01 15:03:12 mclareni Exp $ * * $Log: elin1.F,v $ * Revision 1.1.1.1 1996/04/01 15:03:12 mclareni * Mathlib gen * * #include "sys/CERNLIB_machine.h" #include "_gen/pilot.h" FUNCTION ELIN1(X,M) C C-----COMPUTES ELLIPTIC INTEGRAL 1ST KIND F(X,M) BY LANDEN TRANSFORM. C BASED ON R.BULIRSCH"S PROGRAM EL1 (NUM. MATH.VOL.7,P.78, (1965)) C BUT PARAMETER M=K**2 (K=MODULUS) IS USED INSTEAD OF KC=SQRT(1-M). C IMPROVEMENTS: C 1. CASE F(X,1) IS ALSO ADMITTED (= ATANH(X) ) C 2. PRINTED ERROR MESSAGE AT SINGULARITIES F(1,1), F(-1,1) C 3. IF K IS UNCHANGED IN SUCCESSIVE SUBROUTINE CALLS MODULUS TANSFOR- C MATIONS ARE NOT REPEATED, F(X,M) COMPUTED IN ABOUT HALF THE TIME. REAL M DIMENSION E(50) C THE FOLLOWING CONSTANTS ARE MACHINE-DEPENDANT. EPSI=2**-(M/2) AND C EPSI2=2**-(M+3) WHERE M=NUMBER OF MANTISSA BITS. #if defined(CERNLIB_CDC) DATA CLAST/2./,EPSI/4.2E-8/,EPSI2/4E-16/ #endif #if !defined(CERNLIB_CDC) DATA CLAST/2./,EPSI/1.5E-5/,EPSI2/3E-11/ #endif ELIN1= 0. IF(X.EQ.0.) RETURN C=1.-ABS(M) IF(C.EQ.0.) GO TO 6 IF(C.EQ.CLAST) GO TO 2 CLAST=C C-----ARITHMETICO-GEOMETRIC SCALE.(SKIPPED IF M IS UNCHANGED) A=1. B=SQRT(C) DO 1 I=1,50 E(I)=A*B D=(A-B)/A A=A+B IF(D.LE.EPSI) GO TO 2 1 B=SQRT(E(I))*2. C-----COMPUTATION OF ELIN1= F(X,M) BY LANDEN TRANSFORMATION. 2 Y=SQRT(1.-X**2)/ABS(X) DO 4 J=1,I ELIN1=2.*ELIN1 IF(Y.LT.0.) ELIN1=ELIN1+3.141592653589793 IF(Y.EQ.0.) Y=EPSI2 4 Y=Y-E(J)/Y ELIN1= SIGN((ATAN2(A,Y)+ELIN1)/A,X) RETURN C-----SPECIAL CASE M=1,YIELDS F(X,M)= ATANH(X) 6 IF(ABS(X).GE.1.) WRITE(6,1000)X,M ELIN1=0.5*LOG((1.+X)/(1.-X)) RETURN 1000 FORMAT('0*****ILLEGAL ARGUMENT IN ELLIPTIC INTEGRAL ELIN1(X,M).', * 'X=',F5.2,',M=',F5.2) END