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