]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |