1 DOUBLE PRECISION FUNCTION DBESY0 (X) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1996-03-30 DBESY0 Krogh Added external statement. 6C>> 1995-11-13 DBESY0 Krogh Changes to simplify conversion to C. 7C>> 1994-11-11 DBESY0 Krogh Declared all vars. 8C>> 1994-10-20 DBESY0 Krogh Changes to use M77CON 9C>> 1990-11-29 DBESY0 CLL 10C>> 1985-08-02 DBESY0 Lawson Initial code. 11C JULY 1977 EDITION. W. FULLERTON, C3, LOS ALAMOS SCIENTIFIC LAB. 12C C.L.LAWSON & S.CHAN, JPL, 1984 FEB ADAPTED TO JPL MATH77 LIBRARY. 13c ------------------------------------------------------------------ 14c--D replaces "?": ?BESY0, ?BESJ0, ?BMP0, ?INITS, ?CSEVL, ?ERM1 15c ------------------------------------------------------------------ 16 EXTERNAL D1MACH, DBESJ0, DCSEVL 17 INTEGER NTY0 18 DOUBLE PRECISION X, BY0CS(19), AMPL, THETA, TWODPI, XSML, 19 1 Y, D1MACH, DCSEVL, DBESJ0 20C 21C SERIES FOR BY0 ON THE INTERVAL 0. TO 1.60000D+01 22C WITH WEIGHTED ERROR 8.14D-32 23C LOG WEIGHTED ERROR 31.09 24C SIGNIFICANT FIGURES REQUIRED 30.31 25C DECIMAL PLACES REQUIRED 31.73 26C 27 SAVE NTY0, XSML 28C 29 DATA BY0CS / -.1127783939286557321793980546028D-1, 30 * -.1283452375604203460480884531838D+0, 31 * -.1043788479979424936581762276618D+0, 32 * +.2366274918396969540924159264613D-1, 33 * -.2090391647700486239196223950342D-2, 34 * +.1039754539390572520999246576381D-3, 35 * -.3369747162423972096718775345037D-5, 36 * +.7729384267670667158521367216371D-7, 37 * -.1324976772664259591443476068964D-8, 38 * +.1764823261540452792100389363158D-10, 39 * -.1881055071580196200602823012069D-12, 40 * +.1641865485366149502792237185749D-14, 41 * -.1195659438604606085745991006720D-16, 42 * +.7377296297440185842494112426666D-19, 43 * -.3906843476710437330740906666666D-21, 44 * +.1795503664436157949829120000000D-23, 45 * -.7229627125448010478933333333333D-26, 46 * +.2571727931635168597333333333333D-28, 47 * -.8141268814163694933333333333333D-31 / 48C 49 DATA TWODPI / 0.636619772367581343075535053490057D0 / 50 DATA NTY0, XSML / 0, 0.D0 / 51C ------------------------------------------------------------------ 52 IF (NTY0 .eq. 0) then 53 call DINITS (BY0CS, 19, 0.1D0*D1MACH(3), NTY0) 54 XSML = SQRT (4.0D0*D1MACH(3)) 55 endif 56C 57 IF (X .LE. 0.D0) THEN 58 DBESY0 = 0.D0 59 CALL DERM1 ('DBESY0',1,0,'X IS ZERO OR NEGATIVE','X',X,'.') 60 ELSE IF (X .LE. 4.D0) THEN 61 IF (X .LE. XSML) THEN 62 Y = 0.D0 63 ELSE 64 Y = X * X 65 END IF 66 DBESY0 = TWODPI*LOG(0.5D0*X)*DBESJ0(X) + .375D0 + 67 * DCSEVL (.125D0*Y-1.D0, BY0CS, NTY0) 68 ELSE 69 CALL DBMP0 (X, AMPL, THETA) 70 DBESY0 = AMPL * SIN(THETA) 71 END IF 72C 73 RETURN 74C 75 END 76