]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.2 1997/12/15 16:18:35 mclareni | |
6 | * Changes for the Portland Group f77 compiler inside cpp define CERNLIB_QFPGF77 | |
7 | * | |
8 | * Revision 1.1.1.1 1996/04/01 15:02:01 mclareni | |
9 | * Mathlib gen | |
10 | * | |
11 | * | |
12 | #include "gen/pilot.h" | |
13 | #if defined(CERNLIB_DOUBLE) | |
14 | SUBROUTINE WELFUN(W,AK2,SN,CN,DN) | |
15 | C | |
16 | #include "gen/imp64.inc" | |
17 | C | |
18 | CHARACTER*(*) NAME | |
19 | PARAMETER(NAME='CELFUN/WELFUN') | |
20 | #endif | |
21 | #if !defined(CERNLIB_DOUBLE) | |
22 | SUBROUTINE CELFUN(W,AK2,SN,CN,DN) | |
23 | C | |
24 | CHARACTER*(*) NAME | |
25 | PARAMETER(NAME='CELFUN') | |
26 | #endif | |
27 | C | |
28 | C Jacobian Elliptic Functions SN(W,M),CN(W,M),DN(W,M) for | |
29 | C complex argument w = u+i*v. | |
30 | C Iterates for parameter m = k**2 or mc= 1-m if mc < m to obtain | |
31 | C fastest convergence and uses finally Jacobi's | |
32 | C imaginary transformation to go back to sn, cn, dn of m. | |
33 | C | |
34 | #include "gen/defc64.inc" | |
35 | + W,SN,CN,DN,I | |
36 | DIMENSION C(4) | |
37 | ||
38 | C MACHINE-DEPENDENT: EPS1=2**-(MB/2), EPS2=2**-(MB+3) | |
39 | C Where M = Number of bits in mantissa | |
40 | ||
41 | PARAMETER (MB = 64) | |
42 | PARAMETER (Z0 = 0, Z1 = 1, Z2 = 2, HF = Z1/2, QU = Z1/4) | |
43 | #if !defined(CERNLIB_CRAY) | |
44 | PARAMETER (I = (0D0,1D0)) | |
45 | #endif | |
46 | #if defined(CERNLIB_CRAY) | |
47 | PARAMETER (I = (0E0,1E0)) | |
48 | #endif | |
49 | PARAMETER (PI = 3.14159 26535 89793 24D0, PIH = PI/2) | |
50 | PARAMETER (EPS1 = Z2**(-MB/2), EPS2 = Z2**(-(MB+3))) | |
51 | ||
52 | CHARACTER*80 ERRTXT | |
53 | ||
54 | #if defined(CERNLIB_QFPGF77) | |
55 | DATA AM0 /-1D20/ | |
56 | ||
57 | SAVE AM0,C,A,L,BIGK | |
58 | #endif | |
59 | ||
60 | #if defined(CERNLIB_QF2C) | |
61 | #include "gen/gcmpfun.inc" | |
62 | #endif | |
63 | ||
64 | #if ! defined(CERNLIB_QFPGF77) | |
65 | DATA AM0 /-1D20/ | |
66 | ||
67 | SAVE AM0,C,A,L,BIGK | |
68 | #endif | |
69 | ||
70 | #if !defined(CERNLIB_QF2C) | |
71 | #include "gen/gcmpfun.inc" | |
72 | #endif | |
73 | ||
74 | AM=ABS(AK2) | |
75 | IF(AM .GT. 1) THEN | |
76 | WRITE(ERRTXT,101) AM | |
77 | CALL MTLPRT(NAME,'C320.1',ERRTXT) | |
78 | RETURN | |
79 | ENDIF | |
80 | IF(AM .LE. HF) THEN | |
81 | U=W | |
82 | V=GIMAG(W) | |
83 | IF(AM .EQ. AM0) GO TO 1 | |
84 | XK2=AM | |
85 | B=SQRT(1-AM) | |
86 | ELSE | |
87 | U=GIMAG(W) | |
88 | V=-W | |
89 | IF(AM .EQ. AM0) GO TO 1 | |
90 | XK2=1-AM | |
91 | B=SQRT(AM) | |
92 | ENDIF | |
93 | ||
94 | AM0=AM | |
95 | A=1 | |
96 | L=4 | |
97 | C(L)=QU*XK2 | |
98 | ||
99 | C Gaussian arithmetic-geometric mean. Skipped if previous M. | |
100 | ||
101 | 2 IF(C(L) .GE. EPS1) THEN | |
102 | L=L-1 | |
103 | C(L)=(QU*(A-B))**2 | |
104 | A1=HF*(A+B) | |
105 | B=SQRT(A*B) | |
106 | A=A1 | |
107 | GO TO 2 | |
108 | ENDIF | |
109 | A=HF*(A+B) | |
110 | BIGK=PIH/A | |
111 | ||
112 | C Descending Landen-Gauss Trafo for real U | |
113 | ||
114 | 1 X=SIN(A*U) | |
115 | IF(V .NE. 0 .OR. X .NE. 0) THEN | |
116 | IF(X .NE. 0) THEN | |
117 | X=A/X | |
118 | DO 3 J = L,4 | |
119 | X1=C(J)/X | |
120 | 3 X=X1+X | |
121 | SU=1/X | |
122 | DU=1-2*X1*SU | |
123 | CU=SIGN(SQRT(ABS(1-SU**2)),BIGK-ABS(U)) | |
124 | ENDIF | |
125 | IF(V .NE. 0) THEN | |
126 | Y=A/SINH(A*V) | |
127 | DO 4 J = L,4 | |
128 | Y1=C(J)/Y | |
129 | Y=Y-Y1 | |
130 | IF(Y .EQ. 0) Y=EPS2 | |
131 | 4 CONTINUE | |
132 | SV=1/Y | |
133 | Y1=2*Y1*SV | |
134 | DV=1+Y1 | |
135 | CV=SIGN(SQRT((SV+Y)*SV),Y1) | |
136 | ||
137 | C Evaluation of complex sn, cn, dn from real and imaginary ones and | |
138 | C Jacobi's imaginary argument trafo if necessary | |
139 | ||
140 | IF(U .NE. 0) THEN | |
141 | SS=-SU*SV | |
142 | D=1/(XK2*SS**2+1) | |
143 | DV=DV*D | |
144 | CU=CU*D | |
145 | CC=CU*CV | |
146 | DD=DU*DV | |
147 | DN=GCMPLX(DD,XK2*SS*CC) | |
148 | CN=GCMPLX(CC,SS*DD) | |
149 | SN=GCMPLX(SU*CV*DV,SV*CU*DU) | |
150 | ELSE | |
151 | SN=GCMPLX(Z0,SV) | |
152 | CN=CV | |
153 | DN=DV | |
154 | ENDIF | |
155 | ELSE | |
156 | SN=SU | |
157 | CN=CU | |
158 | DN=DU | |
159 | ENDIF | |
160 | IF(AM .GT. HF) THEN | |
161 | CN=1/CN | |
162 | DN=DN*CN | |
163 | SN=I*SN*CN | |
164 | ENDIF | |
165 | ELSE | |
166 | SN=0 | |
167 | CN=1 | |
168 | DN=1 | |
169 | ENDIF | |
170 | RETURN | |
171 | 101 FORMAT('MODULUS AK2 = ',1P,E15.6,' OUT OF RANGE') | |
172 | END |