lxsd / fvn

Blame view

fvn_fnlib/si.f 1.62 KB
edit raw normal view history
38581db0c   daniau   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
 
        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))
1258bdafa   cwaterkeyn   ChW _ add fvnlib_... 
30
 
        xsml = sqrt (r1mach(3))
38581db0c   daniau   git-svn-id: https... 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 
  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