GNU Octave: libcruft/amos/cbknu.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CBKNU(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
00002 C***BEGIN PROLOGUE  CBKNU
00003 C***REFER TO  CBESI,CBESK,CAIRY,CBESH
00004 C
00005 C     CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE
00006 C
00007 C***ROUTINES CALLED  CKSCL,CSHCH,GAMLN,I1MACH,R1MACH,CUCHK
00008 C***END PROLOGUE  CBKNU
00009 C
00010       COMPLEX CCH, CK, COEF, CONE, CRSC, CS, CSCL, CSH, CSR, CSS, CTWO,
00011      * CZ, CZERO, F, FMU, P, PT, P1, P2, Q, RZ, SMU, ST, S1, S2, Y, Z,
00012      * ZD, CELM, CY
00013       REAL AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, CC, DNU,
00014      * DNU2, ELIM, ETEST, FC, FHS, FK, FKS, FNU, FPI, G1, G2, HPI, PI,
00015      * P2I, P2M, P2R, RK, RTHPI, R1, S, SPI, TM, TOL, TTH, T1, T2, XX,
00016      * YY, GAMLN, R1MACH, HELIM, ELM, XD, YD, ALAS, AS
00017       INTEGER I, IDUM, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N,
00018      * NZ, I1MACH, NW, J, IC, INUB
00019       DIMENSION BRY(3), CC(8), CSS(3), CSR(3), Y(N), CY(2)
00020 C
00021       DATA KMAX / 30 /
00022       DATA R1 / 2.0E0 /
00023       DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/
00024 C
00025       DATA PI, RTHPI, SPI ,HPI, FPI, TTH /
00026      1     3.14159265358979324E0,       1.25331413731550025E0,
00027      2     1.90985931710274403E0,       1.57079632679489662E0,
00028      3     1.89769999331517738E0,       6.66666666666666666E-01/
00029 C
00030       DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
00031      1     5.77215664901532861E-01,    -4.20026350340952355E-02,
00032      2    -4.21977345555443367E-02,     7.21894324666309954E-03,
00033      3    -2.15241674114950973E-04,    -2.01348547807882387E-05,
00034      4     1.13302723198169588E-06,     6.11609510448141582E-09/
00035 C
00036       XX = REAL(Z)
00037       YY = AIMAG(Z)
00038       CAZ = CABS(Z)
00039       CSCL = CMPLX(1.0E0/TOL,0.0E0)
00040       CRSC = CMPLX(TOL,0.0E0)
00041       CSS(1) = CSCL
00042       CSS(2) = CONE
00043       CSS(3) = CRSC
00044       CSR(1) = CRSC
00045       CSR(2) = CONE
00046       CSR(3) = CSCL
00047       BRY(1) = 1.0E+3*R1MACH(1)/TOL
00048       BRY(2) = 1.0E0/BRY(1)
00049       BRY(3) = R1MACH(2)
00050       NZ = 0
00051       IFLAG = 0
00052       KODED = KODE
00053       RZ = CTWO/Z
00054       INU = INT(FNU+0.5E0)
00055       DNU = FNU - FLOAT(INU)
00056       IF (ABS(DNU).EQ.0.5E0) GO TO 110
00057       DNU2 = 0.0E0
00058       IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU
00059       IF (CAZ.GT.R1) GO TO 110
00060 C-----------------------------------------------------------------------
00061 C     SERIES FOR CABS(Z).LE.R1
00062 C-----------------------------------------------------------------------
00063       FC = 1.0E0
00064       SMU = CLOG(RZ)
00065       FMU = SMU*CMPLX(DNU,0.0E0)
00066       CALL CSHCH(FMU, CSH, CCH)
00067       IF (DNU.EQ.0.0E0) GO TO 10
00068       FC = DNU*PI
00069       FC = FC/SIN(FC)
00070       SMU = CSH*CMPLX(1.0E0/DNU,0.0E0)
00071    10 CONTINUE
00072       A2 = 1.0E0 + DNU
00073 C-----------------------------------------------------------------------
00074 C     GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
00075 C-----------------------------------------------------------------------
00076       T2 = EXP(-GAMLN(A2,IDUM))
00077       T1 = 1.0E0/(T2*FC)
00078       IF (ABS(DNU).GT.0.1E0) GO TO 40
00079 C-----------------------------------------------------------------------
00080 C     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
00081 C-----------------------------------------------------------------------
00082       AK = 1.0E0
00083       S = CC(1)
00084       DO 20 K=2,8
00085         AK = AK*DNU2
00086         TM = CC(K)*AK
00087         S = S + TM
00088         IF (ABS(TM).LT.TOL) GO TO 30
00089    20 CONTINUE
00090    30 G1 = -S
00091       GO TO 50
00092    40 CONTINUE
00093       G1 = (T1-T2)/(DNU+DNU)
00094    50 CONTINUE
00095       G2 = 0.5E0*(T1+T2)*FC
00096       G1 = G1*FC
00097       F = CMPLX(G1,0.0E0)*CCH + SMU*CMPLX(G2,0.0E0)
00098       PT = CEXP(FMU)
00099       P = CMPLX(0.5E0/T2,0.0E0)*PT
00100       Q = CMPLX(0.5E0/T1,0.0E0)/PT
00101       S1 = F
00102       S2 = P
00103       AK = 1.0E0
00104       A1 = 1.0E0
00105       CK = CONE
00106       BK = 1.0E0 - DNU2
00107       IF (INU.GT.0 .OR. N.GT.1) GO TO 80
00108 C-----------------------------------------------------------------------
00109 C     GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
00110 C-----------------------------------------------------------------------
00111       IF (CAZ.LT.TOL) GO TO 70
00112       CZ = Z*Z*CMPLX(0.25E0,0.0E0)
00113       T1 = 0.25E0*CAZ*CAZ
00114    60 CONTINUE
00115       F = (F*CMPLX(AK,0.0E0)+P+Q)*CMPLX(1.0E0/BK,0.0E0)
00116       P = P*CMPLX(1.0E0/(AK-DNU),0.0E0)
00117       Q = Q*CMPLX(1.0E0/(AK+DNU),0.0E0)
00118       RK = 1.0E0/AK
00119       CK = CK*CZ*CMPLX(RK,0.0)
00120       S1 = S1 + CK*F
00121       A1 = A1*T1*RK
00122       BK = BK + AK + AK + 1.0E0
00123       AK = AK + 1.0E0
00124       IF (A1.GT.TOL) GO TO 60
00125    70 CONTINUE
00126       Y(1) = S1
00127       IF (KODED.EQ.1) RETURN
00128       Y(1) = S1*CEXP(Z)
00129       RETURN
00130 C-----------------------------------------------------------------------
00131 C     GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
00132 C-----------------------------------------------------------------------
00133    80 CONTINUE
00134       IF (CAZ.LT.TOL) GO TO 100
00135       CZ = Z*Z*CMPLX(0.25E0,0.0E0)
00136       T1 = 0.25E0*CAZ*CAZ
00137    90 CONTINUE
00138       F = (F*CMPLX(AK,0.0E0)+P+Q)*CMPLX(1.0E0/BK,0.0E0)
00139       P = P*CMPLX(1.0E0/(AK-DNU),0.0E0)
00140       Q = Q*CMPLX(1.0E0/(AK+DNU),0.0E0)
00141       RK = 1.0E0/AK
00142       CK = CK*CZ*CMPLX(RK,0.0E0)
00143       S1 = S1 + CK*F
00144       S2 = S2 + CK*(P-F*CMPLX(AK,0.0E0))
00145       A1 = A1*T1*RK
00146       BK = BK + AK + AK + 1.0E0
00147       AK = AK + 1.0E0
00148       IF (A1.GT.TOL) GO TO 90
00149   100 CONTINUE
00150       KFLAG = 2
00151       BK = REAL(SMU)
00152       A1 = FNU + 1.0E0
00153       AK = A1*ABS(BK)
00154       IF (AK.GT.ALIM) KFLAG = 3
00155       P2 = S2*CSS(KFLAG)
00156       S2 = P2*RZ
00157       S1 = S1*CSS(KFLAG)
00158       IF (KODED.EQ.1) GO TO 210
00159       F = CEXP(Z)
00160       S1 = S1*F
00161       S2 = S2*F
00162       GO TO 210
00163 C-----------------------------------------------------------------------
00164 C     IFLAG=0 MEANS NO UNDERFLOW OCCURRED
00165 C     IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
00166 C     KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
00167 C     RECURSION
00168 C-----------------------------------------------------------------------
00169   110 CONTINUE
00170       COEF = CMPLX(RTHPI,0.0E0)/CSQRT(Z)
00171       KFLAG = 2
00172       IF (KODED.EQ.2) GO TO 120
00173       IF (XX.GT.ALIM) GO TO 290
00174 C     BLANK LINE
00175       A1 = EXP(-XX)*REAL(CSS(KFLAG))
00176       PT = CMPLX(A1,0.0E0)*CMPLX(COS(YY),-SIN(YY))
00177       COEF = COEF*PT
00178   120 CONTINUE
00179       IF (ABS(DNU).EQ.0.5E0) GO TO 300
00180 C-----------------------------------------------------------------------
00181 C     MILLER ALGORITHM FOR CABS(Z).GT.R1
00182 C-----------------------------------------------------------------------
00183       AK = COS(PI*DNU)
00184       AK = ABS(AK)
00185       IF (AK.EQ.0.0E0) GO TO 300
00186       FHS = ABS(0.25E0-DNU2)
00187       IF (FHS.EQ.0.0E0) GO TO 300
00188 C-----------------------------------------------------------------------
00189 C     COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
00190 C     DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
00191 C     12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(11))=
00192 C     TOL WHERE B IS THE BASE OF THE ARITHMETIC.
00193 C-----------------------------------------------------------------------
00194       T1 = FLOAT(I1MACH(11)-1)*R1MACH(5)*3.321928094E0
00195       T1 = AMAX1(T1,12.0E0)
00196       T1 = AMIN1(T1,60.0E0)
00197       T2 = TTH*T1 - 6.0E0
00198       IF (XX.NE.0.0E0) GO TO 130
00199       T1 = HPI
00200       GO TO 140
00201   130 CONTINUE
00202       T1 = ATAN(YY/XX)
00203       T1 = ABS(T1)
00204   140 CONTINUE
00205       IF (T2.GT.CAZ) GO TO 170
00206 C-----------------------------------------------------------------------
00207 C     FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
00208 C-----------------------------------------------------------------------
00209       ETEST = AK/(PI*CAZ*TOL)
00210       FK = 1.0E0
00211       IF (ETEST.LT.1.0E0) GO TO 180
00212       FKS = 2.0E0
00213       RK = CAZ + CAZ + 2.0E0
00214       A1 = 0.0E0
00215       A2 = 1.0E0
00216       DO 150 I=1,KMAX
00217         AK = FHS/FKS
00218         BK = RK/(FK+1.0E0)
00219         TM = A2
00220         A2 = BK*A2 - AK*A1
00221         A1 = TM
00222         RK = RK + 2.0E0
00223         FKS = FKS + FK + FK + 2.0E0
00224         FHS = FHS + FK + FK
00225         FK = FK + 1.0E0
00226         TM = ABS(A2)*FK
00227         IF (ETEST.LT.TM) GO TO 160
00228   150 CONTINUE
00229       GO TO 310
00230   160 CONTINUE
00231       FK = FK + SPI*T1*SQRT(T2/CAZ)
00232       FHS = ABS(0.25E0-DNU2)
00233       GO TO 180
00234   170 CONTINUE
00235 C-----------------------------------------------------------------------
00236 C     COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
00237 C-----------------------------------------------------------------------
00238       A2 = SQRT(CAZ)
00239       AK = FPI*AK/(TOL*SQRT(A2))
00240       AA = 3.0E0*T1/(1.0E0+CAZ)
00241       BB = 14.7E0*T1/(28.0E0+CAZ)
00242       AK = (ALOG(AK)+CAZ*COS(AA)/(1.0E0+0.008E0*CAZ))/COS(BB)
00243       FK = 0.12125E0*AK*AK/CAZ + 1.5E0
00244   180 CONTINUE
00245       K = INT(FK)
00246 C-----------------------------------------------------------------------
00247 C     BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
00248 C-----------------------------------------------------------------------
00249       FK = FLOAT(K)
00250       FKS = FK*FK
00251       P1 = CZERO
00252       P2 = CMPLX(TOL,0.0E0)
00253       CS = P2
00254       DO 190 I=1,K
00255         A1 = FKS - FK
00256         A2 = (FKS+FK)/(A1+FHS)
00257         RK = 2.0E0/(FK+1.0E0)
00258         T1 = (FK+XX)*RK
00259         T2 = YY*RK
00260         PT = P2
00261         P2 = (P2*CMPLX(T1,T2)-P1)*CMPLX(A2,0.0E0)
00262         P1 = PT
00263         CS = CS + P2
00264         FKS = A1 - FK + 1.0E0
00265         FK = FK - 1.0E0
00266   190 CONTINUE
00267 C-----------------------------------------------------------------------
00268 C     COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
00269 C     SCALING
00270 C-----------------------------------------------------------------------
00271       TM = CABS(CS)
00272       PT = CMPLX(1.0E0/TM,0.0E0)
00273       S1 = PT*P2
00274       CS = CONJG(CS)*PT
00275       S1 = COEF*S1*CS
00276       IF (INU.GT.0 .OR. N.GT.1) GO TO 200
00277       ZD = Z
00278       IF(IFLAG.EQ.1) GO TO 270
00279       GO TO 240
00280   200 CONTINUE
00281 C-----------------------------------------------------------------------
00282 C     COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
00283 C-----------------------------------------------------------------------
00284       TM = CABS(P2)
00285       PT = CMPLX(1.0E0/TM,0.0E0)
00286       P1 = PT*P1
00287       P2 = CONJG(P2)*PT
00288       PT = P1*P2
00289       S2 = S1*(CONE+(CMPLX(DNU+0.5E0,0.0E0)-PT)/Z)
00290 C-----------------------------------------------------------------------
00291 C     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
00292 C     SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
00293 C-----------------------------------------------------------------------
00294   210 CONTINUE
00295       CK = CMPLX(DNU+1.0E0,0.0E0)*RZ
00296       IF (N.EQ.1) INU = INU - 1
00297       IF (INU.GT.0) GO TO 220
00298       IF (N.EQ.1) S1=S2
00299       ZD = Z
00300       IF(IFLAG.EQ.1) GO TO 270
00301       GO TO 240
00302   220 CONTINUE
00303       INUB = 1
00304       IF (IFLAG.EQ.1) GO TO 261
00305   225 CONTINUE
00306       P1 = CSR(KFLAG)
00307       ASCLE = BRY(KFLAG)
00308       DO 230 I=INUB,INU
00309         ST = S2
00310         S2 = CK*S2 + S1
00311         S1 = ST
00312         CK = CK + RZ
00313         IF (KFLAG.GE.3) GO TO 230
00314         P2 = S2*P1
00315         P2R = REAL(P2)
00316         P2I = AIMAG(P2)
00317         P2R = ABS(P2R)
00318         P2I = ABS(P2I)
00319         P2M = AMAX1(P2R,P2I)
00320         IF (P2M.LE.ASCLE) GO TO 230
00321         KFLAG = KFLAG + 1
00322         ASCLE = BRY(KFLAG)
00323         S1 = S1*P1
00324         S2 = P2
00325         S1 = S1*CSS(KFLAG)
00326         S2 = S2*CSS(KFLAG)
00327         P1 = CSR(KFLAG)
00328   230 CONTINUE
00329       IF (N.EQ.1) S1 = S2
00330   240 CONTINUE
00331       Y(1) = S1*CSR(KFLAG)
00332       IF (N.EQ.1) RETURN
00333       Y(2) = S2*CSR(KFLAG)
00334       IF (N.EQ.2) RETURN
00335       KK = 2
00336   250 CONTINUE
00337       KK = KK + 1
00338       IF (KK.GT.N) RETURN
00339       P1 = CSR(KFLAG)
00340       ASCLE = BRY(KFLAG)
00341       DO 260 I=KK,N
00342         P2 = S2
00343         S2 = CK*S2 + S1
00344         S1 = P2
00345         CK = CK + RZ
00346         P2 = S2*P1
00347         Y(I) = P2
00348         IF (KFLAG.GE.3) GO TO 260
00349         P2R = REAL(P2)
00350         P2I = AIMAG(P2)
00351         P2R = ABS(P2R)
00352         P2I = ABS(P2I)
00353         P2M = AMAX1(P2R,P2I)
00354         IF (P2M.LE.ASCLE) GO TO 260
00355         KFLAG = KFLAG + 1
00356         ASCLE = BRY(KFLAG)
00357         S1 = S1*P1
00358         S2 = P2
00359         S1 = S1*CSS(KFLAG)
00360         S2 = S2*CSS(KFLAG)
00361         P1 = CSR(KFLAG)
00362   260 CONTINUE
00363       RETURN
00364 C-----------------------------------------------------------------------
00365 C     IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
00366 C-----------------------------------------------------------------------
00367   261 CONTINUE
00368       HELIM = 0.5E0*ELIM
00369       ELM = EXP(-ELIM)
00370       CELM = CMPLX(ELM,0.0)
00371       ASCLE = BRY(1)
00372       ZD = Z
00373       XD = XX
00374       YD = YY
00375       IC = -1
00376       J = 2
00377       DO 262 I=1,INU
00378         ST = S2
00379         S2 = CK*S2+S1
00380         S1 = ST
00381         CK = CK+RZ
00382         AS = CABS(S2)
00383         ALAS = ALOG(AS)
00384         P2R = -XD+ALAS
00385         IF(P2R.LT.(-ELIM)) GO TO 263
00386         P2 = -ZD+CLOG(S2)
00387         P2R = REAL(P2)
00388         P2I = AIMAG(P2)
00389         P2M = EXP(P2R)/TOL
00390         P1 = CMPLX(P2M,0.0E0)*CMPLX(COS(P2I),SIN(P2I))
00391         CALL CUCHK(P1,NW,ASCLE,TOL)
00392         IF(NW.NE.0) GO TO 263
00393         J=3-J
00394         CY(J) = P1
00395         IF(IC.EQ.(I-1)) GO TO 264
00396         IC = I
00397         GO TO 262
00398   263   CONTINUE
00399         IF(ALAS.LT.HELIM) GO TO 262
00400         XD = XD-ELIM
00401         S1 = S1*CELM
00402         S2 = S2*CELM
00403         ZD = CMPLX(XD,YD)
00404   262 CONTINUE
00405       IF(N.EQ.1) S1 = S2
00406       GO TO 270
00407   264 CONTINUE
00408       KFLAG = 1
00409       INUB = I+1
00410       S2 = CY(J)
00411       J = 3 - J
00412       S1 = CY(J)
00413       IF(INUB.LE.INU) GO TO 225
00414       IF(N.EQ.1) S1 = S2
00415       GO TO 240
00416   270 CONTINUE
00417       Y(1) = S1
00418       IF (N.EQ.1) GO TO 280
00419       Y(2) = S2
00420   280 CONTINUE
00421       ASCLE = BRY(1)
00422       CALL CKSCL(ZD, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
00423       INU = N - NZ
00424       IF (INU.LE.0) RETURN
00425       KK = NZ + 1
00426       S1 = Y(KK)
00427       Y(KK) = S1*CSR(1)
00428       IF (INU.EQ.1) RETURN
00429       KK = NZ + 2
00430       S2 = Y(KK)
00431       Y(KK) = S2*CSR(1)
00432       IF (INU.EQ.2) RETURN
00433       T2 = FNU + FLOAT(KK-1)
00434       CK = CMPLX(T2,0.0E0)*RZ
00435       KFLAG = 1
00436       GO TO 250
00437   290 CONTINUE
00438 C-----------------------------------------------------------------------
00439 C     SCALE BY EXP(Z), IFLAG = 1 CASES
00440 C-----------------------------------------------------------------------
00441       KODED = 2
00442       IFLAG = 1
00443       KFLAG = 2
00444       GO TO 120
00445 C-----------------------------------------------------------------------
00446 C     FNU=HALF ODD INTEGER CASE, DNU=-0.5
00447 C-----------------------------------------------------------------------
00448   300 CONTINUE
00449       S1 = COEF
00450       S2 = COEF
00451       GO TO 210
00452   310 CONTINUE
00453       NZ=-2
00454       RETURN
00455       END