lxsd / fvn

Blame view

fvn_fnlib/r9atn1.f 2.39 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
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
55
56
57
58
59
60
61
62
63
64
65
66
67
 
        function r9atn1 (x)
  c april 1978 edition.  w. fullerton c3, los alamos scientific lab.
  c
  c evaluate  atan(x)  from first order, that is, evaluate
  c (atan(x)-x)/x**3  with relative error accuracy so that
  c        atan(x) = x + x**3*r9atn1(x).
  c
        dimension atn1cs(21)
        external  csevl, inits, r1mach
  c
  c series for atn1       on the interval  0.          to  1.00000d+00
  c                                        with weighted error   2.21e-17
  c                                         log weighted error  16.66
  c                               significant figures required  15.44
  c                                    decimal places required  17.32
  c
        data atn1cs( 1) /   -.0328399753 5355202e0 /
        data atn1cs( 2) /    .0583343234 3172412e0 /
        data atn1cs( 3) /   -.0074003696 9671964e0 /
        data atn1cs( 4) /    .0010097841 9933728e0 /
        data atn1cs( 5) /   -.0001439787 1635652e0 /
        data atn1cs( 6) /    .0000211451 2648992e0 /
        data atn1cs( 7) /   -.0000031723 2107425e0 /
        data atn1cs( 8) /    .0000004836 6203654e0 /
        data atn1cs( 9) /   -.0000000746 7746546e0 /
        data atn1cs(10) /    .0000000116 4800896e0 /
        data atn1cs(11) /   -.0000000018 3208837e0 /
        data atn1cs(12) /    .0000000002 9019082e0 /
        data atn1cs(13) /   -.0000000000 4623885e0 /
        data atn1cs(14) /    .0000000000 0740552e0 /
        data atn1cs(15) /   -.0000000000 0119135e0 /
        data atn1cs(16) /    .0000000000 0019240e0 /
        data atn1cs(17) /   -.0000000000 0003118e0 /
        data atn1cs(18) /    .0000000000 0000506e0 /
        data atn1cs(19) /   -.0000000000 0000082e0 /
        data atn1cs(20) /    .0000000000 0000013e0 /
        data atn1cs(21) /   -.0000000000 0000002e0 /
  c
        data ntatn1, xsml, xbig, xmax / 0, 3*0.0 /
  c
        if (ntatn1.ne.0) go to 10
        eps = r1mach(3)
        ntatn1 = inits (atn1cs, 21, 0.1*eps)
  c
        xsml = sqrt (0.1*eps)
        xbig = 1.571/sqrt(eps)
        xmax = 1.571/eps
  c
   10   y = abs(x)
        if (y.gt.1.0) go to 20
  c
        if (y.le.xsml) r9atn1 = -1.0/3.0
        if (y.le.xsml) return
  c
        r9atn1 = -0.25 + csevl (2.0*y*y-1., atn1cs, ntatn1)
        return
  c
   20   if (y.gt.xmax) call seteru (
       1  51hr9atn1  no precision in answer because x is too big, 51,2,2)
        if (y.gt.xbig) call seteru (
       1  53hr9atn1  answer lt half precision because x is too big, 53,
       2  1, 1)
  c
        r9atn1 = (atan(x) - x) / x**3
        return
  c
        end