This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / bzejy64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:09  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if !defined(CERNLIB_DOUBLE)
11       SUBROUTINE RBZEJY(A,N,MODE,REL,X)
12 #endif
13 #if defined(CERNLIB_DOUBLE)
14       SUBROUTINE DBZEJY(A,N,MODE,REL,X)
15 #include "gen/imp64.inc"
16 #endif
17 C     Computes the first n positive (in the case Jo'(x) the first n
18 C     non-negative) zeros of the Bessel functions
19 C               Ja(x), Ya(x), Ja'(x), Ya'(x),
20 C     where a >= 0 and ' = d/dx.
21 C
22 C     Based on Algol procedures published in
23 C
24 C     N.M. TEMME, An algorithm with Algol 60 program for the compu-
25 C     tation of the zeros of ordinary Bessel functions and those of
26 C     their derivatives, J. Comput. Phys. 32 (1979) 270-279, and
27 C
28 C     N.M. TEMME, On the numerical evaluation of the ordinary Bessel
29 C     function of the second kind, J. Comput. Phys. 21 (1976) 343-350.
30
31       DIMENSION X(*),FCO(4),ICO(4),H(13)
32       CHARACTER NAME*(*)
33       CHARACTER*80 ERRTXT
34 #if !defined(CERNLIB_DOUBLE)
35       PARAMETER (NAME = 'RBZEJY')
36 #endif
37 #if defined(CERNLIB_DOUBLE)
38       PARAMETER (NAME = 'RBZEJY/DBZEJY')
39 #endif
40
41       PARAMETER (Z1 = 1, Z2 = 2, Z5 = 5, Z6 = 6, Z7 = 7)
42       PARAMETER (C1 = Z5/48, C2 = -Z5/36, C3 = -Z7/48)
43       PARAMETER (C4 = 5*Z7/288, C5 = Z1/64, C6 = Z6/5, C7 = 3*Z6/7)
44       PARAMETER (E1 = Z1/3, E2 = Z2/3, HALF = Z1/2)
45       PARAMETER (PI = 3.14159 26535 89793 24D0)
46       PARAMETER (PIH = PI/2, PI1 = 3*PI/8)
47
48       DATA FCO /-2.33811,-1.17371,-1.01879,-2.29444/, ICO /1,3,3,1/
49
50 #if !defined(CERNLIB_DOUBLE)
51       ENTRY BZEJY(A,N,MODE,REL,X)
52 #endif
53
54       IF(N .LE. 0) RETURN
55       IF(A .LT. 0) THEN
56        WRITE(ERRTXT,101) A
57        CALL MTLPRT(NAME,'C345.1',ERRTXT)
58        RETURN
59       ENDIF
60       NA=NINT(A)
61       A0=A-NA
62       H(1)=A**2
63       H(2)=4*H(1)
64       H(3)=H(2)**2
65       H(4)=H(3)+H(3)
66       H(5)=12*A+6
67       H(6)=H(5)-4
68       H(7)=3*A-8
69       H(8)=HALF*(A-HALF*ICO(MODE))
70       H(9)=COS(A0*PI)/(PI*REL)
71       H(10)=PIH*(A+HALF)
72       H(11)=HALF+A0
73       H(12)=H(11)*(HALF-A0)
74       H(13)=1
75       IF(A .GE. 3) H(13)=A**(-E2)
76       P1=0
77       Q1=0
78       IF(MODE .LT. 3) THEN
79        IF(N .GE. H(7)) THEN
80         P=7*H(2)-31
81         PP=H(2)-1
82         IF(1+P .NE. P) THEN
83          P1=PP*((1012*H(2)-14888)*H(2)+71476)/(15*P)
84          Q1=((664*H(2)-7856)*H(2)+30232)/(5*P)
85         ENDIF
86        ENDIF
87        P2=C1
88        Q2=C2
89        R2=C6
90       ELSE
91        IF(N .GE. H(7)) THEN
92         P=(7*H(2)+82)*H(2)-9
93         PP=H(2)+3
94         IF(1+P .NE. P) THEN
95          P1=((((1012*H(2)+32816)*H(2)-55496)*H(2)-104400)*H(2)+
96      1         253044)/(15*P)
97          Q1=(((664*H(2)+16600)*H(2)-24312)*H(2)+28296)/(5*P)
98         ENDIF
99        ENDIF
100        P2=C3
101        Q2=C4
102        R2=C7
103       ENDIF
104       X(1)=0
105       N0=1
106       IF(A .EQ. 0 .AND. MODE .EQ. 3) N0=2
107
108       DO 1 L = N0,N
109       IF(L .GE. H(7)) THEN
110        B=(L+H(8))*PI
111        C=C5/B**2
112        X0=B-(PP-P1*C)/(8*B*(1-Q1*C))
113       ELSE
114        IF(L .EQ. 1) THEN
115         X0=FCO(MODE)
116        ELSE
117         X0=PI1*(4*L-ICO(MODE))
118         V=1/X0**2
119         X0=-X0**E2*(1+V*(P2+Q2*V))
120        ENDIF
121        U=H(13)*X0
122        Y=E2*SQRT(ABS(U)**3)
123
124        IF(Y .LT. 1) THEN
125         G0=(3*Y)**E1
126         G1=G0**2
127         G0=G0*(1-G1*(210+G1*(G1+G1-27))/1575)
128        ELSE
129         G0=1/(Y+PIH)
130         G1=G0**2
131         G0=PIH-G0*(1+G1*(2310+G1*(3003+G1*(4818+G1*
132      1          (8591+G1*16328))))/3465)
133        ENDIF
134        G2=Y+G0
135        G1=G2**2
136        G=(G0-ATAN(G2))/G1
137        W=1/COS(G0-(1+G1)*(1+G/G2)*G)
138        G=1/(1-W**2)
139        C=SQRT(U*G)
140        X0=W*(A-P2*C*(1/U-C*(G+G-R2))/(A*U))
141       ENDIF
142
143       DO 2 J = 1,5
144       XX=X0**2
145       AX=H(1)-XX
146
147       B=X0+X0
148       E=(H(9)*X0)**2
149       P=1
150       Q=-X0
151       S=1+XX
152       R=S
153       DO 3 K = 2,500
154       IF(R*(K**2) .GE. E) GO TO 9
155       RK=Z1/(K+1)
156       D=((K-1)+H(12)/K)/S
157       P=RK*((K+K)-P*D)
158       Q=RK*(Q*D-B)
159       S=P**2+Q**2
160       R=R*S
161     3 CONTINUE
162     9 G=1/S
163       P=G*P
164       F=P
165       Q=-G*Q
166       G=Q
167       DO 4 M = K+1,2,-1
168       R=M*(2-P)-2
169       S=B+M*Q
170       D=((M-2)+H(12)/(M-1))/(R**2+S**2)
171       P=D*R
172       Q=D*S
173       E=F+1
174       F=P*E-G*Q
175       G=Q*E+P*G
176     4 CONTINUE
177       F=F+1
178       D=1/(F**2+G**2)
179       PA=F*D
180       QA=-G*D
181       D=H(11)-P
182       Q=Q+X0
183       RX=1/X0
184       PA1=RX*(PA*Q-QA*D)
185       QA1=RX*(QA*Q+PA*D)
186       XR=RX+RX
187       B=(1+A0)*XR
188       DO 5 M = 1,NA
189       P0=PA-QA1*B
190       Q0=QA+PA1*B
191       PA=PA1
192       PA1=P0
193       QA=QA1
194       QA1=Q0
195       B=B+XR
196     5 CONTINUE
197
198       SX=SIN(X0-H(10))
199       CX=COS(X0-H(10))
200       GO TO (11,12,13,14), MODE
201    11 R0=(PA*CX-QA*SX)/(PA1*SX+QA1*CX)
202       V=R0/(H(5)*X0)
203       W=R0*(1+V*(1-4*AX))/(1+V*(4*(XX-H(2))-H(6)))
204       GO TO 10
205    12 R0=(PA*SX+QA*CX)/(QA1*SX-PA1*CX)
206       V=R0/(H(5)*X0)
207       W=R0*(1+V*(1-4*AX))/(1+V*(4*(XX-H(2))-H(6)))
208       GO TO 10
209    13 R0=A/X0-(PA1*SX+QA1*CX)/(PA*CX-QA*SX)
210       GO TO 15
211    14 R0=A/X0-(QA1*SX-PA1*CX)/(PA*SX+QA*CX)
212    15 V=E2*R0*X0*AX/(H(1)+XX)
213       W=C5/AX**3
214       W=-XX*R0*(1+V*(1-W*(H(3)-XX*(40*H(2)+48*XX))))/
215      1     (AX*(1+V*(1+W*(H(4)+XX*(64*H(2)+96*XX)))))
216    10 X0=X0+W
217       IF(ABS(W) .LE. REL*X0) GO TO 6
218     2 CONTINUE
219     6 X(L)=X0
220     1 CONTINUE
221       RETURN
222   101 FORMAT('NEGATIVE INDEX A = ',1P,D15.8)
223       END