Blame view

fvn_fnlib/dsi.f 2.32 KB
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
        double precision function dsi (x)
  c december 1980 edition, w. fullerton, bell labs.
        double precision x, sics(18), pi2, xsml, absx, f, g, cosx,
       1  dcsevl, d1mach, dcos, dsin, dsqrt
        external d1mach,  dcsevl, initds
  c
  c series for si   on the interval  0.00000e+00 to  1.60000e+01
  c                                        with weighted error   8.58e-32
  c                                         log weighted error  31.07
  c                               significant figures required  30.53
  c                                    decimal places required  31.69
  c
        data si  cs(  1) / -0.1315646598 1848419289 0427517300 0457d0/
        data si  cs(  2) / -0.2776578526 9736018920 4828766015 7299d0/
        data si  cs(  3) /  0.0354414054 8666591797 4913546471 0086d0/
        data si  cs(  4) / -0.0025631631 4479339776 5875278836 1530d0/
        data si  cs(  5) /  0.0001162365 3904970092 8126492148 2985d0/
        data si  cs(  6) / -0.0000035904 3272416060 4267000434 7148d0/
        data si  cs(  7) /  0.0000000802 3421237057 1016230865 2976d0/
        data si  cs(  8) / -0.0000000013 5629976925 4025064993 1846d0/
        data si  cs(  9) /  0.0000000000 1794407215 9973677556 7759d0/
        data si  cs( 10) / -0.0000000000 0019083873 4308714549 0737d0/
        data si  cs( 11) /  0.0000000000 0000166699 8958682433 0853d0/
        data si  cs( 12) / -0.0000000000 0000001217 3098836850 3042d0/
        data si  cs( 13) /  0.0000000000 0000000007 5418186699 3865d0/
        data si  cs( 14) / -0.0000000000 0000000000 0401417884 2446d0/
        data si  cs( 15) /  0.0000000000 0000000000 0001855369 0716d0/
        data si  cs( 16) / -0.0000000000 0000000000 0000007516 6966d0/
        data si  cs( 17) /  0.0000000000 0000000000 0000000026 9113d0/
        data si  cs( 18) / -0.0000000000 0000000000 0000000000 0858d0/
  c
        data pi2 / 1.5707963267 9489661923 1321691639 75 d0 /
        data nsi, xsml /0, 0.0d0/
  c
        if (nsi.ne.0) go to 10
        nsi = initds (sics, 18, 0.1*sngl(d1mach(3)))
        xsml = dsqrt(d1mach(3))
  c
   10   absx = dabs(x)
        if (absx.gt.4.0d0) go to 20
        dsi = x
        if (absx.lt.xsml) return
  c
        dsi = x*(0.75d0 + dcsevl ((x*x-8.0d0)*.125d0, sics, nsi))
        return
  c
   20   call d9sifg (absx, f, g)
        cosx = dcos (absx)
        call erroff
        dsi = pi2 - f*cosx - g*dsin(x)
        if (x.lt.0.0d0) dsi = -dsi
  c
        return
        end