SUBROUTINE DVK0 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE DVK0 C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order zero (K0) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE DOUBLE PRECISION (VK0-S, DVK0-D) C***KEYWORDS BESSEL FUNCTION,THIRD KIND,HYPERBOLIC BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION, ORDER ZERO, VECTORIZED C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C DVK0 computes the modified (hyperbolic) Bessel function of the C third kind of order zero (K0) for real arguments using uniform C approximation by Chebyshev polynomials. C C C P A R A M E T E R S C C M (Input) Integer C The number of arguments at which the function is to be C evaluated. C C X (Input) Double precision array of length M C The arguments at which the function is to be evaluated are C stored in X(1) to X(M) in any order. C C F (Output) Double precision array of length M C F(i) contains the value of the function at X(i), i=1,..,M. C C WORK (Work) Double precision vector of length 7*M C C IWORK (Work) Integer vector of length M C C INFO (Output) Integer C Indicates status of computed result. The following table C lists possible values and their meanings. If OK=Yes then C all F(i) have been set, otherwise none have been set. C C INFO OK Description C ------------------------------------------------------------ C -1 Yes Warning: Some X(i) so big K0 underflows. C The corresponding F(i) are set to zero. C 0 Yes Successfull execution. C 1 No Error: M .LE. 0 C 2 No Error: Some X(i) is zero or negative. C The index of the first offending argument is C returned in IWORK(1). C C C ********************************************************************* C This routine is a modification of the function DBESK0 developed by C W. Fullerton of LANL. C ********************************************************************* C C***REFERENCES Ronald F. Boisvert and Bonita V. Saunders, Portable C Vectorized Software for Bessel Function Evaluation, C ACM Transactions on Mathematical Software 18 (1992), C pp 456-469. C***ROUTINES CALLED DWK0 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE DVK0 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, IWORK, M DOUBLE PRECISION F, X, WORK C DIMENSION X(M), F(M), WORK(7*M), IWORK(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER IWB0, IWB1, IWB2, IWTC, IWXC, IWYC, IWZC, JIN C C***FIRST EXECUTABLE STATEMENT DVK0 C C ... PARTITION WORK ARRAYS C IWTC = 1 IWXC = IWTC + M IWYC = IWXC + M IWZC = IWYC + M IWB0 = IWZC + M IWB1 = IWB0 + M IWB2 = IWB1 + M C Total = IB2 + M C JIN = 1 C Total = JIN + M C C ... DWK0 DOES ALL THE WORK C CALL DWK0(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE DWK0 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1, + B2, INFO) C***BEGIN PROLOGUE DWK0 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order zero (K0) for a vector of arguments C***LIBRARY VFNLIB C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***ROUTINES CALLED IDWCS, D1MACH, DWNLE, DWGTHR, DWGTLE, DWGT, DWLE, C DWSCTR, DWCS C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE DWK0 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, INDX, M DOUBLE PRECISION B0, B1, B2, F, X, TCMP, XCMP, YCMP, ZCMP C DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), TCMP(M), + XCMP(M), YCMP(M), ZCMP(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER LAK0, LAK02, LBI0, LBK0 PARAMETER ( LAK0=38, LAK02=33, LBI0=18, LBK0=16 ) C INTEGER I, IDWCS, J, KEY, N, NTI0, NTK0, NTAK0, NTAK02, NTOT DOUBLE PRECISION AK0CS, AK02CS, BI0CS, BK0CS, EPMACH, EPS, D1MACH, + XMAX, XSML C DIMENSION AK0CS(LAK0), AK02CS(LAK02), BI0CS(LBI0), BK0CS(LBK0) C SAVE BK0CS, AK0CS, AK02CS, BI0CS, N, NTAK0, NTAK02, NTI0, NTK0, + XMAX, XSML C C---------------------------------------------------------------------- C C Series for BK0 on the interval 0. to 4.00000D+00 C with weighted error 3.08E-33 C log weighted error 32.51 C significant figures required 32.05 C decimal places required 33.11 C DATA BK0 CS( 1) / -.3532739323 3902768720 1140060063 153 D-1 / DATA BK0 CS( 2) / +.3442898999 2462848688 6344927529 213 D+0 / DATA BK0 CS( 3) / +.3597993651 5361501626 5721303687 231 D-1 / DATA BK0 CS( 4) / +.1264615411 4469259233 8479508673 447 D-2 / DATA BK0 CS( 5) / +.2286212103 1194517860 8269830297 585 D-4 / DATA BK0 CS( 6) / +.2534791079 0261494573 0790013428 354 D-6 / DATA BK0 CS( 7) / +.1904516377 2202088589 7214059381 366 D-8 / DATA BK0 CS( 8) / +.1034969525 7633624585 1008317853 089 D-10 / DATA BK0 CS( 9) / +.4259816142 7910825765 2445327170 133 D-13 / DATA BK0 CS( 10) / +.1374465435 8807508969 4238325440 000 D-15 / DATA BK0 CS( 11) / +.3570896528 5083735909 9688597333 333 D-18 / DATA BK0 CS( 12) / +.7631643660 1164373766 7498666666 666 D-21 / DATA BK0 CS( 13) / +.1365424988 4407818590 8053333333 333 D-23 / DATA BK0 CS( 14) / +.2075275266 9066680831 9999999999 999 D-26 / DATA BK0 CS( 15) / +.2712814218 0729856000 0000000000 000 D-29 / DATA BK0 CS( 16) / +.3082593887 9146666666 6666666666 666 D-32 / C C------------------------------------------------------------------- C C Series for AK0 on the interval 1.25000D-01 to 5.00000D-01 C with weighted error 2.85E-32 C log weighted error 31.54 C significant figures required 30.19 C decimal places required 32.33 C DATA AK0 CS( 1) / -.7643947903 3279414240 8297827008 8 D-1 / DATA AK0 CS( 2) / -.2235652605 6998190520 2309555079 1 D-1 / DATA AK0 CS( 3) / +.7734181154 6938582353 0061817404 7 D-3 / DATA AK0 CS( 4) / -.4281006688 8860994644 5214643541 6 D-4 / DATA AK0 CS( 5) / +.3081700173 8629747436 5001482666 0 D-5 / DATA AK0 CS( 6) / -.2639367222 0096649740 6744889272 3 D-6 / DATA AK0 CS( 7) / +.2563713036 4034692062 9408826574 2 D-7 / DATA AK0 CS( 8) / -.2742705549 9002012638 5721191524 4 D-8 / DATA AK0 CS( 9) / +.3169429658 0974995920 8083287340 3 D-9 / DATA AK0 CS( 10) / -.3902353286 9621841416 0106571796 2 D-10 / DATA AK0 CS( 11) / +.5068040698 1885754020 5009212728 6 D-11 / DATA AK0 CS( 12) / -.6889574741 0078706795 4171355798 4 D-12 / DATA AK0 CS( 13) / +.9744978497 8259176913 8820133683 1 D-13 / DATA AK0 CS( 14) / -.1427332841 8845485053 8985534012 2 D-13 / DATA AK0 CS( 15) / +.2156412571 0214630395 5806297652 7 D-14 / DATA AK0 CS( 16) / -.3349654255 1495627721 8878205853 0 D-15 / DATA AK0 CS( 17) / +.5335260216 9529116921 4528039260 1 D-16 / DATA AK0 CS( 18) / -.8693669980 8907538076 3962237883 7 D-17 / DATA AK0 CS( 19) / +.1446404347 8622122278 8776344234 6 D-17 / DATA AK0 CS( 20) / -.2452889825 5001296824 0467875157 3 D-18 / DATA AK0 CS( 21) / +.4233754526 2321715728 2170634240 0 D-19 / DATA AK0 CS( 22) / -.7427946526 4544641956 9534129493 3 D-20 / DATA AK0 CS( 23) / +.1323150529 3926668662 7796746240 0 D-20 / DATA AK0 CS( 24) / -.2390587164 7396494513 3598146559 9 D-21 / DATA AK0 CS( 25) / +.4376827585 9232261401 6571255466 6 D-22 / DATA AK0 CS( 26) / -.8113700607 3451180593 3901141333 3 D-23 / DATA AK0 CS( 27) / +.1521819913 8321729583 1037815466 6 D-23 / DATA AK0 CS( 28) / -.2886041941 4833977702 3595861333 3 D-24 / DATA AK0 CS( 29) / +.5530620667 0547179799 9261013333 3 D-25 / DATA AK0 CS( 30) / -.1070377329 2498987285 9163306666 6 D-25 / DATA AK0 CS( 31) / +.2091086893 1423843002 9632853333 3 D-26 / DATA AK0 CS( 32) / -.4121713723 6462038274 1026133333 3 D-27 / DATA AK0 CS( 33) / +.8193483971 1213076401 3568000000 0 D-28 / DATA AK0 CS( 34) / -.1642000275 4592977267 8075733333 3 D-28 / DATA AK0 CS( 35) / +.3316143281 4802271958 9034666666 6 D-29 / DATA AK0 CS( 36) / -.6746863644 1452959410 8586666666 6 D-30 / DATA AK0 CS( 37) / +.1382429146 3184246776 3541333333 3 D-30 / DATA AK0 CS( 38) / -.2851874167 3598325708 1173333333 3 D-31 / C C------------------------------------------------------------------- C C Series for AK02 on the interval 0. to 1.25000D-01 C with weighted error 2.30E-32 C log weighted error 31.64 C significant figures required 29.68 C decimal places required 32.40 C DATA AK02CS( 1) / -.1201869826 3075922398 3934621245 2 D-1 / DATA AK02CS( 2) / -.9174852691 0256953106 5256107571 3 D-2 / DATA AK02CS( 3) / +.1444550931 7750058210 4884387805 7 D-3 / DATA AK02CS( 4) / -.4013614175 4357097286 7102107787 9 D-5 / DATA AK02CS( 5) / +.1567831810 8523106725 9034899033 3 D-6 / DATA AK02CS( 6) / -.7770110438 5217377103 1579975446 0 D-8 / DATA AK02CS( 7) / +.4611182576 1797178825 3313052958 6 D-9 / DATA AK02CS( 8) / -.3158592997 8605657705 2666580330 9 D-10 / DATA AK02CS( 9) / +.2435018039 3650411278 3588781432 9 D-11 / DATA AK02CS( 10) / -.2074331387 3983478977 0985337350 6 D-12 / DATA AK02CS( 11) / +.1925787280 5899170847 4273650469 3 D-13 / DATA AK02CS( 12) / -.1927554805 8389561036 0034718221 8 D-14 / DATA AK02CS( 13) / +.2062198029 1978182782 8523786964 4 D-15 / DATA AK02CS( 14) / -.2341685117 5792424026 0364019507 1 D-16 / DATA AK02CS( 15) / +.2805902810 6430422468 1517882845 8 D-17 / DATA AK02CS( 16) / -.3530507631 1618079458 1548246357 3 D-18 / DATA AK02CS( 17) / +.4645295422 9351082674 2421633706 6 D-19 / DATA AK02CS( 18) / -.6368625941 3442664739 2205346133 3 D-20 / DATA AK02CS( 19) / +.9069521310 9865155676 2234880000 0 D-21 / DATA AK02CS( 20) / -.1337974785 4236907398 4500531199 9 D-21 / DATA AK02CS( 21) / +.2039836021 8599523155 2208896000 0 D-22 / DATA AK02CS( 22) / -.3207027481 3678405000 6086997333 3 D-23 / DATA AK02CS( 23) / +.5189744413 6623099636 2635946666 6 D-24 / DATA AK02CS( 24) / -.8629501497 5405721929 6460799999 9 D-25 / DATA AK02CS( 25) / +.1472161183 1025598552 0803840000 0 D-25 / DATA AK02CS( 26) / -.2573069023 8670112838 1235199999 9 D-26 / DATA AK02CS( 27) / +.4601774086 6435165873 7664000000 0 D-27 / DATA AK02CS( 28) / -.8411555324 2010937371 3066666666 6 D-28 / DATA AK02CS( 29) / +.1569806306 6353689393 0154666666 6 D-28 / DATA AK02CS( 30) / -.2988226453 0057577889 7919999999 9 D-29 / DATA AK02CS( 31) / +.5796831375 2168365206 1866666666 6 D-30 / DATA AK02CS( 32) / -.1145035994 3476813321 5573333333 3 D-30 / DATA AK02CS( 33) / +.2301266594 2496828020 0533333333 3 D-31 / C C------------------------------------------------------------------- C C Series for BI0 on the interval 0. to 9.00000D+00 C with weighted error 9.51E-34 C log weighted error 33.02 C significant figures required 33.31 C decimal places required 33.65 C DATA BI0 CS( 1) / -.7660547252 8391449510 8189497624 3285 D-1 / DATA BI0 CS( 2) / +.1927337953 9938082699 5240875088 1196 D+1 / DATA BI0 CS( 3) / +.2282644586 9203013389 3702929233 0415 D+0 / DATA BI0 CS( 4) / +.1304891466 7072904280 7933421069 1888 D-1 / DATA BI0 CS( 5) / +.4344270900 8164874513 7868268102 6107 D-3 / DATA BI0 CS( 6) / +.9422657686 0019346639 2317174411 8766 D-5 / DATA BI0 CS( 7) / +.1434006289 5106910799 6209187817 9957 D-6 / DATA BI0 CS( 8) / +.1613849069 6617490699 1541971999 4611 D-8 / DATA BI0 CS( 9) / +.1396650044 5356696994 9509270814 2522 D-10 / DATA BI0 CS( 10) / +.9579451725 5054453446 2752317189 3333 D-13 / DATA BI0 CS( 11) / +.5333981859 8625021310 1510774400 0000 D-15 / DATA BI0 CS( 12) / +.2458716088 4374707746 9678591999 9999 D-17 / DATA BI0 CS( 13) / +.9535680890 2487700269 4434133333 3333 D-20 / DATA BI0 CS( 14) / +.3154382039 7214273367 8933333333 3333 D-22 / DATA BI0 CS( 15) / +.9004564101 0946374314 6666666666 6666 D-25 / DATA BI0 CS( 16) / +.2240647369 1236700160 0000000000 0000 D-27 / DATA BI0 CS( 17) / +.4903034603 2428373333 3333333333 3333 D-30 / DATA BI0 CS( 18) / +.9508172606 1226666666 6666666666 6666 D-33 / C C------------------------------------------------------------------- C DATA NTK0 / 0 / C C***FIRST EXECUTABLE STATEMENT DWK0 C IF (M .LE. 0) GO TO 910 C IF (NTK0 .EQ. 0) THEN EPMACH = D1MACH(3) EPS = 0.10D0*EPMACH NTK0 = IDWCS(BK0CS, LBK0, EPS) NTAK0 = IDWCS(AK0CS, LAK0, EPS) NTAK02 = IDWCS(AK02CS, LAK02, EPS) NTI0 = IDWCS(BI0CS, LBI0, EPS) XSML = SQRT (4.0D0*EPMACH) XMAX = -LOG(D1MACH(1)) XMAX = XMAX - 0.50D0*XMAX*LOG(XMAX)/(XMAX+0.50D0) - 0.010D0 ENDIF C NTOT = 0 C CALL DWNLE(M,X,0.0D0,KEY) IF (KEY .NE. 0) GO TO 920 C C ------------------ C CASE X .GT. XMAX C ------------------ C C NOTE -- K0 UNDERFLOWS FOR X .GT. XMAX C DO 5 I=1,M F(I) = 0.0D0 5 CONTINUE C C ---------------- C CASE X .LE. 2.0 C ---------------- C CALL DWLE(M,X,2.0D0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL DWGTHR(N,X,INDX,XCMP) C C ... COMPUTE I0(X) ... RESULT IN ZCMP C DO 10 J=1,N TCMP(J) = XCMP(J)**2/4.50D0 - 1.0D0 10 CONTINUE CALL DWCS(N,TCMP,BI0CS,NTI0,ZCMP,B0,B1,B2) DO 15 J=1,N ZCMP(J) = 2.750D0 + ZCMP(J) 15 CONTINUE C DO 20 J=1,N TCMP(J) = 0.50D0*XCMP(J)**2 - 1.0D0 20 CONTINUE CALL DWCS(N,TCMP,BK0CS,NTK0,YCMP,B0,B1,B2) DO 30 J=1,N YCMP(J) = -LOG(0.50D0*XCMP(J))*ZCMP(J) - 0.250D0 + YCMP(J) 30 CONTINUE CALL DWSCTR(N,YCMP,INDX,F) ENDIF C C ------------------------- C CASE 2.0 .LT. X .LE. 8.0 C ------------------------- C CALL DWGTLE(M,X,2.0D0,8.0D0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL DWGTHR(N,X,INDX,XCMP) DO 40 J=1,N TCMP(J) = (16.0D0/XCMP(J) - 5.0D0)/3.0D0 40 CONTINUE CALL DWCS(N,TCMP,AK0CS,NTAK0,ZCMP,B0,B1,B2) DO 50 J=1,N YCMP(J) = EXP(-XCMP(J))*(1.250D0 + ZCMP(J))/SQRT(XCMP(J)) 50 CONTINUE CALL DWSCTR(N,YCMP,INDX,F) ENDIF C C -------------------------- C CASE 8.0 .LT. X .LE. XMAX C -------------------------- C CALL DWGTLE(M,X,8.0D0,XMAX,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL DWGTHR(N,X,INDX,XCMP) DO 60 J=1,N TCMP(J) = 16.0D0/XCMP(J) - 1.0D0 60 CONTINUE CALL DWCS(N,TCMP,AK02CS,NTAK02,ZCMP,B0,B1,B2) DO 70 J=1,N YCMP(J) = EXP(-XCMP(J))*(1.250D0 + ZCMP(J))/SQRT(XCMP(J)) 70 CONTINUE CALL DWSCTR(N,YCMP,INDX,F) ENDIF C C ----- C EXITS C ----- C C ... NORMAL C INFO = 0 GO TO 999 C C ... M .LE. 0 C 910 CONTINUE INFO = 1 IF (NTOT .NE. M) INFO = -1 GO TO 999 C C ... X IS ZERO OR NEGATIVE C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END