MINICERN/mathlib/gen/c/religc64.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 DELIGC(AKP,P,A,B)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME
  16       PARAMETER(NAME='RELIGC/DELIGC')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION RELIGC(AKP,P,A,B)
  20 C
  21       CHARACTER*(*) NAME
  22       PARAMETER(NAME='RELIGC')
  23 #endif
  24 C
  25 C     Translation of Algol procedure cel(kc,p,a,b) in
  26 C      R. BULIRSCH Numerical Calculation of Elliptic Integrals and
  27 C      Elliptic Functions III.,  Numer. Math. 13 (1969) 305-315
  28 C
  29       PARAMETER (ID = 16)
  30       PARAMETER (PI = 3.14159 26535 89793 24D0, PIH = PI/2)
  31       PARAMETER (Z10 = 10)
  32       PARAMETER (CA = Z10**(-ID/2))
  33
  34       IF(AKP .EQ. 0) THEN
  35        H=0
  36        CALL MTLPRT(NAME,'C347.4','AKP = 0')
  37       ELSE
  38        PP=P
  39        AA=A
  40        BB=B
  41        YKP=ABS(AKP)
  42        E=YKP
  43        XM=1
  44        IF(PP .GT. 0) THEN
  45         PP=SQRT(PP)
  46         BB=BB/PP
  47        ELSE
  48         F=YKP**2
  49         Q=1-F
  50         G=1-PP
  51         F=F-PP
  52         Q=(BB-AA*PP)*Q
  53         PP=SQRT(F/G)
  54         AA=(AA-BB)/G
  55         BB=-Q/(G**2*PP)+AA*PP
  56        ENDIF
  57     1  F=AA
  58        AA=BB/PP+AA
  59        G=E/PP
  60        BB=2*(F*G+BB)
  61        PP=G+PP
  62        G=XM
  63        XM=YKP+XM
  64        IF(ABS(G-YKP) .GT. G*CA) THEN
  65         YKP=2*SQRT(E)
  66         E=YKP*XM
  67         GO TO 1
  68        ENDIF
  69        H=PIH*(AA*XM+BB)/(XM*(XM+PP))
  70       ENDIF
  71 #if defined(CERNLIB_DOUBLE)
  72       DELIGC=H
  73 #endif
  74 #if !defined(CERNLIB_DOUBLE)
  75       RELIGC=H
  76 #endif
  77       RETURN
  78       END