src/numerical/slatec/fortran/dbsynu.f

   1 *DECK DBSYNU
   2       SUBROUTINE DBSYNU (X, FNU, N, Y)
   3 C***BEGIN PROLOGUE  DBSYNU
   4 C***SUBSIDIARY
   5 C***PURPOSE  Subsidiary to DBESY
   6 C***LIBRARY   SLATEC
   7 C***TYPE      DOUBLE PRECISION (BESYNU-S, DBSYNU-D)
   8 C***AUTHOR  Amos, D. E., (SNLA)
   9 C***DESCRIPTION
  10 C
  11 C     Abstract  **** A DOUBLE PRECISION routine ****
  12 C         DBSYNU computes N member sequences of Y Bessel functions
  13 C         Y/SUB(FNU+I-1)/(X), I=1,N for non-negative orders FNU and
  14 C         positive X. Equations of the references are implemented on
  15 C         small orders DNU for Y/SUB(DNU)/(X) and Y/SUB(DNU+1)/(X).
  16 C         Forward recursion with the three term recursion relation
  17 C         generates higher orders FNU+I-1, I=1,...,N.
  18 C
  19 C         To start the recursion FNU is normalized to the interval
  20 C         -0.5.LE.DNU.LT.0.5. A special form of the power series is
  21 C         implemented on 0.LT.X.LE.X1 while the Miller algorithm for the
  22 C         K Bessel function in terms of the confluent hypergeometric
  23 C         function U(FNU+0.5,2*FNU+1,I*X) is implemented on X1.LT.X.LE.X
  24 C         Here I is the complex number SQRT(-1.).
  25 C         For X.GT.X2, the asymptotic expansion for large X is used.
  26 C         When FNU is a half odd integer, a special formula for
  27 C         DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion.
  28 C
  29 C         The maximum number of significant digits obtainable
  30 C         is the smaller of 14 and the number of digits carried in
  31 C         DOUBLE PRECISION arithmetic.
  32 C
  33 C         DBSYNU assumes that a significant digit SINH function is
  34 C         available.
  35 C
  36 C     Description of Arguments
  37 C
  38 C         INPUT
  39 C           X      - X.GT.0.0D0
  40 C           FNU    - Order of initial Y function, FNU.GE.0.0D0
  41 C           N      - Number of members of the sequence, N.GE.1
  42 C
  43 C         OUTPUT
  44 C           Y      - A vector whose first N components contain values
  45 C                    for the sequence Y(I)=Y/SUB(FNU+I-1), I=1,N.
  46 C
  47 C     Error Conditions
  48 C         Improper input arguments - a fatal error
  49 C         Overflow - a fatal error
  50 C
  51 C***SEE ALSO  DBESY
  52 C***REFERENCES  N. M. Temme, On the numerical evaluation of the ordinary
  53 C                 Bessel function of the second kind, Journal of
  54 C                 Computational Physics 21, (1976), pp. 343-350.
  55 C               N. M. Temme, On the numerical evaluation of the modified
  56 C                 Bessel function of the third kind, Journal of
  57 C                 Computational Physics 19, (1975), pp. 324-337.
  58 C***ROUTINES CALLED  D1MACH, DGAMMA, XERMSG
  59 C***REVISION HISTORY  (YYMMDD)
  60 C   800501  DATE WRITTEN
  61 C   890531  Changed all specific intrinsics to generic.  (WRB)
  62 C   890911  Removed unnecessary intrinsics.  (WRB)
  63 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  64 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  65 C   900326  Removed duplicate information from DESCRIPTION section.
  66 C           (WRB)
  67 C   900328  Added TYPE section.  (WRB)
  68 C   900727  Added EXTERNAL statement.  (WRB)
  69 C   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
  70 C   920501  Reformatted the REFERENCES section.  (WRB)
  71 C***END PROLOGUE  DBSYNU
  72 C
  73       INTEGER I, INU, J, K, KK, N, NN
  74       DOUBLE PRECISION A,AK,ARG,A1,A2,BK,CB,CBK,CC,CCK,CK,COEF,CPT,
  75      1 CP1, CP2, CS, CS1, CS2, CX, DNU, DNU2, ETEST, ETX, F, FC, FHS,
  76      2 FK, FKS, FLRX, FMU, FN, FNU, FX, G, G1, G2, HPI, P, PI, PT, Q,
  77      3 RB, RBK, RCK, RELB, RPT, RP1, RP2, RS, RS1, RS2, RTHPI, RX, S,
  78      4 SA, SB, SMU, SS, ST, S1, S2, TB, TM, TOL, T1, T2, X, X1, X2, Y
  79       DIMENSION A(120), RB(120), CB(120), Y(*), CC(8)
  80       DOUBLE PRECISION DGAMMA, D1MACH
  81       EXTERNAL DGAMMA
  82       SAVE X1, X2,PI, RTHPI, HPI, CC
  83       DATA X1, X2 / 3.0D0, 20.0D0 /
  84       DATA PI,RTHPI        / 3.14159265358979D+00, 7.97884560802865D-01/
  85       DATA HPI             / 1.57079632679490D+00/
  86       DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)
  87      1                     / 5.77215664901533D-01,-4.20026350340952D-02,
  88      2-4.21977345555443D-02, 7.21894324666300D-03,-2.15241674114900D-04,
  89      3-2.01348547807000D-05, 1.13302723200000D-06, 6.11609500000000D-09/
  90 C***FIRST EXECUTABLE STATEMENT  DBSYNU
  91       AK = D1MACH(3)
  92       TOL = MAX(AK,1.0D-15)
  93       IF (X.LE.0.0D0) GO TO 270
  94       IF (FNU.LT.0.0D0) GO TO 280
  95       IF (N.LT.1) GO TO 290
  96       RX = 2.0D0/X
  97       INU = INT(FNU+0.5D0)
  98       DNU = FNU - INU
  99       IF (ABS(DNU).EQ.0.5D0) GO TO 260
 100       DNU2 = 0.0D0
 101       IF (ABS(DNU).LT.TOL) GO TO 10
 102       DNU2 = DNU*DNU
 103    10 CONTINUE
 104       IF (X.GT.X1) GO TO 120
 105 C
 106 C     SERIES FOR X.LE.X1
 107 C
 108       A1 = 1.0D0 - DNU
 109       A2 = 1.0D0 + DNU
 110       T1 = 1.0D0/DGAMMA(A1)
 111       T2 = 1.0D0/DGAMMA(A2)
 112       IF (ABS(DNU).GT.0.1D0) GO TO 40
 113 C     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
 114       S = CC(1)
 115       AK = 1.0D0
 116       DO 20 K=2,8
 117         AK = AK*DNU2
 118         TM = CC(K)*AK
 119         S = S + TM
 120         IF (ABS(TM).LT.TOL) GO TO 30
 121    20 CONTINUE
 122    30 G1 = -(S+S)
 123       GO TO 50
 124    40 CONTINUE
 125       G1 = (T1-T2)/DNU
 126    50 CONTINUE
 127       G2 = T1 + T2
 128       SMU = 1.0D0
 129       FC = 1.0D0/PI
 130       FLRX = LOG(RX)
 131       FMU = DNU*FLRX
 132       TM = 0.0D0
 133       IF (DNU.EQ.0.0D0) GO TO 60
 134       TM = SIN(DNU*HPI)/DNU
 135       TM = (DNU+DNU)*TM*TM
 136       FC = DNU/SIN(DNU*PI)
 137       IF (FMU.NE.0.0D0) SMU = SINH(FMU)/FMU
 138    60 CONTINUE
 139       F = FC*(G1*COSH(FMU)+G2*FLRX*SMU)
 140       FX = EXP(FMU)
 141       P = FC*T1*FX
 142       Q = FC*T2/FX
 143       G = F + TM*Q
 144       AK = 1.0D0
 145       CK = 1.0D0
 146       BK = 1.0D0
 147       S1 = G
 148       S2 = P
 149       IF (INU.GT.0 .OR. N.GT.1) GO TO 90
 150       IF (X.LT.TOL) GO TO 80
 151       CX = X*X*0.25D0
 152    70 CONTINUE
 153       F = (AK*F+P+Q)/(BK-DNU2)
 154       P = P/(AK-DNU)
 155       Q = Q/(AK+DNU)
 156       G = F + TM*Q
 157       CK = -CK*CX/AK
 158       T1 = CK*G
 159       S1 = S1 + T1
 160       BK = BK + AK + AK + 1.0D0
 161       AK = AK + 1.0D0
 162       S = ABS(T1)/(1.0D0+ABS(S1))
 163       IF (S.GT.TOL) GO TO 70
 164    80 CONTINUE
 165       Y(1) = -S1
 166       RETURN
 167    90 CONTINUE
 168       IF (X.LT.TOL) GO TO 110
 169       CX = X*X*0.25D0
 170   100 CONTINUE
 171       F = (AK*F+P+Q)/(BK-DNU2)
 172       P = P/(AK-DNU)
 173       Q = Q/(AK+DNU)
 174       G = F + TM*Q
 175       CK = -CK*CX/AK
 176       T1 = CK*G
 177       S1 = S1 + T1
 178       T2 = CK*(P-AK*G)
 179       S2 = S2 + T2
 180       BK = BK + AK + AK + 1.0D0
 181       AK = AK + 1.0D0
 182       S = ABS(T1)/(1.0D0+ABS(S1)) + ABS(T2)/(1.0D0+ABS(S2))
 183       IF (S.GT.TOL) GO TO 100
 184   110 CONTINUE
 185       S2 = -S2*RX
 186       S1 = -S1
 187       GO TO 160
 188   120 CONTINUE
 189       COEF = RTHPI/SQRT(X)
 190       IF (X.GT.X2) GO TO 210
 191 C
 192 C     MILLER ALGORITHM FOR X1.LT.X.LE.X2
 193 C
 194       ETEST = COS(PI*DNU)/(PI*X*TOL)
 195       FKS = 1.0D0
 196       FHS = 0.25D0
 197       FK = 0.0D0
 198       RCK = 2.0D0
 199       CCK = X + X
 200       RP1 = 0.0D0
 201       CP1 = 0.0D0
 202       RP2 = 1.0D0
 203       CP2 = 0.0D0
 204       K = 0
 205   130 CONTINUE
 206       K = K + 1
 207       FK = FK + 1.0D0
 208       AK = (FHS-DNU2)/(FKS+FK)
 209       PT = FK + 1.0D0
 210       RBK = RCK/PT
 211       CBK = CCK/PT
 212       RPT = RP2
 213       CPT = CP2
 214       RP2 = RBK*RPT - CBK*CPT - AK*RP1
 215       CP2 = CBK*RPT + RBK*CPT - AK*CP1
 216       RP1 = RPT
 217       CP1 = CPT
 218       RB(K) = RBK
 219       CB(K) = CBK
 220       A(K) = AK
 221       RCK = RCK + 2.0D0
 222       FKS = FKS + FK + FK + 1.0D0
 223       FHS = FHS + FK + FK
 224       PT = MAX(ABS(RP1),ABS(CP1))
 225       FC = (RP1/PT)**2 + (CP1/PT)**2
 226       PT = PT*SQRT(FC)*FK
 227       IF (ETEST.GT.PT) GO TO 130
 228       KK = K
 229       RS = 1.0D0
 230       CS = 0.0D0
 231       RP1 = 0.0D0
 232       CP1 = 0.0D0
 233       RP2 = 1.0D0
 234       CP2 = 0.0D0
 235       DO 140 I=1,K
 236         RPT = RP2
 237         CPT = CP2
 238         RP2 = (RB(KK)*RPT-CB(KK)*CPT-RP1)/A(KK)
 239         CP2 = (CB(KK)*RPT+RB(KK)*CPT-CP1)/A(KK)
 240         RP1 = RPT
 241         CP1 = CPT
 242         RS = RS + RP2
 243         CS = CS + CP2
 244         KK = KK - 1
 245   140 CONTINUE
 246       PT = MAX(ABS(RS),ABS(CS))
 247       FC = (RS/PT)**2 + (CS/PT)**2
 248       PT = PT*SQRT(FC)
 249       RS1 = (RP2*(RS/PT)+CP2*(CS/PT))/PT
 250       CS1 = (CP2*(RS/PT)-RP2*(CS/PT))/PT
 251       FC = HPI*(DNU-0.5D0) - X
 252       P = COS(FC)
 253       Q = SIN(FC)
 254       S1 = (CS1*Q-RS1*P)*COEF
 255       IF (INU.GT.0 .OR. N.GT.1) GO TO 150
 256       Y(1) = S1
 257       RETURN
 258   150 CONTINUE
 259       PT = MAX(ABS(RP2),ABS(CP2))
 260       FC = (RP2/PT)**2 + (CP2/PT)**2
 261       PT = PT*SQRT(FC)
 262       RPT = DNU + 0.5D0 - (RP1*(RP2/PT)+CP1*(CP2/PT))/PT
 263       CPT = X - (CP1*(RP2/PT)-RP1*(CP2/PT))/PT
 264       CS2 = CS1*CPT - RS1*RPT
 265       RS2 = RPT*CS1 + RS1*CPT
 266       S2 = (RS2*Q+CS2*P)*COEF/X
 267 C
 268 C     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
 269 C
 270   160 CONTINUE
 271       CK = (DNU+DNU+2.0D0)/X
 272       IF (N.EQ.1) INU = INU - 1
 273       IF (INU.GT.0) GO TO 170
 274       IF (N.GT.1) GO TO 190
 275       S1 = S2
 276       GO TO 190
 277   170 CONTINUE
 278       DO 180 I=1,INU
 279         ST = S2
 280         S2 = CK*S2 - S1
 281         S1 = ST
 282         CK = CK + RX
 283   180 CONTINUE
 284       IF (N.EQ.1) S1 = S2
 285   190 CONTINUE
 286       Y(1) = S1
 287       IF (N.EQ.1) RETURN
 288       Y(2) = S2
 289       IF (N.EQ.2) RETURN
 290       DO 200 I=3,N
 291         Y(I) = CK*Y(I-1) - Y(I-2)
 292         CK = CK + RX
 293   200 CONTINUE
 294       RETURN
 295 C
 296 C     ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2
 297 C
 298   210 CONTINUE
 299       NN = 2
 300       IF (INU.EQ.0 .AND. N.EQ.1) NN = 1
 301       DNU2 = DNU + DNU
 302       FMU = 0.0D0
 303       IF (ABS(DNU2).LT.TOL) GO TO 220
 304       FMU = DNU2*DNU2
 305   220 CONTINUE
 306       ARG = X - HPI*(DNU+0.5D0)
 307       SA = SIN(ARG)
 308       SB = COS(ARG)
 309       ETX = 8.0D0*X
 310       DO 250 K=1,NN
 311         S1 = S2
 312         T2 = (FMU-1.0D0)/ETX
 313         SS = T2
 314         RELB = TOL*ABS(T2)
 315         T1 = ETX
 316         S = 1.0D0
 317         FN = 1.0D0
 318         AK = 0.0D0
 319         DO 230 J=1,13
 320           T1 = T1 + ETX
 321           AK = AK + 8.0D0
 322           FN = FN + AK
 323           T2 = -T2*(FMU-FN)/T1
 324           S = S + T2
 325           T1 = T1 + ETX
 326           AK = AK + 8.0D0
 327           FN = FN + AK
 328           T2 = T2*(FMU-FN)/T1
 329           SS = SS + T2
 330           IF (ABS(T2).LE.RELB) GO TO 240
 331   230   CONTINUE
 332   240   S2 = COEF*(S*SA+SS*SB)
 333         FMU = FMU + 8.0D0*DNU + 4.0D0
 334         TB = SA
 335         SA = -SB
 336         SB = TB
 337   250 CONTINUE
 338       IF (NN.GT.1) GO TO 160
 339       S1 = S2
 340       GO TO 190
 341 C
 342 C     FNU=HALF ODD INTEGER CASE
 343 C
 344   260 CONTINUE
 345       COEF = RTHPI/SQRT(X)
 346       S1 = COEF*SIN(X)
 347       S2 = -COEF*COS(X)
 348       GO TO 160
 349 C
 350 C
 351   270 CALL XERMSG ('SLATEC', 'DBSYNU', 'X NOT GREATER THAN ZERO', 2, 1)
 352       RETURN
 353   280 CALL XERMSG ('SLATEC', 'DBSYNU', 'FNU NOT ZERO OR POSITIVE', 2,
 354      +   1)
 355       RETURN
 356   290 CALL XERMSG ('SLATEC', 'DBSYNU', 'N NOT GREATER THAN 0', 2, 1)
 357       RETURN
 358       END