MINICERN/mathlib/gen/c/reli264.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:11  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DELI2(X,AKP,A,B,MODE)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME
  16       PARAMETER(NAME='RELI2/DELI2')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION RELI2(X,AKP,A,B,MODE)
  20 C
  21       CHARACTER*(*) NAME
  22       PARAMETER(NAME='RELI2')
  23 #endif
  24 C
  25 C     Translation of Algol procedure el2(x,kc,a,b) in
  26 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  27 C      Elliptic Functions,  Numer. Math. 7 (1965) 78-90,
  28 C      extended for (k sub c)**2 < 0  (MODE = -1)
  29 C
  30       PARAMETER (ID = 16)
  31       PARAMETER (PI = 3.14159 26535 89793 24D0)
  32       PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2)
  33       PARAMETER (CA = Z10**(-ID/2), CB = Z10**(-(ID+2)))
  34
  35       CHARACTER*80 ERRTXT
  36
  37       IF(X .EQ. 0 .OR. AKP .EQ. 0) THEN
  38 #if defined(CERNLIB_DOUBLE)
  39        H=B*DASINH(X)+(A-B)*X/SQRT(1+X**2)
  40 #endif
  41 #if !defined(CERNLIB_DOUBLE)
  42        H=B*ASINH(X)+(A-B)*X/SQRT(1+X**2)
  43 #endif
  44       ELSE
  45        IF(MODE .EQ. 1) THEN
  46         XX=X
  47         BKP=AKP
  48         AA=A
  49         BB=B
  50        ELSEIF(MODE .EQ. -1) THEN
  51         W=1-(AKP*X)**2
  52         IF(W .GT. 0) THEN
  53          W1=1/SQRT(1+AKP**2)
  54          XX=X/(W1*SQRT(W))
  55          BKP=AKP*W1
  56          AA=A*W1
  57          BB=(B+A*AKP**2)*W1**3
  58         ELSE
  59          H=0
  60          WRITE(ERRTXT,101) X,AKP
  61          CALL MTLPRT(NAME,'C346.1',ERRTXT)
  62          GO TO 9
  63         ENDIF
  64        ELSE
  65         H=0
  66         WRITE(ERRTXT,102) MODE
  67         CALL MTLPRT(NAME,'C346.2',ERRTXT)
  68         GO TO 9
  69        ENDIF
  70        C=XX**2
  71        D=1+C
  72        PP=SQRT((1+BKP**2*C)/D)
  73        D=XX/D
  74        C=D/(2*PP)
  75        Z=AA-BB
  76        XI=AA
  77        AA=HF*(AA+BB)
  78        Y=ABS(1/XX)
  79        F=0
  80        L=0
  81        XM=1
  82        YKP=ABS(BKP)
  83     1  BB=XI*YKP+BB
  84        E=XM*YKP
  85        G=E/PP
  86        D=F*G+D
  87        F=C
  88        XI=AA
  89        PP=G+PP
  90        C=HF*(D/PP+C)
  91        G=XM
  92        XM=YKP+XM
  93        AA=HF*(BB/XM+AA)
  94        Y=Y-E/Y
  95        IF(Y .EQ. 0) Y=SQRT(E)*CB
  96        IF(ABS(G-YKP) .GT. CA*G) THEN
  97         YKP=2*SQRT(E)
  98         L=2*L
  99         IF(Y .LT. 0) L=L+1
 100         GO TO 1
 101        ENDIF
 102        IF(Y .LT. 0) L=L+1
 103        E=(ATAN(XM/Y)+PI*L)*AA/XM
 104        IF(XX .LT. 0) E=-E
 105        H=E+C*Z
 106       ENDIF
 107 #if defined(CERNLIB_DOUBLE)
 108     9 DELI2=H
 109 #endif
 110 #if !defined(CERNLIB_DOUBLE)
 111     9 RELI2=H
 112 #endif
 113       RETURN
 114   101 FORMAT('X = ',1P,D15.8,' AKP = ',D15.8,' ILLEGAL',
 115      1       '[(AKP * X)**2 >= 1]')
 116   102 FORMAT('MODE = ',I5,' ILLEGAL')
 117       END