]>
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 WELINT(Z,AKC,A,B) | |
12 | C | |
13 | #include "gen/imp64.inc" | |
14 | C | |
15 | CHARACTER*(*) NAME | |
16 | PARAMETER(NAME='CELINT/WELINT') | |
17 | COMPLEX*16 WELINT | |
18 | #endif | |
19 | #if !defined(CERNLIB_DOUBLE) | |
20 | FUNCTION CELINT(Z,AKC,A,B) | |
21 | C | |
22 | CHARACTER*(*) NAME | |
23 | PARAMETER(NAME='CELINT') | |
24 | COMPLEX CELINT | |
25 | #endif | |
26 | C | |
27 | C Translation of Algol procedure elco2(x,y,kc,a,b,u,v) in | |
28 | C R. BULIRSCH Numerical Calculation of Elliptic Integrals and | |
29 | C Elliptic Functions, Numer. Math. 7 (1965) 78-90 | |
30 | C | |
31 | #include "gen/defc64.inc" | |
32 | + Z,I,W,C0 | |
33 | DIMENSION D1(0:12),D2(0:12) | |
34 | ||
35 | PARAMETER (ID = 16) | |
36 | #if !defined(CERNLIB_CRAY) | |
37 | PARAMETER (I = (0D0,1D0)) | |
38 | #endif | |
39 | #if defined(CERNLIB_CRAY) | |
40 | PARAMETER (I = (0E0,1E0)) | |
41 | #endif | |
42 | PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2) | |
43 | PARAMETER (CC = Z10**(-ID), C0 = I-I) | |
44 | ||
45 | CHARACTER*80 ERRTXT | |
46 | ||
47 | X=Z | |
48 | Y=-I*Z | |
49 | IF(Z .EQ. C0 .OR. AKC .EQ. 0) THEN | |
50 | W=SQRT(1+Z**2) | |
51 | W=B*LOG(Z+W)+(A-B)*Z/W | |
52 | ELSEIF(X .LT. 0) THEN | |
53 | W=0 | |
54 | WRITE(ERRTXT,101) X | |
55 | CALL MTLPRT(NAME,'C348.1',ERRTXT) | |
56 | ELSE | |
57 | AA=A | |
58 | BB=B | |
59 | SY=SIGN(Z1,Y) | |
60 | Y=ABS(Y) | |
61 | C=X**2-Y**2 | |
62 | E2=2*X*Y | |
63 | D=AKC**2 | |
64 | XK=1-D | |
65 | E1=1+C | |
66 | F2=1/(E1**2+E2**2) | |
67 | F1=((1+D*C)*E1+D*E2**2)*F2 | |
68 | F2=-2*XK*X*Y*F2 | |
69 | DN1=SQRT(HF*(SQRT(F1**2+F2**2)+ABS(F1))) | |
70 | DN2=HF*F2/DN1 | |
71 | IF(F1 .LT. 0) THEN | |
72 | F1=DN1 | |
73 | DN1=-DN2 | |
74 | DN2=-F1 | |
75 | ENDIF | |
76 | IF(XK .LT. 0) THEN | |
77 | DN1=ABS(DN1) | |
78 | DN2=ABS(DN2) | |
79 | ENDIF | |
80 | C=1+DN1 | |
81 | F1=E1*C-E2*DN2 | |
82 | F2=E1*DN2+E2*C | |
83 | D2(0)=1/(F1**2+F2**2) | |
84 | D1(0)=(X*F1+Y*F2)*D2(0) | |
85 | D2(0)=(Y*F1-X*F2)*D2(0) | |
86 | H=AA-BB | |
87 | N=1 | |
88 | XM=1 | |
89 | F=1 | |
90 | D=1 | |
91 | YKC=ABS(AKC) | |
92 | E=AA | |
93 | AA=BB+AA | |
94 | L=4 | |
95 | 1 XM1=HF*(YKC+XM) | |
96 | XM2=XM1**2 | |
97 | XK=F*XK/(4*XM2) | |
98 | BB=E*YKC+BB | |
99 | E=AA | |
100 | F2=1/(C**2+DN2**2) | |
101 | F1=((YKC+XM*DN1)*C+XM*DN2**2)*F2 | |
102 | E1=F1/XM1 | |
103 | E2=XK*DN2*F2 | |
104 | DN1=SQRT(HF*(SQRT(E1**2+(2*E2)**2)+ABS(E1))) | |
105 | DN2=E2/DN1 | |
106 | F1=DN1*X-DN2*Y | |
107 | F2=DN1*Y+DN2*X | |
108 | X=ABS(F1) | |
109 | Y=ABS(F2) | |
110 | AA=BB/XM1+AA | |
111 | L=2*L | |
112 | C=1+DN1 | |
113 | D=HF*XK*D | |
114 | E1=1+(X**2-Y**2)*XM2 | |
115 | E2=2*X*Y*XM2 | |
116 | F1=C*E1-DN2*E2 | |
117 | F2=C*E2+DN2*E1 | |
118 | E1=D/(F1**2+F2**2) | |
119 | D1(N)=(X*F1+Y*F2)*E1 | |
120 | D2(N)=(Y*F1-X*F2)*E1 | |
121 | XK=XK**2 | |
122 | IF(XK .GT. CC) THEN | |
123 | YKC=SQRT(XM*YKC) | |
124 | F=XM2 | |
125 | XM=XM1 | |
126 | N=N+1 | |
127 | GO TO 1 | |
128 | ENDIF | |
129 | F2=0 | |
130 | F1=0 | |
131 | DO 2 J = N,0,-1 | |
132 | F1=D1(J)+F1 | |
133 | 2 F2=D2(J)+F2 | |
134 | X=XM1*X | |
135 | Y=XM1*Y | |
136 | C=X**2+Y**2 | |
137 | E2=1/(1+2*Y+C) | |
138 | E1=(1-C)*E2 | |
139 | E2=2*X*E2 | |
140 | D=AA/(XM1*L) | |
141 | W=D*ATAN2(E2,E1)+H*F1+I*SY*(H*F2-LOG(E1**2+E2**2)*HF*D) | |
142 | ENDIF | |
143 | #if defined(CERNLIB_DOUBLE) | |
144 | WELINT=W | |
145 | #endif | |
146 | #if !defined(CERNLIB_DOUBLE) | |
147 | CELINT=W | |
148 | #endif | |
149 | RETURN | |
150 | 101 FORMAT('RE(Z) = ',1P,D15.8,' < 0') | |
151 | END |