Blame view
fvn_fnlib/si.f
1.62 KB
38581db0c
git-svn-id: https...
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
function si (x)
c december 1980 edition, w. fullerton, bell labs.
dimension sics(12)
external csevl, inits, r1mach
c
c series for si on the interval 0.00000e+00 to 1.60000e+01
c with weighted error 1.22e-17
c log weighted error 16.91
c significant figures required 16.37
c decimal places required 17.45
c
data si cs( 1) / -0.1315646598 184841929 e0/
data si cs( 2) / -0.2776578526 973601892 e0/
data si cs( 3) / 0.0354414054 866659180 e0/
data si cs( 4) / -0.0025631631 447933978 e0/
data si cs( 5) / 0.0001162365 390497009 e0/
data si cs( 6) / -0.0000035904 327241606 e0/
data si cs( 7) / 0.0000000802 342123706 e0/
data si cs( 8) / -0.0000000013 562997693 e0/
data si cs( 9) / 0.0000000000 179440722 e0/
data si cs( 10) / -0.0000000000 001908387 e0/
data si cs( 11) / 0.0000000000 000016670 e0/
data si cs( 12) / -0.0000000000 000000122 e0/
c
data pi2 / 1.570796326 7948966e0 /
data nsi, xsml /0, 0.0/
c
if (nsi.ne.0) go to 10
nsi = inits (sics, 12, 0.1*r1mach(3))
xsml = sqrt (eps)
c
10 absx = abs(x)
if (absx.gt.4.0) go to 20
si = x
if (absx.lt.xsml) return
c
si = x*(0.75 + csevl ((x*x-8.0)*0.125, sics, nsi))
return
c
20 call r9sifg (absx, f, g)
cosx = cos (absx)
call erroff
si = pi2 - f*cosx - g*sin(x)
if (x.lt.0.0) si = -si
c
return
end
|