Blame view

fvn_fnlib/dai.f 3.66 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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
        double precision function dai (x)
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
        double precision x, aifcs(13), aigcs(13), theta, xm, xmax, x3sml,
       1  z, d1mach, dcsevl, daie, dcos, dexp, dlog, dsqrt
        external d1mach, daie, dcsevl, initds
  c
  c series for aif        on the interval -1.00000e+00 to  1.00000e+00
  c                                        with weighted error   8.37e-33
  c                                         log weighted error  32.08
  c                               significant figures required  30.87
  c                                    decimal places required  32.63
  c
        data aif cs(  1) / -.3797135849 6669997496 1970894694 14 d-1     /
        data aif cs(  2) / +.5919188853 7263638574 3197280137 77 d-1     /
        data aif cs(  3) / +.9862928057 7279975365 6038910440 60 d-3     /
        data aif cs(  4) / +.6848843819 0765667554 8548301824 12 d-5     /
        data aif cs(  5) / +.2594202596 2194713019 4892790814 03 d-7     /
        data aif cs(  6) / +.6176612774 0813750329 4457496972 36 d-10    /
        data aif cs(  7) / +.1009245417 2466117901 4295562246 01 d-12    /
        data aif cs(  8) / +.1201479251 1179938141 2880332253 33 d-15    /
        data aif cs(  9) / +.1088294558 8716991878 5252954666 66 d-18    /
        data aif cs( 10) / +.7751377219 6684887039 2384000000 00 d-22    /
        data aif cs( 11) / +.4454811203 7175638391 4666666666 66 d-25    /
        data aif cs( 12) / +.2109284523 1692343466 6666666666 66 d-28    /
        data aif cs( 13) / +.8370173591 0741333333 3333333333 33 d-32    /
  c
  c series for aig        on the interval -1.00000e+00 to  1.00000e+00
  c                                        with weighted error   7.47e-34
  c                                         log weighted error  33.13
  c                               significant figures required  31.50
  c                                    decimal places required  33.68
  c
        data aig cs(  1) / +.1815236558 1161273011 5562099578 64 d-1     /
        data aig cs(  2) / +.2157256316 6010755534 0306388199 68 d-1     /
        data aig cs(  3) / +.2567835698 7483249659 0524280901 33 d-3     /
        data aig cs(  4) / +.1426521411 9792403898 8294969217 21 d-5     /
        data aig cs(  5) / +.4572114920 0180426070 4340975581 91 d-8     /
        data aig cs(  6) / +.9525170843 5647098607 3922788405 92 d-11    /
        data aig cs(  7) / +.1392563460 5771399051 1504206861 90 d-13    /
        data aig cs(  8) / +.1507099914 2762379592 3069911386 66 d-16    /
        data aig cs(  9) / +.1255914831 2567778822 7032053333 33 d-19    /
        data aig cs( 10) / +.8306307377 0821340343 8293333333 33 d-23    /
        data aig cs( 11) / +.4465753849 3718567445 3333333333 33 d-26    /
        data aig cs( 12) / +.1990085503 4518869333 3333333333 33 d-29    /
        data aig cs( 13) / +.7470288525 6533333333 3333333333 33 d-33    /
  c
        data naif, naig, x3sml, xmax / 2*0, 2*0.d0 /
  c
        if (naif.ne.0) go to 10
        naif = initds (aifcs, 13, 0.1*sngl(d1mach(3)))
        naig = initds (aigcs, 13, 0.1*sngl(d1mach(3)))
  c
        x3sml = d1mach(3)**0.3334d0
        xmax = (-1.5d0*dlog(d1mach(1)))**0.6667d0
        xmax = xmax - xmax*dlog(xmax)/(4.0d0*xmax*dsqrt(xmax)+1.0d0)
       1  - 0.01d0
  c
   10   if (x.ge.(-1.d0)) go to 20
        call d9aimp (x, xm, theta)
        dai = xm * dcos(theta)
        return
  c
   20   if (x.gt.1.0d0) go to 30
        z = 0.0d0
        if (dabs(x).gt.x3sml) z = x**3
        dai = 0.375d0 + (dcsevl (z, aifcs, naif) - x*(0.25d0 +
       1  dcsevl (z, aigcs, naig)) )
        return
  c
   30   if (x.gt.xmax) go to 40
        dai = daie(x) * dexp(-2.0d0*x*dsqrt(x)/3.0d0)
        return
  c
   40   dai = 0.0d0
        call seteru (30hdai     x so big ai underflows, 30, 1, 0)
        return
  c
        end