si.f
1.62 KB
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