MINICERN/mathlib/gen/c/reli164.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:10  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DELI1(X,AKP)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15 #endif
  16 #if !defined(CERNLIB_DOUBLE)
  17       FUNCTION RELI1(X,AKP)
  18 C
  19 #endif
  20 C
  21 C     Translation of Algol procedure el1(x,kc) in
  22 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  23 C      Elliptic Functions,  Numer. Math. 7 (1965) 78-90
  24 C
  25       PARAMETER (ID = 16)
  26       PARAMETER (PI = 3.14159 26535 89793 24D0)
  27       PARAMETER (Z10 = 10)
  28       PARAMETER (CA = Z10**(-ID/2), CB = Z10**(-(ID+2)))
  29
  30       IF(X .EQ. 0) THEN
  31        H=0
  32       ELSEIF(AKP .EQ. 0) THEN
  33 #if defined(CERNLIB_DOUBLE)
  34        H=DASINH(X)
  35 #endif
  36 #if !defined(CERNLIB_DOUBLE)
  37        H=ASINH(X)
  38 #endif
  39       ELSE
  40        Y=ABS(1/X)
  41        YKP=ABS(AKP)
  42        XM=1
  43        L=0
  44     1  E=XM*YKP
  45        G=XM
  46        XM=YKP+XM
  47        Y=-E/Y+Y
  48        IF(Y .EQ. 0) Y=SQRT(E)*CB
  49        IF(ABS(G-YKP) .GT. CA*G) THEN
  50         YKP=2*SQRT(E)
  51         L=2*L
  52         IF(Y .LT. 0) L=L+1
  53         GO TO 1
  54        ENDIF
  55        IF(Y .LT. 0) L=L+1
  56        E=(ATAN(XM/Y)+PI*L)/XM
  57        IF(X .LT. 0) E=-E
  58        H=E
  59       ENDIF
  60 #if defined(CERNLIB_DOUBLE)
  61       DELI1=H
  62 #endif
  63 #if !defined(CERNLIB_DOUBLE)
  64       RELI1=H
  65 #endif
  66       RETURN
  67       END