This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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