MINICERN/mathlib/gen/c/reli1c64.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 DELI1C(AKP)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME1,NAME2,NAME3
  16       PARAMETER(NAME1='RELI1C/DELI1C')
  17       PARAMETER(NAME2='RELI2C/DELI2C')
  18       PARAMETER(NAME3='RELI3C/DELI3C')
  19 #endif
  20 #if !defined(CERNLIB_DOUBLE)
  21       FUNCTION RELI1C(AKP)
  22 C
  23       CHARACTER*(*) NAME1,NAME2,NAME3
  24       PARAMETER(NAME1='RELI1C')
  25       PARAMETER(NAME2='RELI2C')
  26       PARAMETER(NAME3='RELI3C')
  27 #endif
  28 C
  29 C     Translation of Algol procedures cel1(kc), cel2(kc,a,b) in
  30 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  31 C      Elliptic Functions,  Numer. Math. 7 (1965) 78-90
  32 C     and of Algol procedure cel3(kc,m,p) in
  33 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  34 C      Elliptic Functions II.,  Numer. Math. 7 (1965) 353-354
  35 C
  36       PARAMETER (ID = 16)
  37       PARAMETER (PI = 3.14159 26535 89793 24D0)
  38       PARAMETER (PIH = PI/2, PIQ = PI/4)
  39       PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2)
  40       PARAMETER (CA = Z10**(-ID/2))
  41
  42       IF(AKP .EQ. 0) THEN
  43        H=0
  44        CALL MTLPRT(NAME1,'C347.1','AKP = 0')
  45       ELSE
  46        YKP=ABS(AKP)
  47        XM=1
  48     1  G=XM
  49        XM=YKP+XM
  50        IF(ABS(G-YKP) .GT. CA*G) THEN
  51         YKP=SQRT(G*YKP)
  52         XM=HF*XM
  53         GO TO 1
  54        ENDIF
  55        H=PI/XM
  56       ENDIF
  57       GO TO 9
  58
  59 #if defined(CERNLIB_DOUBLE)
  60       ENTRY DELI2C(AKP,A,B)
  61 #endif
  62 #if !defined(CERNLIB_DOUBLE)
  63       ENTRY RELI2C(AKP,A,B)
  64 #endif
  65
  66       IF(AKP .EQ. 0) THEN
  67        IF(B .EQ. 0) THEN
  68         H=A
  69        ELSE
  70         H=0
  71         CALL MTLPRT(NAME2,'C347.2','AKP = 0, B NE 0')
  72        ENDIF
  73       ELSE
  74        AA=A
  75        BB=B
  76        XM=1
  77        C=AA
  78        AA=BB+AA
  79        YKP=ABS(AKP)
  80     2  BB=2*(C*YKP+BB)
  81        C=AA
  82        XM0=XM
  83        XM=YKP+XM
  84        AA=BB/XM+AA
  85        IF(ABS(XM0-YKP) .GT. CA*XM0) THEN
  86         YKP=2*SQRT(YKP*XM0)
  87         GO TO 2
  88        ENDIF
  89        H=PIQ*AA/XM
  90       ENDIF
  91       GO TO 9
  92
  93 #if defined(CERNLIB_DOUBLE)
  94       ENTRY DELI3C(AKP,AK2,P)
  95 #endif
  96 #if !defined(CERNLIB_DOUBLE)
  97       ENTRY RELI3C(AKP,AK2,P)
  98 #endif
  99
 100       IF(AKP*P .EQ. 0) THEN
 101        H=0
 102        CALL MTLPRT(NAME3,'C347.3','AKP * P = 0')
 103       ELSE
 104        YKP=ABS(AKP)
 105        E=YKP
 106        AM0=1
 107        PP=P
 108        IF(PP .GT. 0) THEN
 109         C=1
 110         PP=SQRT(PP)
 111         D=1/PP
 112        ELSE
 113         G=1-PP
 114         F=YKP**2-PP
 115         PP=SQRT(F/G)
 116         D=-AK2/(G*PP)
 117         C=0
 118        ENDIF
 119     3  F=C
 120        C=D/PP+C
 121        G=E/PP
 122        D=2*(F*G+D)
 123        PP=G+PP
 124        G=AM0
 125        AM0=YKP+AM0
 126        IF(ABS(G-YKP) .GT. CA*G) THEN
 127         YKP=2*SQRT(E)
 128         E=YKP*AM0
 129         GO TO 3
 130        ENDIF
 131        H=PIH*(C*AM0+D)/(AM0*(AM0+PP))
 132       ENDIF
 133 #if defined(CERNLIB_DOUBLE)
 134     9 DELI1C=H
 135 #endif
 136 #if !defined(CERNLIB_DOUBLE)
 137     9 RELI1C=H
 138 #endif
 139       RETURN
 140       END