*
* $Id$
*
* $Log$
* Revision 1.1  1996/04/01 15:03:11  mclareni
* Initial revision
*
*
#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
