C C --------------------------------------- C Compare scalar fnlib routines to vfnlib C --------------------------------------- C C Ron Boisvert and Bonita Saunders C Computing and Applied Mathematics Laboratory C National Institute of Standards and Technology C Gaithersburg, MD 20899 C C boisvert@cam.nist.gov C saunders@cam.nist.gov C C 21 Feb 91 C C---------------------------------------------------------------------- C INTEGER M PARAMETER ( M=1000 ) C DOUBLE PRECISION XLO, XHI, X, F, FX, WORK INTEGER I, IWORK, J C DIMENSION XLO(3,8), XHI(3,8), X(M), F(M), FX(M), WORK(7*M), + IWORK(M) C SAVE XLO, XHI C EXTERNAL DVI0, DBESI0, DVI1, DBESI1, DVJ0, DBESJ0, DVJ1, DBESJ1, + DVK0, DBESK0, DVK1, DBESK1, DVY0, DBESY0, DVY1, DBESY1 C DATA ((XLO(I,J),XHI(I,J),I=1,3),J=1,8) / + 0.0D0, 3.0D0 , -8.0D0, -3.0D0, 10.0D0, 50.0D0, + 0.0D0, 3.0D0 , -8.0D0, -3.0D0, 10.0D0, 50.0D0, + 0.0D0, 4.0D0 , -8.0D0, -4.0D0, 8.0D0, 50.0D0, + 0.0D0, 4.0D0 , -8.0D0, -4.0D0, 8.0D0, 50.0D0, + 0.0D0, 2.0D0 , 2.0D0, 8.0D0, 10.0D0, 50.0D0, + 0.0D0, 2.0D0 , 2.0D0, 8.0D0, 10.0D0, 50.0D0, + 0.0D0, 1.0D-3, 1.0D0, 4.0D0, 4.0D0, 50.0D0, + 0.0D0, 1.0D-3, 1.0D0, 4.0D0, 4.0D0, 50.0D0 / C C---------------------------------------------------------------------- C C ... TEST EACH USER-CALLABLE ROUTINE C CALL TEST('I0',DBESI0,DVI0,XLO(1,1),XHI(1,1),M,X,F,FX,WORK,IWORK) CALL TEST('I1',DBESI1,DVI1,XLO(1,2),XHI(1,2),M,X,F,FX,WORK,IWORK) CALL TEST('J0',DBESJ0,DVJ0,XLO(1,3),XHI(1,3),M,X,F,FX,WORK,IWORK) CALL TEST('J1',DBESJ1,DVJ1,XLO(1,4),XHI(1,4),M,X,F,FX,WORK,IWORK) CALL TEST('K0',DBESK0,DVK0,XLO(1,5),XHI(1,5),M,X,F,FX,WORK,IWORK) CALL TEST('K1',DBESK1,DVK1,XLO(1,6),XHI(1,6),M,X,F,FX,WORK,IWORK) CALL TEST('Y0',DBESY0,DVY0,XLO(1,7),XHI(1,7),M,X,F,FX,WORK,IWORK) CALL TEST('Y1',DBESY1,DVY1,XLO(1,8),XHI(1,8),M,X,F,FX,WORK,IWORK) C END SUBROUTINE TEST (NAME, SFUN, VSUB, XLOW, XHIGH, M, X, F, FX, + WORK, IWORK) C C---------------------------------------------------------------------- C C compare scalar function sfun to vectorized subroutine vsub for m C arguments distributed in various ways in the intervals C (xlo(i),xhi(i)), i=1,2,3. C C---------------------------------------------------------------------- C C ---------- C PARAMETERS C ---------- C DOUBLE PRECISION SFUN, XLOW, XHIGH, X, F, FX, WORK INTEGER IWORK, M CHARACTER*(*) NAME C DIMENSION XLOW(3), XHIGH(3), X(M), F(M), FX(M), WORK(6*M), + IWORK(M) C EXTERNAL SFUN, VSUB C C --------------- C LOCAL VARIABLES C --------------- C LOGICAL DEBUG INTEGER NREP, NTRY PARAMETER ( DEBUG=.FALSE., NREP=5, NTRY=19 ) C INTEGER I, INFO, IREP, IOFF, IPC(NTRY), JPC(NTRY), K, KPC, + MA, MB, MC, NFAIL, NREPV DOUBLE PRECISION DIFF, DMAX, DX, EPS, D1MACH REAL SECOND, SU, SUMAX, SUMIN, SUAVE, T0, T1, TS, TSAVE, TV, + TVAVE C SAVE IPC, JPC C DATA IPC / 0, 0, 0, 0, 0, 25, 25, 25, 25, 33, 50, 50, 50, + 75, 75,100, 1, 1, 98/ DATA JPC / 0, 25, 50, 75,100, 0, 25, 50, 75, 33, 0, 25, 50, + 0, 25, 0, 1, 98, 1/ C C---------------------------------------------------------------------- C C ntry = number of trials (each with different argument distribution) C nrep = number of times test repeated to get reliable timings C C the distribution of arguments in trial k is C C ipc(k)% in (xlow(1),xhigh(1)) C jpc(k)% in (xlow(2),xhigh(2)) C 100-ipc(k)-jpc(k)% in (xlow(3),xhigh(3)) C C scalar results (fx) and vector results (f) are compared and accepted C if abs((f - fx)/fx) < eps for fx > 1 (relative difference) C or abs(f - fx) < eps for fx < 1 (absolute difference) C C---------------------------------------------------------------------- C EPS = 5.0D0*D1MACH(4) C ... TRIGGER INITIALIZATIONS C X(1) = 1.0D0 CALL VSUB(1,X,F,WORK,IWORK,INFO) FX(1) = SFUN(X(1)) C PRINT * PRINT * PRINT *,' ----------' PRINT *,' TEST OF ',NAME PRINT *,' ----------' PRINT * PRINT *,' Vector Length = ',M PRINT *,' Repetitions = ',NREP PRINT *,' Eps = ',EPS PRINT * PRINT *,' Range A = (',XLOW(1),',',XHIGH(1),')' PRINT *,' Range B = (',XLOW(2),',',XHIGH(2),')' PRINT *,' Range C = (',XLOW(3),',',XHIGH(3),')' PRINT * PRINT *,' %A = percent arguments in range A' PRINT *,' %B = percent arguments in range B' PRINT *,' %C = percent arguments in range C' PRINT *,' Scalar = time for scalar code (FNLIB)' PRINT *,' Vector = time for vector code (VFNLIB)' PRINT *,' SU = speedup (Scalar/Vector)' PRINT *,' Nfail = number of differences .gt. eps' PRINT *,' Dmax = maximum difference' PRINT * PRINT *,' Times given in seconds/argument.' PRINT * PRINT *,' Differences are abs(s-v)/max(abs(s),1))' PRINT *,' where s = FNLIB result and v = VFNLIB result.' PRINT * PRINT *,' %A %B %C Scalar Vector SU Nfail Dmax' PRINT *,' -------------------------------------------------------' C SUMAX = 0.0 SUMIN = 1000.0 SUAVE = 0.0 TSAVE = 0.0 TVAVE = 0.0 C C ... PERFORM EACH TEST C DO 100 K=1,NTRY C MA = IPC(K)*M/100 MB = JPC(K)*M/100 KPC = 100 - IPC(K) - JPC(K) MC = M - MA - MB C C ... SELECT X VECTOR C IOFF = 0 C IF (MA .GT. 0) THEN DX = (XHIGH(1) - XLOW(1))/MA DO 10 I=1,MA X(I) = XLOW(1) + I*DX 10 CONTINUE IOFF = IOFF + MA ENDIF C IF (MB .GT. 0) THEN DX = (XHIGH(2) - XLOW(2))/MB DO 20 I=1,MB X(I+IOFF) = XLOW(2) + I*DX 20 CONTINUE IOFF = IOFF + MB ENDIF C IF (MC .GT. 0) THEN DX = (XHIGH(3) - XLOW(3))/MC DO 30 I=1,MC X(I+IOFF) = XLOW(3) + I*DX 30 CONTINUE ENDIF C C ... RUN VECTOR CODE C NREPV = 5*NREP T0 = SECOND() DO 40 IREP=1,NREPV CALL VSUB(M,X,F,WORK,IWORK,INFO) 40 CONTINUE T1 = SECOND() TV = (T1 - T0)/NREPV IF (INFO .LT. 0) PRINT *,'Warning in ',NAME,' : Info = ',INFO IF (INFO .GT. 0) PRINT *,'Failure in ',NAME,' : Info = ',INFO, + ' Position = ',IWORK(1) C C ... RUN SCALAR CODE C T0 = SECOND() DO 60 IREP=1,NREP DO 50 I=1,M FX(I) = SFUN(X(I)) 50 CONTINUE 60 CONTINUE T1 = SECOND() TS = (T1 - T0)/NREP C C ... COMPUTE TIMINGS C SU = TS/TV SUMIN = MIN(SU,SUMIN) SUMAX = MAX(SU,SUMAX) SUAVE = SUAVE + SU TSAVE = TSAVE + TS TVAVE = TVAVE + TV C C ... EVALUATE RESULTS C DMAX = 0.0D0 NFAIL = 0 DO 70 I=1,M DIFF = ABS(FX(I)-F(I))/MAX(ABS(FX(I)),1.0D0) DMAX = MAX(DIFF,DMAX) IF (DIFF .GT. EPS) THEN NFAIL = NFAIL + 1 IF (DEBUG) PRINT *,'Failure ',I,X(I),F(I),FX(I),DIFF ENDIF 70 CONTINUE C C ... PRINT RESULTS C PRINT 1000, IPC(K),JPC(K),KPC,TS/M,TV/M,SU,NFAIL,DMAX C 100 CONTINUE C SUAVE = SUAVE/NTRY TSAVE = TSAVE/M/NTRY TVAVE = TVAVE/M/NTRY C PRINT * PRINT 1002,'Ave Scalar = ',TSAVE,' sec/argument' PRINT 1002,'Ave Vector = ',TVAVE,' sec/argument' PRINT * PRINT 1001,'Max speed up = ',SUMAX PRINT 1001,'Min speed up = ',SUMIN PRINT * PRINT 1001,'Ave speed up = ',SUAVE C RETURN 1000 FORMAT(3(1X,I3),1P,2(3X,E7.1),0P,1X,F5.1,1X,I4,3X,1P,E7.1) 1001 FORMAT(1X,A,F5.1,A) 1002 FORMAT(1X,A,1P,E7.1,A) END double precision function dbesi0 (x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bi0cs(18), xmax, xsml, y, d1mach, 1 dcsevl, dbsi0e, dexp, dlog, dsqrt c external d1mach, dbsi0e, dcsevl, dexp, dlog, dsqrt, initds c c series for bi0 on the interval 0. to 9.00000e+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 data nti0, xsml, xmax / 0, 2*0.d0 / c if (nti0.ne.0) go to 10 nti0 = initds (bi0cs, 18, 0.1*sngl(d1mach(3))) xsml = dsqrt (8.0d0*d1mach(3)) xmax = dlog (d1mach(2)) c 10 y = dabs(x) if (y.gt.3.0d0) go to 20 c dbesi0 = 1.0d0 if (y.gt.xsml) dbesi0 = 2.75d0 + dcsevl (y*y/4.5d0-1.d0, bi0cs, 1 nti0) return c 20 if (y.gt.xmax) call seteru ( 1 35hdbesi0 dabs(x) so big i0 overflows, 35, 2, 2) c dbesi0 = dexp(y) * dbsi0e(x) c return end double precision function dbesi1 (x) c oct 1983 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bi1cs(17), xmax, xmin, xsml, y, d1mach, 1 dcsevl, dbsi1e, dexp, dlog, dsqrt c external d1mach, dbsi1e, dcsevl, dexp, dlog, dsqrt, initds c c series for bi1 on the interval 0. to 9.00000e+00 c with weighted error 1.44e-32 c log weighted error 31.84 c significant figures required 31.45 c decimal places required 32.46 c data bi1 cs( 1) / -.1971713261 0998597316 1385032181 49 d-2 / data bi1 cs( 2) / +.4073488766 7546480608 1553936520 14 d+0 / data bi1 cs( 3) / +.3483899429 9959455866 2450377837 87 d-1 / data bi1 cs( 4) / +.1545394556 3001236038 5984010584 89 d-2 / data bi1 cs( 5) / +.4188852109 8377784129 4588320041 20 d-4 / data bi1 cs( 6) / +.7649026764 8362114741 9597039660 69 d-6 / data bi1 cs( 7) / +.1004249392 4741178689 1798080372 38 d-7 / data bi1 cs( 8) / +.9932207791 9238106481 3712980548 63 d-10 / data bi1 cs( 9) / +.7663801791 8447637275 2001716813 49 d-12 / data bi1 cs( 10) / +.4741418923 8167394980 3880919481 60 d-14 / data bi1 cs( 11) / +.2404114404 0745181799 8631720320 00 d-16 / data bi1 cs( 12) / +.1017150500 7093713649 1211007999 99 d-18 / data bi1 cs( 13) / +.3645093565 7866949458 4917333333 33 d-21 / data bi1 cs( 14) / +.1120574950 2562039344 8106666666 66 d-23 / data bi1 cs( 15) / +.2987544193 4468088832 0000000000 00 d-26 / data bi1 cs( 16) / +.6973231093 9194709333 3333333333 33 d-29 / data bi1 cs( 17) / +.1436794822 0620800000 0000000000 00 d-31 / c data nti1, xmin, xsml, xmax / 0, 3*0.d0 / c if (nti1.ne.0) go to 10 nti1 = initds (bi1cs, 17, 0.1*sngl(d1mach(3))) xmin = 2.0d0*d1mach(1) xsml = dsqrt (8.0d0*d1mach(3)) xmax = dlog (d1mach(2)) c 10 y = dabs(x) if (y.gt.3.0d0) go to 20 c dbesi1 = 0.0d0 if (y.eq.0.0d0) return c if (y.le.xmin) call seteru ( 1 38hdbesi1 dabs(x) so small i1 underflows, 38, 1, 0) if (y.gt.xmin) dbesi1 = 0.5d0*x if (y.gt.xsml) dbesi1 = x*(0.875d0 + dcsevl (y*y/4.5d0-1.d0, 1 bi1cs, nti1)) return c 20 if (y.gt.xmax) call seteru ( 1 35hdbesi1 dabs(x) so big i1 overflows, 35, 2, 2) c dbesi1 = dexp(y) * dbsi1e(x) c return end double precision function dbesj0 (x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bj0cs(19), ampl, theta, xsml, y, 1 d1mach, dcsevl, dcos, dsqrt c external d1mach, dcos, dcsevl, dsqrt, initds c c series for bj0 on the interval 0. to 1.60000e+01 c with weighted error 4.39e-32 c log weighted error 31.36 c significant figures required 31.21 c decimal places required 32.00 c data bj0 cs( 1) / +.1002541619 6893913701 0731272640 74 d+0 / data bj0 cs( 2) / -.6652230077 6440513177 6787578311 24 d+0 / data bj0 cs( 3) / +.2489837034 9828131370 4604687266 80 d+0 / data bj0 cs( 4) / -.3325272317 0035769653 8843415038 54 d-1 / data bj0 cs( 5) / +.2311417930 4694015462 9049241177 29 d-2 / data bj0 cs( 6) / -.9911277419 9508092339 0485193365 49 d-4 / data bj0 cs( 7) / +.2891670864 3998808884 7339037470 78 d-5 / data bj0 cs( 8) / -.6121085866 3032635057 8184074815 16 d-7 / data bj0 cs( 9) / +.9838650793 8567841324 7687486364 15 d-9 / data bj0 cs( 10) / -.1242355159 7301765145 5158970068 36 d-10 / data bj0 cs( 11) / +.1265433630 2559045797 9158272103 63 d-12 / data bj0 cs( 12) / -.1061945649 5287244546 9148175129 59 d-14 / data bj0 cs( 13) / +.7470621075 8024567437 0989155840 00 d-17 / data bj0 cs( 14) / -.4469703227 4412780547 6270079999 99 d-19 / data bj0 cs( 15) / +.2302428158 4337436200 5230933333 33 d-21 / data bj0 cs( 16) / -.1031914479 4166698148 5226666666 66 d-23 / data bj0 cs( 17) / +.4060817827 4873322700 8000000000 00 d-26 / data bj0 cs( 18) / -.1414383600 5240913919 9999999999 99 d-28 / data bj0 cs( 19) / +.4391090549 6698880000 0000000000 00 d-31 / c data ntj0, xsml / 0, 0.d0 / c if (ntj0.ne.0) go to 10 ntj0 = initds (bj0cs, 19, 0.1*sngl(d1mach(3))) xsml = dsqrt (4.0d0*d1mach(3)) c 10 y = dabs(x) if (y.gt.4.0d0) go to 20 c dbesj0 = 1.0d0 if (y.gt.xsml) dbesj0 = dcsevl (.125d0*y*y-1.d0, bj0cs, ntj0) return c 20 call d9b0mp (y, ampl, theta) dbesj0 = ampl * dcos(theta) c return end double precision function dbesj1 (x) c sept 1983 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bj1cs(19), ampl, theta, xsml, xmin, y, 1 d1mach, dcsevl, dcos, dsqrt c external d1mach, dcos, dcsevl, dsqrt, initds c c series for bj1 on the interval 0. to 1.60000e+01 c with weighted error 1.16e-33 c log weighted error 32.93 c significant figures required 32.36 c decimal places required 33.57 c data bj1 cs( 1) / -.1172614151 3332786560 6240574524 003 d+0 / data bj1 cs( 2) / -.2536152183 0790639562 3030884554 698 d+0 / data bj1 cs( 3) / +.5012708098 4469568505 3656363203 743 d-1 / data bj1 cs( 4) / -.4631514809 6250819184 2619728789 772 d-2 / data bj1 cs( 5) / +.2479962294 1591402453 9124064592 364 d-3 / data bj1 cs( 6) / -.8678948686 2788258452 1246435176 416 d-5 / data bj1 cs( 7) / +.2142939171 4379369150 2766250991 292 d-6 / data bj1 cs( 8) / -.3936093079 1831797922 9322764073 061 d-8 / data bj1 cs( 9) / +.5591182317 9468800401 8248059864 032 d-10 / data bj1 cs( 10) / -.6327616404 6613930247 7695274014 880 d-12 / data bj1 cs( 11) / +.5840991610 8572470032 6945563268 266 d-14 / data bj1 cs( 12) / -.4482533818 7012581903 9135059199 999 d-16 / data bj1 cs( 13) / +.2905384492 6250246630 6018688000 000 d-18 / data bj1 cs( 14) / -.1611732197 8414416541 2118186666 666 d-20 / data bj1 cs( 15) / +.7739478819 3927463729 8346666666 666 d-23 / data bj1 cs( 16) / -.3248693782 1119984114 3466666666 666 d-25 / data bj1 cs( 17) / +.1202237677 2274102272 0000000000 000 d-27 / data bj1 cs( 18) / -.3952012212 6513493333 3333333333 333 d-30 / data bj1 cs( 19) / +.1161678082 2664533333 3333333333 333 d-32 / c data ntj1, xsml, xmin / 0, 2*0.d0 / c if (ntj1.ne.0) go to 10 ntj1 = initds (bj1cs, 19, 0.1*sngl(d1mach(3))) c xsml = dsqrt (4.0d0*d1mach(3)) xmin = 2.0d0*d1mach(1) c 10 y = dabs(x) if (y.gt.4.0d0) go to 20 c dbesj1 = 0.0d0 if (y.eq.0.0d0) return c if (y.le.xmin) call seteru ( 1 38hdbesj1 dabs(x) so small j1 underflows, 38, 1, 0) if (y.gt.xmin) dbesj1 = 0.5d0*x if (y.gt.xsml) dbesj1 = x*(.25d0 + dcsevl (.125d0*y*y-1.d0, 1 bj1cs, ntj1) ) return c 20 call d9b1mp (y, ampl, theta) C**Next line changed from***************** C dbesj1 = ampl * dcos(theta) C to dbesj1 = sign(ampl,x) * dcos(theta) C**by Ron Boisvert, 19 Mar 91 ************ c return end double precision function dbesk0 (x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bk0cs(16), xmax, xsml, y, 1 d1mach, dcsevl, dbesi0, dbsk0e, dexp, dlog, dsqrt c external d1mach, dbesi0, dbsk0e, dcsevl, dexp, dlog, dsqrt, initds c c series for bk0 on the interval 0. to 4.00000e+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 data ntk0, xsml, xmax / 0, 2*0.d0 / c if (ntk0.ne.0) go to 10 ntk0 = initds (bk0cs, 16, 0.1*sngl(d1mach(3))) xsml = dsqrt (4.0d0*d1mach(3)) xmax = -dlog(d1mach(1)) xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0) c 10 if (x.le.0.d0) call seteru (29hdbesk0 x is zero or negative, 1 29, 2, 2) if (x.gt.2.0d0) go to 20 c y = 0.d0 if (x.gt.xsml) y = x*x dbesk0 = -dlog(0.5d0*x)*dbesi0(x) - 0.25d0 + dcsevl (.5d0*y-1.d0, 1 bk0cs, ntk0) return c 20 dbesk0 = 0.d0 if (x.gt.xmax) call seteru (30hdbesk0 x so big k0 underflows, 1 30, 1, 0) if (x.gt.xmax) return c dbesk0 = dexp(-x) * dbsk0e(x) c return end double precision function dbesk1 (x) c july 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bk1cs(16), xmax, xmin, xsml, y, 1 d1mach, dcsevl, dbesi1, dbsk1e, dexp, dlog, dsqrt c external d1mach, dbesi1, dbsk1e, dcsevl, dexp, dlog, dsqrt, initds c c series for bk1 on the interval 0. to 4.00000e+00 c with weighted error 9.16e-32 c log weighted error 31.04 c significant figures required 30.61 c decimal places required 31.64 c data bk1 cs( 1) / +.2530022733 8947770532 5311208685 33 d-1 / data bk1 cs( 2) / -.3531559607 7654487566 7238316918 01 d+0 / data bk1 cs( 3) / -.1226111808 2265714823 4790679300 42 d+0 / data bk1 cs( 4) / -.6975723859 6398643501 8129202960 83 d-2 / data bk1 cs( 5) / -.1730288957 5130520630 1765073689 79 d-3 / data bk1 cs( 6) / -.2433406141 5659682349 6007350301 64 d-5 / data bk1 cs( 7) / -.2213387630 7347258558 3152525451 26 d-7 / data bk1 cs( 8) / -.1411488392 6335277610 9583302126 08 d-9 / data bk1 cs( 9) / -.6666901694 1993290060 8537512643 73 d-12 / data bk1 cs( 10) / -.2427449850 5193659339 2631968648 53 d-14 / data bk1 cs( 11) / -.7023863479 3862875971 7837971200 00 d-17 / data bk1 cs( 12) / -.1654327515 5100994675 4910293333 33 d-19 / data bk1 cs( 13) / -.3233834745 9944491991 8933333333 33 d-22 / data bk1 cs( 14) / -.5331275052 9265274999 4666666666 66 d-25 / data bk1 cs( 15) / -.7513040716 2157226666 6666666666 66 d-28 / data bk1 cs( 16) / -.9155085717 6541866666 6666666666 66 d-31 / c data ntk1, xmin, xsml, xmax / 0, 3*0.d0 / c if (ntk1.ne.0) go to 10 ntk1 = initds (bk1cs, 16, 0.1*sngl(d1mach(3))) xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0) xsml = dsqrt (4.0d0*d1mach(3)) xmax = -dlog(d1mach(1)) xmax = xmax - 0.5d0*xmax*dlog(xmax)/(xmax+0.5d0) - 0.01d0 c 10 if (x.le.0.d0) call seteru (29hdbesk1 x is zero or negative, 1 29, 2, 2) if (x.gt.2.0d0) go to 20 c if (x.lt.xmin) call seteru (31hdbesk1 x so small k1 overflows, 1 31, 3, 2) y = 0.d0 if (x.gt.xsml) y = x*x dbesk1 = dlog(0.5d0*x)*dbesi1(x) + (0.75d0 + dcsevl (.5d0*y-1.d0, 1 bk1cs, ntk1))/x return c 20 dbesk1 = 0.d0 if (x.gt.xmax) call seteru (30hdbesk1 x so big k1 underflows, 1 30, 1, 0) if (x.gt.xmax) return c dbesk1 = dexp(-x) * dbsk1e(x) c return end double precision function dbesy0 (x) c august 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision x, by0cs(19), ampl, theta, twodpi, xsml, 1 y, alnhaf, d1mach, dcsevl, dbesj0, dlog, dsin, dsqrt c external d1mach, dbesj0, dcsevl, dlog, dsin, dsqrt, initds c c series for by0 on the interval 0. to 1.60000e+01 c with weighted error 8.14e-32 c log weighted error 31.09 c significant figures required 30.31 c decimal places required 31.73 c data by0 cs( 1) / -.1127783939 2865573217 9398054602 8 d-1 / data by0 cs( 2) / -.1283452375 6042034604 8088453183 8 d+0 / data by0 cs( 3) / -.1043788479 9794249365 8176227661 8 d+0 / data by0 cs( 4) / +.2366274918 3969695409 2415926461 3 d-1 / data by0 cs( 5) / -.2090391647 7004862391 9622395034 2 d-2 / data by0 cs( 6) / +.1039754539 3905725209 9924657638 1 d-3 / data by0 cs( 7) / -.3369747162 4239720967 1877534503 7 d-5 / data by0 cs( 8) / +.7729384267 6706671585 2136721637 1 d-7 / data by0 cs( 9) / -.1324976772 6642595914 4347606896 4 d-8 / data by0 cs( 10) / +.1764823261 5404527921 0038936315 8 d-10 / data by0 cs( 11) / -.1881055071 5801962006 0282301206 9 d-12 / data by0 cs( 12) / +.1641865485 3661495027 9223718574 9 d-14 / data by0 cs( 13) / -.1195659438 6046060857 4599100672 0 d-16 / data by0 cs( 14) / +.7377296297 4401858424 9411242666 6 d-19 / data by0 cs( 15) / -.3906843476 7104373307 4090666666 6 d-21 / data by0 cs( 16) / +.1795503664 4361579498 2912000000 0 d-23 / data by0 cs( 17) / -.7229627125 4480104789 3333333333 3 d-26 / data by0 cs( 18) / +.2571727931 6351685973 3333333333 3 d-28 / data by0 cs( 19) / -.8141268814 1636949333 3333333333 3 d-31 / c data twodpi / 0.6366197723 6758134307 5535053490 057 d0 / data alnhaf /-0.6931471805 5994530941 7232121458 18d0 / data nty0, xsml / 0, 0.d0 / c if (nty0.ne.0) go to 10 nty0 = initds (by0cs, 19, 0.1*sngl(d1mach(3))) xsml = dsqrt (4.0d0*d1mach(3)) c 10 if (x.le.0.d0) call seteru (29hdbesy0 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0d0) go to 20 c y = 0.d0 if (x.gt.xsml) y = x*x dbesy0 = twodpi*(alnhaf+dlog(x))*dbesj0(x) + .375d0 1 + dcsevl (.125d0*y-1.d0, by0cs, nty0) return c 20 call d9b0mp (x, ampl, theta) dbesy0 = ampl * dsin(theta) return c end double precision function dbesy1 (x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, by1cs(20), ampl, theta, twodpi, xmin, xsml, 1 y, d1mach, dcsevl, dbesj1, dexp, dlog, dsin, dsqrt c external d1mach, dbesj1, dcsevl, dexp, dlog, dsin, dsqrt, initds c c series for by1 on the interval 0. to 1.60000e+01 c with weighted error 8.65e-33 c log weighted error 32.06 c significant figures required 32.17 c decimal places required 32.71 c data by1 cs( 1) / +.3208047100 6119086293 2352018628 015 d-1 / data by1 cs( 2) / +.1262707897 4335004495 3431725999 727 d+1 / data by1 cs( 3) / +.6499961899 9231750009 7490637314 144 d-2 / data by1 cs( 4) / -.8936164528 8605041165 3144160009 712 d-1 / data by1 cs( 5) / +.1325088122 1757095451 2375510370 043 d-1 / data by1 cs( 6) / -.8979059119 6483523775 3039508298 105 d-3 / data by1 cs( 7) / +.3647361487 9583067824 2287368165 349 d-4 / data by1 cs( 8) / -.1001374381 6660005554 9075523845 295 d-5 / data by1 cs( 9) / +.1994539657 3901739703 1159372421 243 d-7 / data by1 cs( 10) / -.3023065601 8033816728 4799332520 743 d-9 / data by1 cs( 11) / +.3609878156 9478119611 6252914242 474 d-11 / data by1 cs( 12) / -.3487488297 2875824241 4552947409 066 d-13 / data by1 cs( 13) / +.2783878971 5591766581 3507698517 333 d-15 / data by1 cs( 14) / -.1867870968 6194876876 6825352533 333 d-17 / data by1 cs( 15) / +.1068531533 9116825975 7070336000 000 d-19 / data by1 cs( 16) / -.5274721956 6844822894 3872000000 000 d-22 / data by1 cs( 17) / +.2270199403 1556641437 0133333333 333 d-24 / data by1 cs( 18) / -.8595390353 9452310869 3333333333 333 d-27 / data by1 cs( 19) / +.2885404379 8337945600 0000000000 000 d-29 / data by1 cs( 20) / -.8647541138 9371733333 3333333333 333 d-32 / c data twodpi / 0.6366197723 6758134307 5535053490 057 d0 / data nty1, xmin, xsml / 0, 2*0.d0 / c if (nty1.ne.0) go to 10 nty1 = initds (by1cs, 20, 0.1*sngl(d1mach(3))) c xmin = 1.571d0 * dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 1 0.01d0) xsml = dsqrt (4.0d0*d1mach(3)) c 10 if (x.le.0.d0) call seteru (29hdbesy1 x is zero or negative, 29, 1 1, 2) if (x.gt.4.0d0) go to 20 c if (x.lt.xmin) call seteru (31hdbesy1 x so small y1 overflows, 1 31, 3, 2) y = 0.d0 if (x.gt.xsml) y = x*x dbesy1 = twodpi * dlog(0.5d0*x)*dbesj1(x) + (0.5d0 + 1 dcsevl (.125d0*y-1.d0, by1cs, nty1))/x return c 20 call d9b1mp (x, ampl, theta) dbesy1 = ampl * dsin(theta) return c end function initds (dos, nos, eta) c june 1977 edition. w. fullerton, c3, los alamos scientific lab. c c initialize the double precision orthogonal series dos so that initds c is the number of terms needed to insure the error is no larger than c eta. ordinarily eta will be chosen to be one-tenth machine precision. c c input arguments -- c dos dble prec array of nos coefficients in an orthogonal series. c nos number of coefficients in dos. c eta requested accuracy of series. c double precision dos(nos) c if (nos.lt.1) call seteru ( 1 35hinitds number of coefficients lt 1, 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(sngl(dos(i))) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call seteru (28hinitds eta may be too small, 28, 1 1, 2) initds = i c return end double precision function dcsevl (x, a, n) c c evaluate the n-term chebyshev series a at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). c c input arguments -- c x dble prec value at which the series is to be evaluated. c a dble prec array of n terms of a chebyshev series. in eval- c uating a, only half the first coef is summed. c n number of terms in array a. c double precision a(n), x, twox, b0, b1, b2 c if (n.lt.1) call seteru (28hdcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hdcsevl number of terms gt 1000, 1 31, 3, 2) if (x.lt.(-1.1d0) .or. x.gt.1.1d0) call seteru ( 1 25hdcsevl x outside (-1,+1), 25, 1, 1) c twox = 2.0d0*x b1 = 0.d0 b0 = 0.d0 do 10 i=1,n b2 = b1 b1 = b0 ni = n - i + 1 b0 = twox*b1 - b2 + a(ni) 10 continue c dcsevl = 0.5d0 * (b0-b2) c return end subroutine seteru (messg, nmessg, nerr, iopt) common /cseter/ iunflo integer messg(1) data iunflo / 0 / c if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt) if (iopt.ne.0) return c if (iunflo.le.0) return call seterr (messg, nmessg, nerr, 1) c return end double precision function dbsi0e (x) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bi0cs(18), ai0cs(46), ai02cs(69), 1 xsml, y, d1mach, dcsevl, dexp, dsqrt c external d1mach, dcsevl, dexp, dsqrt, initds c c series for bi0 on the interval 0. to 9.00000e+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 series for ai0 on the interval 1.25000e-01 to 3.33333e-01 c with weighted error 2.74e-32 c log weighted error 31.56 c significant figures required 30.15 c decimal places required 32.39 c data ai0 cs( 1) / +.7575994494 0237959427 2987203743 8 d-1 / data ai0 cs( 2) / +.7591380810 8233455072 9297873320 4 d-2 / data ai0 cs( 3) / +.4153131338 9237505018 6319749138 2 d-3 / data ai0 cs( 4) / +.1070076463 4390730735 8242970217 0 d-4 / data ai0 cs( 5) / -.7901179979 2128946607 5031948573 0 d-5 / data ai0 cs( 6) / -.7826143501 4387522697 8898980690 9 d-6 / data ai0 cs( 7) / +.2783849942 9488708063 8118538985 7 d-6 / data ai0 cs( 8) / +.8252472600 6120271919 6682913319 8 d-8 / data ai0 cs( 9) / -.1204463945 5201991790 5496089110 3 d-7 / data ai0 cs( 10) / +.1559648598 5060764436 1228752792 8 d-8 / data ai0 cs( 11) / +.2292556367 1033165434 7725480285 7 d-9 / data ai0 cs( 12) / -.1191622884 2790646036 7777423447 8 d-9 / data ai0 cs( 13) / +.1757854916 0324098302 1833124774 3 d-10 / data ai0 cs( 14) / +.1128224463 2189005171 4441135682 4 d-11 / data ai0 cs( 15) / -.1146848625 9272988777 2963387698 2 d-11 / data ai0 cs( 16) / +.2715592054 8036628726 4365192160 6 d-12 / data ai0 cs( 17) / -.2415874666 5626878384 4247572028 1 d-13 / data ai0 cs( 18) / -.6084469888 2551250646 0609963922 4 d-14 / data ai0 cs( 19) / +.3145705077 1754772937 0836026730 3 d-14 / data ai0 cs( 20) / -.7172212924 8711877179 6217505917 6 d-15 / data ai0 cs( 21) / +.7874493403 4541033960 8390960332 7 d-16 / data ai0 cs( 22) / +.1004802753 0094624023 4524457183 9 d-16 / data ai0 cs( 23) / -.7566895365 3505348534 2843588881 0 d-17 / data ai0 cs( 24) / +.2150380106 8761198878 1205128784 5 d-17 / data ai0 cs( 25) / -.3754858341 8308744291 5158445260 8 d-18 / data ai0 cs( 26) / +.2354065842 2269925769 0075710532 2 d-19 / data ai0 cs( 27) / +.1114667612 0479285302 2637335511 0 d-19 / data ai0 cs( 28) / -.5398891884 3969903786 9677932270 9 d-20 / data ai0 cs( 29) / +.1439598792 2407526770 4285840452 2 d-20 / data ai0 cs( 30) / -.2591916360 1110934064 6081840196 2 d-21 / data ai0 cs( 31) / +.2238133183 9985839074 3409229824 0 d-22 / data ai0 cs( 32) / +.5250672575 3647711727 7221683199 9 d-23 / data ai0 cs( 33) / -.3249904138 5332307841 7343228586 6 d-23 / data ai0 cs( 34) / +.9924214103 2050379278 5728471040 0 d-24 / data ai0 cs( 35) / -.2164992254 2446695231 4655429973 3 d-24 / data ai0 cs( 36) / +.3233609471 9435940839 7333299199 9 d-25 / data ai0 cs( 37) / -.1184620207 3967424898 2473386666 6 d-26 / data ai0 cs( 38) / -.1281671853 9504986505 4833868799 9 d-26 / data ai0 cs( 39) / +.5827015182 2793905116 0556885333 3 d-27 / data ai0 cs( 40) / -.1668222326 0261097193 6450150399 9 d-27 / data ai0 cs( 41) / +.3625309510 5415699757 0068480000 0 d-28 / data ai0 cs( 42) / -.5733627999 0557135899 4595839999 9 d-29 / data ai0 cs( 43) / +.3736796722 0630982296 4258133333 3 d-30 / data ai0 cs( 44) / +.1602073983 1568519633 6551253333 3 d-30 / data ai0 cs( 45) / -.8700424864 0572298845 2249599999 9 d-31 / data ai0 cs( 46) / +.2741320937 9374811456 0341333333 3 d-31 / c c series for ai02 on the interval 0. to 1.25000e-01 c with weighted error 1.97e-32 c log weighted error 31.71 c significant figures required 30.15 c decimal places required 32.63 c data ai02cs( 1) / +.5449041101 4108831607 8960962268 0 d-1 / data ai02cs( 2) / +.3369116478 2556940898 9785662979 9 d-2 / data ai02cs( 3) / +.6889758346 9168239842 6263914301 1 d-4 / data ai02cs( 4) / +.2891370520 8347564829 6692402323 2 d-5 / data ai02cs( 5) / +.2048918589 4690637418 2760534093 1 d-6 / data ai02cs( 6) / +.2266668990 4981780645 9327743136 1 d-7 / data ai02cs( 7) / +.3396232025 7083863451 5084396952 3 d-8 / data ai02cs( 8) / +.4940602388 2249695891 0482449783 5 d-9 / data ai02cs( 9) / +.1188914710 7846438342 4084525196 3 d-10 / data ai02cs( 10) / -.3149916527 9632413645 3864862961 9 d-10 / data ai02cs( 11) / -.1321581184 0447713118 7540739926 7 d-10 / data ai02cs( 12) / -.1794178531 5068061177 7943574026 9 d-11 / data ai02cs( 13) / +.7180124451 3836662336 7106429346 9 d-12 / data ai02cs( 14) / +.3852778382 7421427011 4089801777 6 d-12 / data ai02cs( 15) / +.1540086217 5214098269 1325823339 7 d-13 / data ai02cs( 16) / -.4150569347 2872220866 2689972015 6 d-13 / data ai02cs( 17) / -.9554846698 8283076487 0214494312 5 d-14 / data ai02cs( 18) / +.3811680669 3526224207 4605535511 8 d-14 / data ai02cs( 19) / +.1772560133 0565263836 0493266675 8 d-14 / data ai02cs( 20) / -.3425485619 6772191346 1924790328 2 d-15 / data ai02cs( 21) / -.2827623980 5165834849 4205593759 4 d-15 / data ai02cs( 22) / +.3461222867 6974610930 9706250813 4 d-16 / data ai02cs( 23) / +.4465621420 2967599990 1042054284 3 d-16 / data ai02cs( 24) / -.4830504485 9441820712 5525403795 4 d-17 / data ai02cs( 25) / -.7233180487 8747539545 6227240924 5 d-17 / data ai02cs( 26) / +.9921475412 1736985988 8046093981 0 d-18 / data ai02cs( 27) / +.1193650890 8459820855 0439949924 2 d-17 / data ai02cs( 28) / -.2488709837 1508072357 2054491660 2 d-18 / data ai02cs( 29) / -.1938426454 1609059289 8469781132 6 d-18 / data ai02cs( 30) / +.6444656697 3734438687 8301949394 9 d-19 / data ai02cs( 31) / +.2886051596 2892243264 8171383073 4 d-19 / data ai02cs( 32) / -.1601954907 1749718070 6167156200 7 d-19 / data ai02cs( 33) / -.3270815010 5923147208 9193567485 9 d-20 / data ai02cs( 34) / +.3686932283 8264091811 4600723939 3 d-20 / data ai02cs( 35) / +.1268297648 0309501530 1359529710 9 d-22 / data ai02cs( 36) / -.7549825019 3772739076 9636664410 1 d-21 / data ai02cs( 37) / +.1502133571 3778353496 3712789053 4 d-21 / data ai02cs( 38) / +.1265195883 5096485349 3208799248 3 d-21 / data ai02cs( 39) / -.6100998370 0836807086 2940891600 2 d-22 / data ai02cs( 40) / -.1268809629 2601282643 6872095924 2 d-22 / data ai02cs( 41) / +.1661016099 8907414578 4038487490 5 d-22 / data ai02cs( 42) / -.1585194335 7658855793 7970504881 4 d-23 / data ai02cs( 43) / -.3302645405 9682178009 5381766755 6 d-23 / data ai02cs( 44) / +.1313580902 8392397817 4039623117 4 d-23 / data ai02cs( 45) / +.3689040246 6711567933 1425637280 4 d-24 / data ai02cs( 46) / -.4210141910 4616891492 1978247249 9 d-24 / data ai02cs( 47) / +.4791954591 0828657806 3171401373 0 d-25 / data ai02cs( 48) / +.8459470390 2218217952 9971707412 4 d-25 / data ai02cs( 49) / -.4039800940 8728324931 4607937181 0 d-25 / data ai02cs( 50) / -.6434714653 6504313473 0100850469 5 d-26 / data ai02cs( 51) / +.1225743398 8756659903 4464736990 5 d-25 / data ai02cs( 52) / -.2934391316 0257089231 9879821175 4 d-26 / data ai02cs( 53) / -.1961311309 1949829262 0371205728 9 d-26 / data ai02cs( 54) / +.1503520374 8221934241 6229900309 8 d-26 / data ai02cs( 55) / -.9588720515 7448265520 3386388206 9 d-28 / data ai02cs( 56) / -.3483339380 8170454863 9441108511 4 d-27 / data ai02cs( 57) / +.1690903610 2630436730 6244960725 6 d-27 / data ai02cs( 58) / +.1982866538 7356030438 9400115718 8 d-28 / data ai02cs( 59) / -.5317498081 4918162145 7583002528 4 d-28 / data ai02cs( 60) / +.1803306629 8883929462 3501450390 1 d-28 / data ai02cs( 61) / +.6213093341 4548931758 8405311242 2 d-29 / data ai02cs( 62) / -.7692189292 7721618632 0072806673 0 d-29 / data ai02cs( 63) / +.1858252826 1117025426 2556016596 3 d-29 / data ai02cs( 64) / +.1237585142 2813957248 9927154554 1 d-29 / data ai02cs( 65) / -.1102259120 4092238032 1779478779 2 d-29 / data ai02cs( 66) / +.1886287118 0397044900 7787447943 1 d-30 / data ai02cs( 67) / +.2160196872 2436589131 4903141406 0 d-30 / data ai02cs( 68) / -.1605454124 9197432005 8446594965 5 d-30 / data ai02cs( 69) / +.1965352984 5942906039 3884807331 8 d-31 / c data nti0, ntai0, ntai02, xsml / 3*0, 0.d0 / c if (nti0.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) nti0 = initds (bi0cs, 18, eta) ntai0 = initds (ai0cs, 46, eta) ntai02 = initds (ai02cs, 69, eta) xsml = dsqrt (8.0d0*d1mach(3)) c 10 y = dabs(x) if (y.gt.3.0d0) go to 20 c dbsi0e = 1.0d0 if (y.gt.xsml) dbsi0e = dexp(-y) * (2.75d0 + 1 dcsevl (y*y/4.5d0-1.d0, bi0cs, nti0) ) return c 20 if (y.le.8.d0) dbsi0e = (0.375d0 + dcsevl ((48.d0/y-11.d0)/5.d0, 1 ai0cs, ntai0))/dsqrt(y) if (y.gt.8.d0) dbsi0e = (0.375d0 + dcsevl (16.d0/y-1.d0, ai02cs, 1 ntai02))/dsqrt(y) c return end double precision function dbsi1e (x) c oct 1983 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bi1cs(17), ai1cs(46), ai12cs(69), xmin, 1 xsml, y, d1mach, dcsevl, dexp, dsqrt c external d1mach, dcsevl, dexp, dsqrt, initds c c series for bi1 on the interval 0. to 9.00000e+00 c with weighted error 1.44e-32 c log weighted error 31.84 c significant figures required 31.45 c decimal places required 32.46 c data bi1 cs( 1) / -.1971713261 0998597316 1385032181 49 d-2 / data bi1 cs( 2) / +.4073488766 7546480608 1553936520 14 d+0 / data bi1 cs( 3) / +.3483899429 9959455866 2450377837 87 d-1 / data bi1 cs( 4) / +.1545394556 3001236038 5984010584 89 d-2 / data bi1 cs( 5) / +.4188852109 8377784129 4588320041 20 d-4 / data bi1 cs( 6) / +.7649026764 8362114741 9597039660 69 d-6 / data bi1 cs( 7) / +.1004249392 4741178689 1798080372 38 d-7 / data bi1 cs( 8) / +.9932207791 9238106481 3712980548 63 d-10 / data bi1 cs( 9) / +.7663801791 8447637275 2001716813 49 d-12 / data bi1 cs( 10) / +.4741418923 8167394980 3880919481 60 d-14 / data bi1 cs( 11) / +.2404114404 0745181799 8631720320 00 d-16 / data bi1 cs( 12) / +.1017150500 7093713649 1211007999 99 d-18 / data bi1 cs( 13) / +.3645093565 7866949458 4917333333 33 d-21 / data bi1 cs( 14) / +.1120574950 2562039344 8106666666 66 d-23 / data bi1 cs( 15) / +.2987544193 4468088832 0000000000 00 d-26 / data bi1 cs( 16) / +.6973231093 9194709333 3333333333 33 d-29 / data bi1 cs( 17) / +.1436794822 0620800000 0000000000 00 d-31 / c c series for ai1 on the interval 1.25000e-01 to 3.33333e-01 c with weighted error 2.81e-32 c log weighted error 31.55 c significant figures required 29.93 c decimal places required 32.38 c data ai1 cs( 1) / -.2846744181 8814786741 0037246830 7 d-1 / data ai1 cs( 2) / -.1922953231 4432206510 4444877497 9 d-1 / data ai1 cs( 3) / -.6115185857 9437889822 5624991778 5 d-3 / data ai1 cs( 4) / -.2069971253 3502277088 8282377797 9 d-4 / data ai1 cs( 5) / +.8585619145 8107255655 3694467313 8 d-5 / data ai1 cs( 6) / +.1049498246 7115908625 1745399786 0 d-5 / data ai1 cs( 7) / -.2918338918 4479022020 9343232669 7 d-6 / data ai1 cs( 8) / -.1559378146 6317390001 6068096907 7 d-7 / data ai1 cs( 9) / +.1318012367 1449447055 2530287390 9 d-7 / data ai1 cs( 10) / -.1448423418 1830783176 3913446781 5 d-8 / data ai1 cs( 11) / -.2908512243 9931420948 2504099301 0 d-9 / data ai1 cs( 12) / +.1266388917 8753823873 1115969040 3 d-9 / data ai1 cs( 13) / -.1664947772 9192206706 2417839858 0 d-10 / data ai1 cs( 14) / -.1666653644 6094329760 9593715499 9 d-11 / data ai1 cs( 15) / +.1242602414 2907682652 3216847201 7 d-11 / data ai1 cs( 16) / -.2731549379 6724323972 5146142863 3 d-12 / data ai1 cs( 17) / +.2023947881 6458037807 0026268898 1 d-13 / data ai1 cs( 18) / +.7307950018 1168836361 9869812612 3 d-14 / data ai1 cs( 19) / -.3332905634 4046749438 1377861713 3 d-14 / data ai1 cs( 20) / +.7175346558 5129537435 4225466567 0 d-15 / data ai1 cs( 21) / -.6982530324 7962563558 5062922365 6 d-16 / data ai1 cs( 22) / -.1299944201 5627607600 6044608058 7 d-16 / data ai1 cs( 23) / +.8120942864 2427988920 5467834286 0 d-17 / data ai1 cs( 24) / -.2194016207 4107368981 5626664378 3 d-17 / data ai1 cs( 25) / +.3630516170 0296548482 7986093233 4 d-18 / data ai1 cs( 26) / -.1695139772 4391041663 0686679039 9 d-19 / data ai1 cs( 27) / -.1288184829 8979078071 1688253822 2 d-19 / data ai1 cs( 28) / +.5694428604 9670527801 0999107310 9 d-20 / data ai1 cs( 29) / -.1459597009 0904800565 4550990028 7 d-20 / data ai1 cs( 30) / +.2514546010 6757173140 8469133448 5 d-21 / data ai1 cs( 31) / -.1844758883 1391248181 6040002901 3 d-22 / data ai1 cs( 32) / -.6339760596 2279486419 2860979199 9 d-23 / data ai1 cs( 33) / +.3461441102 0310111111 0814662656 0 d-23 / data ai1 cs( 34) / -.1017062335 3713935475 9654102357 3 d-23 / data ai1 cs( 35) / +.2149877147 0904314459 6250077866 6 d-24 / data ai1 cs( 36) / -.3045252425 2386764017 4620617386 6 d-25 / data ai1 cs( 37) / +.5238082144 7212859821 7763498666 6 d-27 / data ai1 cs( 38) / +.1443583107 0893824464 1678950399 9 d-26 / data ai1 cs( 39) / -.6121302074 8900427332 0067071999 9 d-27 / data ai1 cs( 40) / +.1700011117 4678184183 4918980266 6 d-27 / data ai1 cs( 41) / -.3596589107 9842441585 3521578666 6 d-28 / data ai1 cs( 42) / +.5448178578 9484185766 5051306666 6 d-29 / data ai1 cs( 43) / -.2731831789 6890849891 6256426666 6 d-30 / data ai1 cs( 44) / -.1858905021 7086007157 7190399999 9 d-30 / data ai1 cs( 45) / +.9212682974 5139334411 2776533333 3 d-31 / data ai1 cs( 46) / -.2813835155 6535611063 7083306666 6 d-31 / c c series for ai12 on the interval 0. to 1.25000e-01 c with weighted error 1.83e-32 c log weighted error 31.74 c significant figures required 29.97 c decimal places required 32.66 c data ai12cs( 1) / +.2857623501 8280120474 4984594846 9 d-1 / data ai12cs( 2) / -.9761097491 3614684077 6516445730 2 d-2 / data ai12cs( 3) / -.1105889387 6262371629 1256921277 5 d-3 / data ai12cs( 4) / -.3882564808 8776903934 5654477627 4 d-5 / data ai12cs( 5) / -.2512236237 8702089252 9452002212 1 d-6 / data ai12cs( 6) / -.2631468846 8895195068 3705236523 2 d-7 / data ai12cs( 7) / -.3835380385 9642370220 4500678796 8 d-8 / data ai12cs( 8) / -.5589743462 1965838068 6811252222 9 d-9 / data ai12cs( 9) / -.1897495812 3505412344 9892503323 8 d-10 / data ai12cs( 10) / +.3252603583 0154882385 5508067994 9 d-10 / data ai12cs( 11) / +.1412580743 6613781331 6336633284 6 d-10 / data ai12cs( 12) / +.2035628544 1470895072 2452613684 0 d-11 / data ai12cs( 13) / -.7198551776 2459085120 9258989044 6 d-12 / data ai12cs( 14) / -.4083551111 0921973182 2849963969 1 d-12 / data ai12cs( 15) / -.2101541842 7726643130 1984572746 2 d-13 / data ai12cs( 16) / +.4272440016 7119513542 9778833699 7 d-13 / data ai12cs( 17) / +.1042027698 4128802764 1741449994 8 d-13 / data ai12cs( 18) / -.3814403072 4370078047 6707253539 6 d-14 / data ai12cs( 19) / -.1880354775 5107824485 1273453396 3 d-14 / data ai12cs( 20) / +.3308202310 9209282827 3190335240 5 d-15 / data ai12cs( 21) / +.2962628997 6459501390 6854654205 2 d-15 / data ai12cs( 22) / -.3209525921 9934239587 7837353288 7 d-16 / data ai12cs( 23) / -.4650305368 4893583255 7128281897 9 d-16 / data ai12cs( 24) / +.4414348323 0717079499 4611375964 1 d-17 / data ai12cs( 25) / +.7517296310 8421048054 2545808029 5 d-17 / data ai12cs( 26) / -.9314178867 3268833756 8484784515 7 d-18 / data ai12cs( 27) / -.1242193275 1948909561 1678448869 7 d-17 / data ai12cs( 28) / +.2414276719 4548484690 0515390217 6 d-18 / data ai12cs( 29) / +.2026944384 0532851789 7192286069 2 d-18 / data ai12cs( 30) / -.6394267188 2690977870 4391988681 1 d-19 / data ai12cs( 31) / -.3049812452 3730958960 8488450357 1 d-19 / data ai12cs( 32) / +.1612841851 6514802251 3462230769 1 d-19 / data ai12cs( 33) / +.3560913964 3099250545 1027090462 0 d-20 / data ai12cs( 34) / -.3752017947 9364390796 6682800324 6 d-20 / data ai12cs( 35) / -.5787037427 0747993459 5198231074 1 d-22 / data ai12cs( 36) / +.7759997511 6481619619 8236963209 2 d-21 / data ai12cs( 37) / -.1452790897 2022333940 6445987408 5 d-21 / data ai12cs( 38) / -.1318225286 7390367021 2192275337 4 d-21 / data ai12cs( 39) / +.6116654862 9030707018 7999133171 7 d-22 / data ai12cs( 40) / +.1376279762 4271264277 3024338363 4 d-22 / data ai12cs( 41) / -.1690837689 9593478849 1983938230 6 d-22 / data ai12cs( 42) / +.1430596088 5954331539 8720108538 5 d-23 / data ai12cs( 43) / +.3409557828 0905940204 0536772990 2 d-23 / data ai12cs( 44) / -.1309457666 2707602278 4573872642 4 d-23 / data ai12cs( 45) / -.3940706411 2402574360 9352141755 7 d-24 / data ai12cs( 46) / +.4277137426 9808765808 0616679735 2 d-24 / data ai12cs( 47) / -.4424634830 9826068819 0028312302 9 d-25 / data ai12cs( 48) / -.8734113196 2307149721 1530978874 7 d-25 / data ai12cs( 49) / +.4045401335 6835333921 4340414242 8 d-25 / data ai12cs( 50) / +.7067100658 0946894656 5160771780 6 d-26 / data ai12cs( 51) / -.1249463344 5651052230 0286451860 5 d-25 / data ai12cs( 52) / +.2867392244 4034370329 7948339142 6 d-26 / data ai12cs( 53) / +.2044292892 5042926702 8177957421 0 d-26 / data ai12cs( 54) / -.1518636633 8204625683 7134680291 1 d-26 / data ai12cs( 55) / +.8110181098 1875758861 3227910703 7 d-28 / data ai12cs( 56) / +.3580379354 7735860911 2717370327 0 d-27 / data ai12cs( 57) / -.1692929018 9279025095 9305717544 8 d-27 / data ai12cs( 58) / -.2222902499 7024276390 6775852777 4 d-28 / data ai12cs( 59) / +.5424535127 1459696550 4860040112 8 d-28 / data ai12cs( 60) / -.1787068401 5780186887 6491299330 4 d-28 / data ai12cs( 61) / -.6565479068 7228149388 2392943788 0 d-29 / data ai12cs( 62) / +.7807013165 0611452809 2206770683 9 d-29 / data ai12cs( 63) / -.1816595260 6689797173 7933315222 1 d-29 / data ai12cs( 64) / -.1287704952 6600848203 7687559895 9 d-29 / data ai12cs( 65) / +.1114548172 9881645474 1370927369 4 d-29 / data ai12cs( 66) / -.1808343145 0393369391 5936887668 7 d-30 / data ai12cs( 67) / -.2231677718 2037719522 3244822893 9 d-30 / data ai12cs( 68) / +.1619029596 0803415106 1790980361 4 d-30 / data ai12cs( 69) / -.1834079908 8049414139 0130843921 0 d-31 / c data nti1, ntai1, ntai12, xmin, xsml / 3*0, 2*0.d0 / c if (nti1.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) nti1 = initds (bi1cs, 17, eta) ntai1 = initds (ai1cs, 46, eta) ntai12 = initds (ai12cs, 69, eta) c xmin = 2.0d0*d1mach(1) xsml = dsqrt (8.0d0*d1mach(3)) c 10 y = dabs(x) if (y.gt.3.0d0) go to 20 c dbsi1e = 0.0d0 if (y.eq.0.0d0) return c if (y.le.xmin) call seteru ( 1 38hdbsi1e dabs(x) so small i1 underflows, 37, 1, 0) if (y.gt.xmin) dbsi1e = 0.5d0*x if (y.gt.xsml) dbsi1e = x*(0.875d0 + dcsevl (y*y/4.5d0-1.d0, 1 bi1cs, nti1) ) dbsi1e = dexp(-y) * dbsi1e return c 20 if (y.le.8.d0) dbsi1e = (0.375d0 + dcsevl ((48.d0/y-11.d0)/5.d0, 1 ai1cs, ntai1))/dsqrt(y) if (y.gt.8.d0) dbsi1e = (0.375d0 + dcsevl (16.d0/y-1.d0, ai12cs, 1 ntai12))/dsqrt(y) dbsi1e = dsign (dbsi1e, x) c return end subroutine d9b0mp (x, ampl, theta) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. c c evaluate the modulus and phase for the bessel j0 and y0 functions. c double precision x, ampl, theta, bm0cs(37), bt02cs(39), 1 bm02cs(40), bth0cs(44), xmax, pi4, z, d1mach, dcsevl, 2 dsqrt c external d1mach, dcsevl, dsqrt, initds c c series for bm0 on the interval 1.56250e-02 to 6.25000e-02 c with weighted error 4.40e-32 c log weighted error 31.36 c significant figures required 30.02 c decimal places required 32.14 c data bm0 cs( 1) / +.9211656246 8277427125 7376773018 2 d-1 / data bm0 cs( 2) / -.1050590997 2719051024 8071637175 5 d-2 / data bm0 cs( 3) / +.1470159840 7687597540 5639285095 2 d-4 / data bm0 cs( 4) / -.5058557606 0385542233 4792932770 2 d-6 / data bm0 cs( 5) / +.2787254538 6324441766 3035613788 1 d-7 / data bm0 cs( 6) / -.2062363611 7809148026 1884101897 3 d-8 / data bm0 cs( 7) / +.1870214313 1388796751 3817259626 1 d-9 / data bm0 cs( 8) / -.1969330971 1356362002 4173077782 5 d-10 / data bm0 cs( 9) / +.2325973793 9992754440 1250881805 2 d-11 / data bm0 cs( 10) / -.3009520344 9382502728 5122473448 2 d-12 / data bm0 cs( 11) / +.4194521333 8506691814 7120676864 6 d-13 / data bm0 cs( 12) / -.6219449312 1884458259 7326742956 4 d-14 / data bm0 cs( 13) / +.9718260411 3360684696 0176588526 9 d-15 / data bm0 cs( 14) / -.1588478585 7010752073 6663596693 7 d-15 / data bm0 cs( 15) / +.2700072193 6713088900 8621732445 8 d-16 / data bm0 cs( 16) / -.4750092365 2340089924 7750478677 3 d-17 / data bm0 cs( 17) / +.8615128162 6043708731 9170374656 0 d-18 / data bm0 cs( 18) / -.1605608686 9561448157 4560270335 9 d-18 / data bm0 cs( 19) / +.3066513987 3144829751 8853980159 9 d-19 / data bm0 cs( 20) / -.5987764223 1939564306 9650561706 6 d-20 / data bm0 cs( 21) / +.1192971253 7482483064 8906984106 6 d-20 / data bm0 cs( 22) / -.2420969142 0448054894 8468258133 3 d-21 / data bm0 cs( 23) / +.4996751760 5106164533 7100287999 9 d-22 / data bm0 cs( 24) / -.1047493639 3511585100 9504051199 9 d-22 / data bm0 cs( 25) / +.2227786843 7974681010 4818346666 6 d-23 / data bm0 cs( 26) / -.4801813239 3981628623 7054293333 3 d-24 / data bm0 cs( 27) / +.1047962723 4709599564 7699626666 6 d-24 / data bm0 cs( 28) / -.2313858165 6786153251 0126080000 0 d-25 / data bm0 cs( 29) / +.5164823088 4626742116 3519999999 9 d-26 / data bm0 cs( 30) / -.1164691191 8500653895 2540159999 9 d-26 / data bm0 cs( 31) / +.2651788486 0433192829 5833600000 0 d-27 / data bm0 cs( 32) / -.6092559503 8257284976 9130666666 6 d-28 / data bm0 cs( 33) / +.1411804686 1442593080 3882666666 6 d-28 / data bm0 cs( 34) / -.3298094961 2317372457 5061333333 3 d-29 / data bm0 cs( 35) / +.7763931143 0740650317 1413333333 3 d-30 / data bm0 cs( 36) / -.1841031343 6614584784 2133333333 3 d-30 / data bm0 cs( 37) / +.4395880138 5943107371 0079999999 9 d-31 / c c series for bth0 on the interval 0. to 1.56250e-02 c with weighted error 2.66e-32 c log weighted error 31.57 c significant figures required 30.67 c decimal places required 32.40 c data bth0cs( 1) / -.2490178086 2128936717 7097937899 67 d+0 / data bth0cs( 2) / +.4855029960 9623749241 0486155354 85 d-3 / data bth0cs( 3) / -.5451183734 5017204950 6562735635 05 d-5 / data bth0cs( 4) / +.1355867305 9405964054 3774459299 03 d-6 / data bth0cs( 5) / -.5569139890 2227626227 5832184149 20 d-8 / data bth0cs( 6) / +.3260903182 4994335304 0042057194 68 d-9 / data bth0cs( 7) / -.2491880786 2461341125 2379038779 93 d-10 / data bth0cs( 8) / +.2344937742 0882520554 3524135648 91 d-11 / data bth0cs( 9) / -.2609653444 4310387762 1775747661 36 d-12 / data bth0cs( 10) / +.3335314042 0097395105 8699550149 23 d-13 / data bth0cs( 11) / -.4789000044 0572684646 7507705574 09 d-14 / data bth0cs( 12) / +.7595617843 6192215972 6425685452 48 d-15 / data bth0cs( 13) / -.1313155601 6891440382 7733974876 33 d-15 / data bth0cs( 14) / +.2448361834 5240857495 4268207383 55 d-16 / data bth0cs( 15) / -.4880572981 0618777683 2567619183 31 d-17 / data bth0cs( 16) / +.1032728502 9786316149 2237563612 04 d-17 / data bth0cs( 17) / -.2305763381 5057217157 0047445270 25 d-18 / data bth0cs( 18) / +.5404444300 1892693993 0171084837 65 d-19 / data bth0cs( 19) / -.1324069519 4366572724 1550328823 85 d-19 / data bth0cs( 20) / +.3378079562 1371970203 4247921247 22 d-20 / data bth0cs( 21) / -.8945762915 7111779003 0269262922 99 d-21 / data bth0cs( 22) / +.2451990688 9219317090 8999086514 05 d-21 / data bth0cs( 23) / -.6938842287 6866318680 1399331576 57 d-22 / data bth0cs( 24) / +.2022827871 4890138392 9463033377 91 d-22 / data bth0cs( 25) / -.6062850000 2335483105 7941953717 64 d-23 / data bth0cs( 26) / +.1864974896 4037635381 8237883962 70 d-23 / data bth0cs( 27) / -.5878373238 4849894560 2450365308 67 d-24 / data bth0cs( 28) / +.1895859144 7999563485 5311795035 13 d-24 / data bth0cs( 29) / -.6248197937 2258858959 2916207285 65 d-25 / data bth0cs( 30) / +.2101790168 4551024686 6386335290 74 d-25 / data bth0cs( 31) / -.7208430093 5209253690 8139339924 46 d-26 / data bth0cs( 32) / +.2518136389 2474240867 1564059767 46 d-26 / data bth0cs( 33) / -.8951804225 8785778806 1439459536 43 d-27 / data bth0cs( 34) / +.3235723747 9762298533 2562358685 87 d-27 / data bth0cs( 35) / -.1188301051 9855353657 0471441137 96 d-27 / data bth0cs( 36) / +.4430628690 7358104820 5792319417 31 d-28 / data bth0cs( 37) / -.1676100964 8834829495 7920101356 81 d-28 / data bth0cs( 38) / +.6429294692 1207466972 5323939660 88 d-29 / data bth0cs( 39) / -.2499226116 6978652421 2072136827 63 d-29 / data bth0cs( 40) / +.9839979429 9521955672 8282603553 18 d-30 / data bth0cs( 41) / -.3922037524 2408016397 9891316261 58 d-30 / data bth0cs( 42) / +.1581810703 0056522138 5906188456 92 d-30 / data bth0cs( 43) / -.6452550614 4890715944 3440983654 26 d-31 / data bth0cs( 44) / +.2661111136 9199356137 1770183463 67 d-31 / c c series for bm02 on the interval 0. to 1.56250e-02 c with weighted error 4.72e-32 c log weighted error 31.33 c significant figures required 30.00 c decimal places required 32.13 c data bm02cs( 1) / +.9500415145 2283813693 3086133556 0 d-1 / data bm02cs( 2) / -.3801864682 3656709917 4808156685 1 d-3 / data bm02cs( 3) / +.2258339301 0314811929 5182992722 4 d-5 / data bm02cs( 4) / -.3895725802 3722287647 3062141260 5 d-7 / data bm02cs( 5) / +.1246886416 5120816979 3099052972 5 d-8 / data bm02cs( 6) / -.6065949022 1025037798 0383505838 7 d-10 / data bm02cs( 7) / +.4008461651 4217469910 1527597104 5 d-11 / data bm02cs( 8) / -.3350998183 3980942184 6729879457 4 d-12 / data bm02cs( 9) / +.3377119716 5174173670 6326434199 6 d-13 / data bm02cs( 10) / -.3964585901 6350127005 6935629582 3 d-14 / data bm02cs( 11) / +.5286111503 8838572173 8793974473 5 d-15 / data bm02cs( 12) / -.7852519083 4508523136 5464024349 3 d-16 / data bm02cs( 13) / +.1280300573 3866822010 1163407344 9 d-16 / data bm02cs( 14) / -.2263996296 3914297762 8709924488 4 d-17 / data bm02cs( 15) / +.4300496929 6567903886 4641029047 7 d-18 / data bm02cs( 16) / -.8705749805 1325870797 4753545145 5 d-19 / data bm02cs( 17) / +.1865862713 9620951411 8144277205 0 d-19 / data bm02cs( 18) / -.4210482486 0930654573 4508697230 1 d-20 / data bm02cs( 19) / +.9956676964 2284009915 8162741784 2 d-21 / data bm02cs( 20) / -.2457357442 8053133596 0592147854 7 d-21 / data bm02cs( 21) / +.6307692160 7620315680 8735370705 9 d-22 / data bm02cs( 22) / -.1678773691 4407401426 9333117238 8 d-22 / data bm02cs( 23) / +.4620259064 6739044337 7087813608 7 d-23 / data bm02cs( 24) / -.1311782266 8603087322 3769340249 6 d-23 / data bm02cs( 25) / +.3834087564 1163028277 4792244027 6 d-24 / data bm02cs( 26) / -.1151459324 0777412710 7261329357 6 d-24 / data bm02cs( 27) / +.3547210007 5233385230 7697134521 3 d-25 / data bm02cs( 28) / -.1119218385 8150046462 6435594217 6 d-25 / data bm02cs( 29) / +.3611879427 6298378316 9840499425 7 d-26 / data bm02cs( 30) / -.1190687765 9133331500 9264176246 3 d-26 / data bm02cs( 31) / +.4005094059 4039681318 0247644953 6 d-27 / data bm02cs( 32) / -.1373169422 4522123905 9519391601 7 d-27 / data bm02cs( 33) / +.4794199088 7425315859 9649152643 7 d-28 / data bm02cs( 34) / -.1702965627 6241095840 0699447645 2 d-28 / data bm02cs( 35) / +.6149512428 9363300715 0357516132 4 d-29 / data bm02cs( 36) / -.2255766896 5818283499 4430023724 2 d-29 / data bm02cs( 37) / +.8399707509 2942994860 6165835320 0 d-30 / data bm02cs( 38) / -.3172997595 5626023555 6742393615 2 d-30 / data bm02cs( 39) / +.1215205298 8812985545 8333302651 4 d-30 / data bm02cs( 40) / -.4715852749 7544386930 1321056804 5 d-31 / c c series for bt02 on the interval 1.56250e-02 to 6.25000e-02 c with weighted error 2.99e-32 c log weighted error 31.52 c significant figures required 30.61 c decimal places required 32.32 c data bt02cs( 1) / -.2454829521 3424597462 0504672493 24 d+0 / data bt02cs( 2) / +.1254412103 9084615780 7853317782 99 d-2 / data bt02cs( 3) / -.3125395041 4871522854 9734467095 71 d-4 / data bt02cs( 4) / +.1470977824 9940831164 4534269693 14 d-5 / data bt02cs( 5) / -.9954348893 7950033643 4688503511 58 d-7 / data bt02cs( 6) / +.8549316673 3203041247 5787113977 51 d-8 / data bt02cs( 7) / -.8698975952 6554334557 9855121791 92 d-9 / data bt02cs( 8) / +.1005209953 3559791084 5401010821 53 d-9 / data bt02cs( 9) / -.1282823060 1708892903 4836236855 44 d-10 / data bt02cs( 10) / +.1773170078 1805131705 6557504510 23 d-11 / data bt02cs( 11) / -.2617457456 9485577488 6362841809 25 d-12 / data bt02cs( 12) / +.4082835138 9972059621 9664812211 03 d-13 / data bt02cs( 13) / -.6675166823 9742720054 6067495542 61 d-14 / data bt02cs( 14) / +.1136576139 3071629448 3924695499 51 d-14 / data bt02cs( 15) / -.2005118962 0647160250 5592664121 17 d-15 / data bt02cs( 16) / +.3649797879 4766269635 7205914641 06 d-16 / data bt02cs( 17) / -.6830963756 4582303169 3558437888 00 d-17 / data bt02cs( 18) / +.1310758314 5670756620 0571042679 46 d-17 / data bt02cs( 19) / -.2572336310 1850607778 7571306495 99 d-18 / data bt02cs( 20) / +.5152165744 1863959925 2677809493 33 d-19 / data bt02cs( 21) / -.1051301756 3758802637 9407414613 33 d-19 / data bt02cs( 22) / +.2182038199 1194813847 3010845013 33 d-20 / data bt02cs( 23) / -.4600470121 0362160577 2259054933 33 d-21 / data bt02cs( 24) / +.9840700692 5466818520 9536511999 99 d-22 / data bt02cs( 25) / -.2133403803 5728375844 7359863466 66 d-22 / data bt02cs( 26) / +.4683103642 3973365296 0662869333 33 d-23 / data bt02cs( 27) / -.1040021369 1985747236 5133823999 99 d-23 / data bt02cs( 28) / +.2334910567 7301510051 7777408000 00 d-24 / data bt02cs( 29) / -.5295682532 3318615788 0497493333 33 d-25 / data bt02cs( 30) / +.1212634195 2959756829 1962879999 99 d-25 / data bt02cs( 31) / -.2801889708 2289428760 2756266666 66 d-26 / data bt02cs( 32) / +.6529267898 7012873342 5937066666 66 d-27 / data bt02cs( 33) / -.1533798006 1873346427 8357333333 33 d-27 / data bt02cs( 34) / +.3630588430 6364536682 3594666666 66 d-28 / data bt02cs( 35) / -.8656075571 3629122479 1722666666 66 d-29 / data bt02cs( 36) / +.2077990997 2536284571 2383999999 99 d-29 / data bt02cs( 37) / -.5021117022 1417221674 3253333333 33 d-30 / data bt02cs( 38) / +.1220836027 9441714184 1919999999 99 d-30 / data bt02cs( 39) / -.2986005626 7039913454 2506666666 66 d-31 / c data pi4 / 0.7853981633 9744830961 5660845819 876 d0 / data nbm0, nbt02, nbm02, nbth0, xmax / 4*0, 0.d0 / c if (nbm0.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) nbm0 = initds (bm0cs, 37, eta) nbt02 = initds (bt02cs, 39, eta) nbm02 = initds (bm02cs, 40, eta) nbth0 = initds (bth0cs, 44, eta) c xmax = 1.0d0/d1mach(4) c 10 if (x.lt.4.d0) call seteru (22hd9b0mp x must be ge 4, 22, 1, 2) c if (x.gt.8.d0) go to 20 z = (128.d0/(x*x) - 5.d0)/3.d0 ampl = (.75d0 + dcsevl (z, bm0cs, nbm0))/dsqrt(x) theta = x - pi4 + dcsevl (z, bt02cs, nbt02)/x return c 20 if (x.gt.xmax) call seteru ( 1 37hd9b0mp no precision because x is big, 37, 2, 2) c z = 128.d0/(x*x) - 1.d0 ampl = (.75d0 + dcsevl (z, bm02cs, nbm02))/dsqrt(x) theta = x - pi4 + dcsevl (z, bth0cs, nbth0)/x return c end subroutine d9b1mp (x, ampl, theta) c july 1977 edition. w. fullerton, c3, los alamos scientific lab. c c evaluate the modulus and phase for the bessel j1 and y1 functions. c double precision x, ampl, theta, bm1cs(37), bt12cs(39), 1 bm12cs(40), bth1cs(44), xmax, pi4, z, d1mach, dcsevl, 2 dsqrt c external d1mach, dcsevl, dsqrt, initds c c series for bm1 on the interval 1.56250e-02 to 6.25000e-02 c with weighted error 4.91e-32 c log weighted error 31.31 c significant figures required 30.04 c decimal places required 32.09 c data bm1 cs( 1) / +.1069845452 6180630149 6998530853 8 d+0 / data bm1 cs( 2) / +.3274915039 7159649007 2905514344 5 d-2 / data bm1 cs( 3) / -.2987783266 8316985920 3044577793 8 d-4 / data bm1 cs( 4) / +.8331237177 9919745313 9322266902 3 d-6 / data bm1 cs( 5) / -.4112665690 3020073048 9638172549 8 d-7 / data bm1 cs( 6) / +.2855344228 7892152207 1975766316 1 d-8 / data bm1 cs( 7) / -.2485408305 4156238780 6002659605 5 d-9 / data bm1 cs( 8) / +.2543393338 0725824427 4248439717 4 d-10 / data bm1 cs( 9) / -.2941045772 8229675234 8975082790 9 d-11 / data bm1 cs( 10) / +.3743392025 4939033092 6505615362 6 d-12 / data bm1 cs( 11) / -.5149118293 8211672187 2054824352 7 d-13 / data bm1 cs( 12) / +.7552535949 8651439080 3404076419 9 d-14 / data bm1 cs( 13) / -.1169409706 8288464441 6629062246 4 d-14 / data bm1 cs( 14) / +.1896562449 4347915717 2182460506 0 d-15 / data bm1 cs( 15) / -.3201955368 6932864206 6477531639 4 d-16 / data bm1 cs( 16) / +.5599548399 3162041144 8416990549 3 d-17 / data bm1 cs( 17) / -.1010215894 7304324431 1939044454 4 d-17 / data bm1 cs( 18) / +.1873844985 7275629833 0204271957 3 d-18 / data bm1 cs( 19) / -.3563537470 3285802192 7430143999 9 d-19 / data bm1 cs( 20) / +.6931283819 9712383304 2276351999 9 d-20 / data bm1 cs( 21) / -.1376059453 4065001522 5140893013 3 d-20 / data bm1 cs( 22) / +.2783430784 1070802205 9977932799 9 d-21 / data bm1 cs( 23) / -.5727595364 3205616893 4866943999 9 d-22 / data bm1 cs( 24) / +.1197361445 9188926725 3575679999 9 d-22 / data bm1 cs( 25) / -.2539928509 8918719766 4144042666 6 d-23 / data bm1 cs( 26) / +.5461378289 6572959730 6961919999 9 d-24 / data bm1 cs( 27) / -.1189211341 7733202889 8628949333 3 d-24 / data bm1 cs( 28) / +.2620150977 3400815949 5782400000 0 d-25 / data bm1 cs( 29) / -.5836810774 2556859019 2093866666 6 d-26 / data bm1 cs( 30) / +.1313743500 0805957734 2361599999 9 d-26 / data bm1 cs( 31) / -.2985814622 5103803553 3277866666 6 d-27 / data bm1 cs( 32) / +.6848390471 3346049376 2559999999 9 d-28 / data bm1 cs( 33) / -.1584401568 2224767211 9296000000 0 d-28 / data bm1 cs( 34) / +.3695641006 5709380543 0101333333 3 d-29 / data bm1 cs( 35) / -.8687115921 1446682430 1226666666 6 d-30 / data bm1 cs( 36) / +.2057080846 1587634629 2906666666 6 d-30 / data bm1 cs( 37) / -.4905225761 1162255185 2373333333 3 d-31 / c c series for bt12 on the interval 1.56250e-02 to 6.25000e-02 c with weighted error 3.33e-32 c log weighted error 31.48 c significant figures required 31.05 c decimal places required 32.27 c data bt12cs( 1) / +.7382386012 8742974662 6208397927 64 d+0 / data bt12cs( 2) / -.3336111317 4483906384 4701476811 89 d-2 / data bt12cs( 3) / +.6146345488 8046964698 5148994201 86 d-4 / data bt12cs( 4) / -.2402458516 1602374264 9776354695 68 d-5 / data bt12cs( 5) / +.1466355557 7509746153 2105919972 04 d-6 / data bt12cs( 6) / -.1184191730 5589180567 0051475049 83 d-7 / data bt12cs( 7) / +.1157419896 3919197052 1254663030 55 d-8 / data bt12cs( 8) / -.1300116112 9439187449 3660077945 71 d-9 / data bt12cs( 9) / +.1624539114 1361731937 7421662736 67 d-10 / data bt12cs( 10) / -.2208963682 1403188752 1554417701 28 d-11 / data bt12cs( 11) / +.3218030425 8553177090 4743586537 78 d-12 / data bt12cs( 12) / -.4965314793 2768480785 5520211353 81 d-13 / data bt12cs( 13) / +.8043890043 2847825985 5588826393 17 d-14 / data bt12cs( 14) / -.1358912131 0161291384 6947126822 82 d-14 / data bt12cs( 15) / +.2381050439 7147214869 6765296059 73 d-15 / data bt12cs( 16) / -.4308146636 3849106724 4712414207 99 d-16 / data bt12cs( 17) / +.8020254403 2771002434 9935125504 00 d-17 / data bt12cs( 18) / -.1531631064 2462311864 2300274687 99 d-17 / data bt12cs( 19) / +.2992860635 2715568924 0730405546 66 d-18 / data bt12cs( 20) / -.5970996465 8085443393 8156366506 66 d-19 / data bt12cs( 21) / +.1214028966 9415185024 1608526506 66 d-19 / data bt12cs( 22) / -.2511511469 6612948901 0069777066 66 d-20 / data bt12cs( 23) / +.5279056717 0328744850 7383807999 99 d-21 / data bt12cs( 24) / -.1126050922 7550498324 3611613866 66 d-21 / data bt12cs( 25) / +.2434827735 9576326659 6634624000 00 d-22 / data bt12cs( 26) / -.5331726123 6931800130 0384426666 66 d-23 / data bt12cs( 27) / +.1181361505 9707121039 2059903999 99 d-23 / data bt12cs( 28) / -.2646536828 3353523514 8567893333 33 d-24 / data bt12cs( 29) / +.5990339404 1361503945 5778133333 33 d-25 / data bt12cs( 30) / -.1369085463 0829503109 1363839999 99 d-25 / data bt12cs( 31) / +.3157679015 4380228326 4136533333 33 d-26 / data bt12cs( 32) / -.7345791508 2084356491 4005333333 33 d-27 / data bt12cs( 33) / +.1722808148 0722747930 7059200000 00 d-27 / data bt12cs( 34) / -.4071690796 1286507941 0688000000 00 d-28 / data bt12cs( 35) / +.9693474513 6779622700 3733333333 33 d-29 / data bt12cs( 36) / -.2323763633 7765716765 3546666666 66 d-29 / data bt12cs( 37) / +.5607451067 3522029406 8906666666 66 d-30 / data bt12cs( 38) / -.1361646539 1539005860 5226666666 66 d-30 / data bt12cs( 39) / +.3326310923 3894654388 9066666666 66 d-31 / c c series for bm12 on the interval 0. to 1.56250e-02 c with weighted error 5.01e-32 c log weighted error 31.30 c significant figures required 29.99 c decimal places required 32.10 c data bm12cs( 1) / +.9807979156 2330500272 7209354693 7 d-1 / data bm12cs( 2) / +.1150961189 5046853061 7548348460 2 d-2 / data bm12cs( 3) / -.4312482164 3382054098 8935809773 2 d-5 / data bm12cs( 4) / +.5951839610 0888163078 1302980183 2 d-7 / data bm12cs( 5) / -.1704844019 8269098574 0070158647 8 d-8 / data bm12cs( 6) / +.7798265413 6111095086 5817382740 1 d-10 / data bm12cs( 7) / -.4958986126 7664158094 9175495186 5 d-11 / data bm12cs( 8) / +.4038432416 4211415168 3820226514 4 d-12 / data bm12cs( 9) / -.3993046163 7251754457 6548384664 5 d-13 / data bm12cs( 10) / +.4619886183 1189664943 1334243277 5 d-14 / data bm12cs( 11) / -.6089208019 0953833013 4547261933 3 d-15 / data bm12cs( 12) / +.8960930916 4338764821 5704804124 9 d-16 / data bm12cs( 13) / -.1449629423 9420231229 1651891892 5 d-16 / data bm12cs( 14) / +.2546463158 5377760561 6514964806 8 d-17 / data bm12cs( 15) / -.4809472874 6478364442 5926371862 0 d-18 / data bm12cs( 16) / +.9687684668 2925990490 8727583912 4 d-19 / data bm12cs( 17) / -.2067213372 2779660232 4503811755 1 d-19 / data bm12cs( 18) / +.4646651559 1503847318 0276780959 0 d-20 / data bm12cs( 19) / -.1094966128 8483341382 4135132833 9 d-20 / data bm12cs( 20) / +.2693892797 2886828609 0570761278 5 d-21 / data bm12cs( 21) / -.6894992910 9303744778 1897002685 7 d-22 / data bm12cs( 22) / +.1830268262 7520629098 9066855474 0 d-22 / data bm12cs( 23) / -.5025064246 3519164281 5611355322 4 d-23 / data bm12cs( 24) / +.1423545194 4548060396 3169363419 4 d-23 / data bm12cs( 25) / -.4152191203 6164503880 6888676980 1 d-24 / data bm12cs( 26) / +.1244609201 5039793258 8233007654 7 d-24 / data bm12cs( 27) / -.3827336370 5693042994 3191866128 6 d-25 / data bm12cs( 28) / +.1205591357 8156175353 7472398183 5 d-25 / data bm12cs( 29) / -.3884536246 3764880764 3185936112 4 d-26 / data bm12cs( 30) / +.1278689528 7204097219 0489528346 1 d-26 / data bm12cs( 31) / -.4295146689 4479462720 6193691591 2 d-27 / data bm12cs( 32) / +.1470689117 8290708864 5680270798 3 d-27 / data bm12cs( 33) / -.5128315665 1060731281 8037401779 6 d-28 / data bm12cs( 34) / +.1819509585 4711693854 8143737328 6 d-28 / data bm12cs( 35) / -.6563031314 8419808676 1863505037 3 d-29 / data bm12cs( 36) / +.2404898976 9199606531 9891487583 4 d-29 / data bm12cs( 37) / -.8945966744 6906124732 3495824297 9 d-30 / data bm12cs( 38) / +.3376085160 6572310266 3714897824 0 d-30 / data bm12cs( 39) / -.1291791454 6206563609 1309991696 6 d-30 / data bm12cs( 40) / +.5008634462 9588105206 8495150125 4 d-31 / c c series for bth1 on the interval 0. to 1.56250e-02 c with weighted error 2.82e-32 c log weighted error 31.55 c significant figures required 31.12 c decimal places required 32.37 c data bth1cs( 1) / +.7474995720 3587276055 4434839696 95 d+0 / data bth1cs( 2) / -.1240077714 4651711252 5457775413 84 d-2 / data bth1cs( 3) / +.9925244240 4424527376 6414976895 92 d-5 / data bth1cs( 4) / -.2030369073 7159711052 4193753756 08 d-6 / data bth1cs( 5) / +.7535961770 5690885712 1840175836 29 d-8 / data bth1cs( 6) / -.4166161271 5343550107 6300238562 28 d-9 / data bth1cs( 7) / +.3070161807 0834890481 2451020912 16 d-10 / data bth1cs( 8) / -.2817849963 7605213992 3240088839 24 d-11 / data bth1cs( 9) / +.3079069673 9040295476 0281468216 47 d-12 / data bth1cs( 10) / -.3880330026 2803434112 7873475547 81 d-13 / data bth1cs( 11) / +.5509603960 8630904934 5617262085 62 d-14 / data bth1cs( 12) / -.8659006076 8383779940 1033989539 94 d-15 / data bth1cs( 13) / +.1485604914 1536749003 4236890606 83 d-15 / data bth1cs( 14) / -.2751952981 5904085805 3712121250 09 d-16 / data bth1cs( 15) / +.5455079609 0481089625 0362236409 23 d-17 / data bth1cs( 16) / -.1148653450 1983642749 5436310271 77 d-17 / data bth1cs( 17) / +.2553521337 7973900223 1990525335 22 d-18 / data bth1cs( 18) / -.5962149019 7413450395 7682879078 49 d-19 / data bth1cs( 19) / +.1455662290 2372718620 2883020058 33 d-19 / data bth1cs( 20) / -.3702218542 2450538201 5797760195 93 d-20 / data bth1cs( 21) / +.9776307412 5345357664 1684345179 24 d-21 / data bth1cs( 22) / -.2672682163 9668488468 7237753930 52 d-21 / data bth1cs( 23) / +.7545330038 4983271794 0381906557 64 d-22 / data bth1cs( 24) / -.2194789991 9802744897 8923833716 47 d-22 / data bth1cs( 25) / +.6564839462 3955262178 9069998174 93 d-23 / data bth1cs( 26) / -.2015560429 8370207570 7840768695 19 d-23 / data bth1cs( 27) / +.6341776855 6776143492 1446671856 70 d-24 / data bth1cs( 28) / -.2041927788 5337895634 8137699555 91 d-24 / data bth1cs( 29) / +.6719146422 0720567486 6589800185 51 d-25 / data bth1cs( 30) / -.2256907911 0207573595 7090036873 36 d-25 / data bth1cs( 31) / +.7729771989 2989706370 9269598719 29 d-26 / data bth1cs( 32) / -.2696744451 2294640913 2114240809 20 d-26 / data bth1cs( 33) / +.9574934451 8502698072 2955219336 27 d-27 / data bth1cs( 34) / -.3456916844 8890113000 1756808276 27 d-27 / data bth1cs( 35) / +.1268123481 7398436504 2119862383 74 d-27 / data bth1cs( 36) / -.4723253663 0722639860 4649937134 45 d-28 / data bth1cs( 37) / +.1785000847 8186376177 8586197964 17 d-28 / data bth1cs( 38) / -.6840436100 4510395406 2152235667 46 d-29 / data bth1cs( 39) / +.2656602867 1720419358 2934226722 12 d-29 / data bth1cs( 40) / -.1045040252 7914452917 7141614846 70 d-29 / data bth1cs( 41) / +.4161829082 5377144306 8619171970 64 d-30 / data bth1cs( 42) / -.1677163920 3643714856 5013478828 87 d-30 / data bth1cs( 43) / +.6836199777 6664389173 5359280285 28 d-31 / data bth1cs( 44) / -.2817224786 1233641166 7395746228 10 d-31 / c data pi4 / 0.7853981633 9744830961 5660845819 876 d0 / data nbm1, nbt12, nbm12, nbth1, xmax / 4*0, 0.d0 / c if (nbm1.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) nbm1 = initds (bm1cs, 37, eta) nbt12 = initds (bt12cs, 39, eta) nbm12 = initds (bm12cs, 40, eta) nbth1 = initds (bth1cs, 44, eta) c xmax = 1.0d0/d1mach(4) c 10 if (x.lt.4.d0) call seteru (22hd9b1mp x must be ge 4, 22, 1, 2) c if (x.gt.8.d0) go to 20 z = (128.d0/(x*x) - 5.d0)/3.d0 ampl = (0.75d0 + dcsevl (z, bm1cs, nbm1))/dsqrt(x) theta = x - 3.d0*pi4 + dcsevl (z, bt12cs, nbt12)/x return c 20 if (x.gt.xmax) call seteru ( 1 37hd9b1mp no precision because x is big, 37, 2, 2) c z = 128.d0/(x*x) - 1.d0 ampl = (0.75d0 + dcsevl (z, bm12cs, nbm12))/dsqrt(x) theta = x - 3.d0*pi4 + dcsevl (z, bth1cs, nbth1)/x return c end double precision function dbsk0e (x) c july 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bk0cs(16), ak0cs(38), ak02cs(33), 1 xsml, y, d1mach, dcsevl, dbesi0, dexp, dlog, dsqrt c external d1mach, dbesi0, dcsevl, dexp, dlog, dsqrt, initds c c series for bk0 on the interval 0. to 4.00000e+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 series for ak0 on the interval 1.25000e-01 to 5.00000e-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 series for ak02 on the interval 0. to 1.25000e-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 / data ntk0, ntak0, ntak02, xsml / 3*0, 0.0d0 / c if (ntk0.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) ntk0 = initds (bk0cs, 16, eta) ntak0 = initds (ak0cs, 38, eta) ntak02 = initds (ak02cs, 33, eta) xsml = dsqrt (4.0d0*d1mach(3)) c 10 if (x.le.0.d0) call seteru (29hdbsk0e x is zero or negative, 29, 1 1, 2) if (x.gt.2.0d0) go to 20 c y = 0.d0 if (x.gt.xsml) y = x*x dbsk0e = dexp(x)*(-dlog(0.5d0*x)*dbesi0(x) - 0.25d0 + 1 dcsevl (.5d0*y-1.d0, bk0cs, ntk0)) return c 20 if (x.le.8.d0) dbsk0e = (1.25d0 + dcsevl ((16.d0/x-5.d0)/3.d0, 1 ak0cs, ntak0))/dsqrt(x) if (x.gt.8.d0) dbsk0e = (1.25d0 + 1 dcsevl (16.d0/x-1.d0, ak02cs, ntak02))/dsqrt(x) c return end double precision function dbsk1e (x) c july 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision x, bk1cs(16), ak1cs(38), ak12cs(33), xmin, 1 xsml, y, d1mach, dcsevl, dbesi1, dexp, dlog, dsqrt c external d1mach, dbesi1, dcsevl, dexp, dlog, dsqrt, initds c c series for bk1 on the interval 0. to 4.00000e+00 c with weighted error 9.16e-32 c log weighted error 31.04 c significant figures required 30.61 c decimal places required 31.64 c data bk1 cs( 1) / +.2530022733 8947770532 5311208685 33 d-1 / data bk1 cs( 2) / -.3531559607 7654487566 7238316918 01 d+0 / data bk1 cs( 3) / -.1226111808 2265714823 4790679300 42 d+0 / data bk1 cs( 4) / -.6975723859 6398643501 8129202960 83 d-2 / data bk1 cs( 5) / -.1730288957 5130520630 1765073689 79 d-3 / data bk1 cs( 6) / -.2433406141 5659682349 6007350301 64 d-5 / data bk1 cs( 7) / -.2213387630 7347258558 3152525451 26 d-7 / data bk1 cs( 8) / -.1411488392 6335277610 9583302126 08 d-9 / data bk1 cs( 9) / -.6666901694 1993290060 8537512643 73 d-12 / data bk1 cs( 10) / -.2427449850 5193659339 2631968648 53 d-14 / data bk1 cs( 11) / -.7023863479 3862875971 7837971200 00 d-17 / data bk1 cs( 12) / -.1654327515 5100994675 4910293333 33 d-19 / data bk1 cs( 13) / -.3233834745 9944491991 8933333333 33 d-22 / data bk1 cs( 14) / -.5331275052 9265274999 4666666666 66 d-25 / data bk1 cs( 15) / -.7513040716 2157226666 6666666666 66 d-28 / data bk1 cs( 16) / -.9155085717 6541866666 6666666666 66 d-31 / c c series for ak1 on the interval 1.25000e-01 to 5.00000e-01 c with weighted error 3.07e-32 c log weighted error 31.51 c significant figures required 30.71 c decimal places required 32.30 c data ak1 cs( 1) / +.2744313406 9738829695 2576662272 66 d+0 / data ak1 cs( 2) / +.7571989953 1993678170 8923781492 90 d-1 / data ak1 cs( 3) / -.1441051556 4754061229 8531161756 25 d-2 / data ak1 cs( 4) / +.6650116955 1257479394 2513854770 36 d-4 / data ak1 cs( 5) / -.4369984709 5201407660 5808450891 67 d-5 / data ak1 cs( 6) / +.3540277499 7630526799 4171390085 34 d-6 / data ak1 cs( 7) / -.3311163779 2932920208 9826882457 04 d-7 / data ak1 cs( 8) / +.3445977581 9010534532 3114997709 92 d-8 / data ak1 cs( 9) / -.3898932347 4754271048 9819374927 58 d-9 / data ak1 cs( 10) / +.4720819750 4658356400 9474493390 05 d-10 / data ak1 cs( 11) / -.6047835662 8753562345 3735915628 90 d-11 / data ak1 cs( 12) / +.8128494874 8658747888 1938379856 63 d-12 / data ak1 cs( 13) / -.1138694574 7147891428 9239159510 42 d-12 / data ak1 cs( 14) / +.1654035840 8462282325 9729482050 90 d-13 / data ak1 cs( 15) / -.2480902567 7068848221 5160104405 33 d-14 / data ak1 cs( 16) / +.3829237890 7024096948 4292272991 57 d-15 / data ak1 cs( 17) / -.6064734104 0012418187 7682103773 86 d-16 / data ak1 cs( 18) / +.9832425623 2648616038 1940046506 66 d-17 / data ak1 cs( 19) / -.1628416873 8284380035 6666201156 26 d-17 / data ak1 cs( 20) / +.2750153649 6752623718 2841203370 66 d-18 / data ak1 cs( 21) / -.4728966646 3953250924 2810695680 00 d-19 / data ak1 cs( 22) / +.8268150002 8109932722 3920503466 66 d-20 / data ak1 cs( 23) / -.1468140513 6624956337 1939648853 33 d-20 / data ak1 cs( 24) / +.2644763926 9208245978 0858948266 66 d-21 / data ak1 cs( 25) / -.4829015756 4856387897 9698688000 00 d-22 / data ak1 cs( 26) / +.8929302074 3610130180 6563327999 99 d-23 / data ak1 cs( 27) / -.1670839716 8972517176 9977514666 66 d-23 / data ak1 cs( 28) / +.3161645603 4040694931 3686186666 66 d-24 / data ak1 cs( 29) / -.6046205531 2274989106 5064106666 66 d-25 / data ak1 cs( 30) / +.1167879894 2042732700 7184213333 33 d-25 / data ak1 cs( 31) / -.2277374158 2653996232 8678400000 00 d-26 / data ak1 cs( 32) / +.4481109730 0773675795 3058133333 33 d-27 / data ak1 cs( 33) / -.8893288476 9020194062 3360000000 00 d-28 / data ak1 cs( 34) / +.1779468001 8850275131 3920000000 00 d-28 / data ak1 cs( 35) / -.3588455596 7329095821 9946666666 66 d-29 / data ak1 cs( 36) / +.7290629049 2694257991 6799999999 99 d-30 / data ak1 cs( 37) / -.1491844984 5546227073 0240000000 00 d-30 / data ak1 cs( 38) / +.3073657387 2934276300 7999999999 99 d-31 / c c series for ak12 on the interval 0. to 1.25000e-01 c with weighted error 2.41e-32 c log weighted error 31.62 c significant figures required 30.25 c decimal places required 32.38 c data ak12cs( 1) / +.6379308343 7390010366 0048853410 2 d-1 / data ak12cs( 2) / +.2832887813 0497209358 3503028470 8 d-1 / data ak12cs( 3) / -.2475370673 9052503454 1454556673 2 d-3 / data ak12cs( 4) / +.5771972451 6072488204 7097662576 3 d-5 / data ak12cs( 5) / -.2068939219 5365483027 4553319655 2 d-6 / data ak12cs( 6) / +.9739983441 3818041803 0921309788 7 d-8 / data ak12cs( 7) / -.5585336140 3806249846 8889551112 9 d-9 / data ak12cs( 8) / +.3732996634 0461852402 2121285473 1 d-10 / data ak12cs( 9) / -.2825051961 0232254451 3506575492 8 d-11 / data ak12cs( 10) / +.2372019002 4841441736 4349695548 6 d-12 / data ak12cs( 11) / -.2176677387 9917539792 6830166793 8 d-13 / data ak12cs( 12) / +.2157914161 6160324539 3956268970 6 d-14 / data ak12cs( 13) / -.2290196930 7182692759 9155133815 4 d-15 / data ak12cs( 14) / +.2582885729 8232749619 1993956522 6 d-16 / data ak12cs( 15) / -.3076752641 2684631876 2109817344 0 d-17 / data ak12cs( 16) / +.3851487721 2804915970 9489684479 9 d-18 / data ak12cs( 17) / -.5044794897 6415289771 1728250880 0 d-19 / data ak12cs( 18) / +.6888673850 4185442370 1829222399 9 d-20 / data ak12cs( 19) / -.9775041541 9501183030 0213248000 0 d-21 / data ak12cs( 20) / +.1437416218 5238364610 0165973333 3 d-21 / data ak12cs( 21) / -.2185059497 3443473734 9973333333 3 d-22 / data ak12cs( 22) / +.3426245621 8092206316 4538880000 0 d-23 / data ak12cs( 23) / -.5531064394 2464082325 0124800000 0 d-24 / data ak12cs( 24) / +.9176601505 6859954037 8282666666 6 d-25 / data ak12cs( 25) / -.1562287203 6180249114 4874666666 6 d-25 / data ak12cs( 26) / +.2725419375 4843331323 4943999999 9 d-26 / data ak12cs( 27) / -.4865674910 0748279923 7802666666 6 d-27 / data ak12cs( 28) / +.8879388552 7235025873 5786666666 6 d-28 / data ak12cs( 29) / -.1654585918 0392575489 3653333333 3 d-28 / data ak12cs( 30) / +.3145111321 3578486743 0399999999 9 d-29 / data ak12cs( 31) / -.6092998312 1931276124 1600000000 0 d-30 / data ak12cs( 32) / +.1202021939 3698158346 2399999999 9 d-30 / data ak12cs( 33) / -.2412930801 4594088413 8666666666 6 d-31 / c data ntk1, ntak1, ntak12, xmin, xsml / 3*0, 2*0.d0 / c if (ntk1.ne.0) go to 10 eta = 0.1*sngl(d1mach(3)) ntk1 = initds (bk1cs, 16, eta) ntak1 = initds (ak1cs, 38, eta) ntak12 = initds (ak12cs, 33, eta) c xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0) xsml = dsqrt (4.0d0*d1mach(3)) c 10 if (x.le.0.d0) call seteru (29hdbsk1e x is zero or negative, 29, 1 1, 2) if (x.gt.2.0d0) go to 20 c if (x.lt.xmin) call seteru (31hdbsk1e x so small k1 overflows, 1 31, 2, 2) y = 0.d0 if (x.gt.xsml) y = x*x dbsk1e = dexp(x)*(dlog(0.5d0*x)*dbesi1(x) + (0.75d0 + 1 dcsevl (0.5d0*y-1.d0, bk1cs, ntk1))/x ) return c 20 if (x.le.8.d0) dbsk1e = (1.25d0 + dcsevl ((16.d0/x-5.d0)/3.d0, 1 ak1cs, ntak1))/dsqrt(x) if (x.gt.8.d0) dbsk1e = (1.25d0 + 1 dcsevl (16.d0/x-1.d0, ak12cs, ntak12))/dsqrt(x) c return end subroutine seterr (messg, nmessg, nerr, iopt) c c this version modified by w. fullerton to dump if iopt = 1 and c not recovering. c seterr sets lerror = nerr, optionally prints the message and dumps c according to the following rules... c c if iopt = 1 and recovering - just remember the error. c if iopt = 1 and not recovering - print, dump and stop. c if iopt = 2 - print, dump and stop. c c input c c messg - the error message. c nmessg - the length of the message, in characters. c nerr - the error number. must have nerr non-zero. c iopt - the option. must have iopt=1 or 2. c c error states - c c 1 - message length not positive. c 2 - cannot have nerr=0. c 3 - an unrecovered error followed by another error. c 4 - bad value for iopt. c c only the first 72 characters of the message are printed. c c the error handler calls a subroutine named fdump to produce a c symbolic dump. to complete the package, a dummy version of fdump c is supplied, but it should be replaced by a locally written version c which at least gives a trace-back. c integer messg(1) c external i1mach, i8save c c the unit for error messages. c iwunit=i1mach(4) c if (nmessg.ge.1) go to 10 c c a message of non-positive length is fatal. c write(iwunit,9000) 9000 format(52h1error 1 in seterr - message length not positive.) go to 60 c c nw is the number of words the message occupies. c 10 nw=(min0(nmessg,72)-1)/i1mach(6)+1 c if (nerr.ne.0) go to 20 c c cannot turn the error state off using seterr. c write(iwunit,9001) 9001 format(42h1error 2 in seterr - cannot have nerr=0// 1 34h the current error message follows///) call e9rint(messg,nw,nerr,.true.) itemp=i8save(1,1,.true.) go to 50 c c set lerror and test for a previous unrecovered error. c 20 if (i8save(1,nerr,.true.).eq.0) go to 30 c write(iwunit,9002) 9002 format(23h1error 3 in seterr -, 1 48h an unrecovered error followed by another error.// 2 48h the previous and current error messages follow.///) call eprint call e9rint(messg,nw,nerr,.true.) go to 50 c c save this message in case it is not recovered from properly. c 30 call e9rint(messg,nw,nerr,.true.) c if (iopt.eq.1 .or. iopt.eq.2) go to 40 c c must have iopt = 1 or 2. c write(iwunit,9003) 9003 format(42h1error 4 in seterr - bad value for iopt// 1 34h the current error message follows///) go to 50 c c test for recovery. c 40 if (iopt.eq.2) go to 50 c if (i8save(2,0,.false.).eq.1) return c c call eprint c stop c 50 call eprint 60 call fdump stop c end subroutine e9rint(messg,nw,nerr,save) c c this routine stores the current error message or prints the old one, c if any, depending on whether or not save = .true. . c integer messg(nw) logical save c external i1mach, i8save c c messgp stores at least the first 72 characters of the previous c message. its length is machine dependent and must be at least c c 1 + 71/(the number of characters stored per integer word). c integer messgp(36),fmt(14),ccplus c c start with no previous message. c data messgp(1)/1h1/, nwp/0/, nerrp/0/ c c set up the format for printing the error message. c the format is simply (a1,14x,72axx) where xx=i1mach(6) is the c number of characters stored per integer word. c data ccplus / 1h+ / c data fmt( 1) / 1h( / data fmt( 2) / 1ha / data fmt( 3) / 1h1 / data fmt( 4) / 1h, / data fmt( 5) / 1h1 / data fmt( 6) / 1h4 / data fmt( 7) / 1hx / data fmt( 8) / 1h, / data fmt( 9) / 1h7 / data fmt(10) / 1h2 / data fmt(11) / 1ha / data fmt(12) / 1hx / data fmt(13) / 1hx / data fmt(14) / 1h) / c if (.not.save) go to 20 c c save the message. c nwp=nw nerrp=nerr do 10 i=1,nw 10 messgp(i)=messg(i) c go to 30 c 20 if (i8save(1,0,.false.).eq.0) go to 30 c c print the message. c iwunit=i1mach(4) write(iwunit,9000) nerrp 9000 format(7h error ,i4,4h in ) c call s88fmt(2,i1mach(6),fmt(12)) write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp) c 30 return c end integer function i8save(isw,ivalue,set) c c if (isw = 1) i8save returns the current error number and c sets it to ivalue if set = .true. . c c if (isw = 2) i8save returns the current recovery switch and c sets it to ivalue if set = .true. . c logical set c integer iparam(2) c iparam(1) is the error number and iparam(2) is the recovery switch. c c start execution error free and with recovery turned off. c data iparam(1) /0/, iparam(2) /2/ c i8save=iparam(isw) if (set) iparam(isw)=ivalue c return c end subroutine eprint c c this subroutine prints the last error message, if any. c integer messg(1) c call e9rint(messg,1,1,.false.) return c end subroutine s88fmt( n, w, ifmt ) c c s88fmt replaces ifmt(1), ... , ifmt(n) with c the characters corresponding to the n least significant c digits of w. c integer n,w,ifmt(n) c integer nt,wt,digits(10) c data digits( 1) / 1h0 / data digits( 2) / 1h1 / data digits( 3) / 1h2 / data digits( 4) / 1h3 / data digits( 5) / 1h4 / data digits( 6) / 1h5 / data digits( 7) / 1h6 / data digits( 8) / 1h7 / data digits( 9) / 1h8 / data digits(10) / 1h9 / c nt = n wt = w c 10 if (nt .le. 0) return idigit = mod( wt, 10 ) ifmt(nt) = digits(idigit+1) wt = wt/10 nt = nt - 1 go to 10 c end subroutine fdump return end REAL FUNCTION SECOND () C C OBTAIN ELAPSED CPU TIME SINCE START OF RUN. C C ON CRAY -- THIS ROUTINE NOT NEEDED. C C ON CONVEX -- MORE ACCURATE TIMINGS ARE OBTAINED BY CALLING C THE VECLIB ROUTINE CPUTIME AS FOLLOWS C C SECOND = CPUTIME(0.0E0) C C ON UNIX -- (SUN, FOR EXAMPLE) THE FOLLOWING MAY BE USED C REAL TARRAY(2), ETIME, TOTAL TOTAL = ETIME(TARRAY) SECOND = TARRAY(1) + TARRAY(2) C RETURN END