Blame view
fvn_fnlib/aide.f
5.9 KB
38581db0c
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 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 |
function aide (x)
c july 1980 edition. w. fullerton, bell labs.
c
c evaluate the derivative of the airy function for x .le. 0
c and evaluate aid(x)*exp(zeta) for x .ge. 0 where
c zeta = 2/3 * x**(3/2)
c
dimension aifcs(8), aigcs(9), aip1cs(25), aip2cs(15)
external csevl, inits, r1mach
c
c series for aif on the interval -1.00000e+00 to 1.00000e+00
c with weighted error 5.22e-18
c log weighted error 17.28
c significant figures required 16.01
c decimal places required 17.73
c
data aifcs( 1) / 0.1052746122 6531408809 e0/
data aifcs( 2) / 0.0118361362 8152997844 e0/
data aifcs( 3) / 0.0001232810 4173225664 e0/
data aifcs( 4) / 0.0000006226 1225638140 e0/
data aifcs( 5) / 0.0000000018 5298887844 e0/
data aifcs( 6) / 0.0000000000 0363328873 e0/
data aifcs( 7) / 0.0000000000 0000504622 e0/
data aifcs( 8) / 0.0000000000 0000000522 e0/
c
c series for aig on the interval -1.00000e+00 to 1.00000e+00
c with weighted error 3.14e-19
c log weighted error 18.50
c significant figures required 17.44
c decimal places required 18.98
c
data aigcs( 1) / 0.0212338781 5091866685 2e0/
data aigcs( 2) / 0.0863159303 3521440675 2e0/
data aigcs( 3) / 0.0017975947 2038323135 8e0/
data aigcs( 4) / 0.0000142654 9987555069 3e0/
data aigcs( 5) / 0.0000000594 3799528368 3e0/
data aigcs( 6) / 0.0000000001 5240336647 9e0/
data aigcs( 7) / 0.0000000000 0026458766 0e0/
data aigcs( 8) / 0.0000000000 0000033156 2e0/
data aigcs( 9) / 0.0000000000 0000000031 4e0/
c
c series for aip2 on the interval 0.00000e+00 to 1.25000e-01
c with weighted error 2.15e-17
c log weighted error 16.67
c significant figures required 14.27
c decimal places required 17.26
c
data aip2cs( 1) / 0.0065457691 989713757 e0/
data aip2cs( 2) / 0.0023833724 120774592 e0/
data aip2cs( 3) / -0.0000430700 770220586 e0/
data aip2cs( 4) / 0.0000015629 125858629 e0/
data aip2cs( 5) / -0.0000000815 417186163 e0/
data aip2cs( 6) / 0.0000000054 103738057 e0/
data aip2cs( 7) / -0.0000000004 284130883 e0/
data aip2cs( 8) / 0.0000000000 389497963 e0/
data aip2cs( 9) / -0.0000000000 039623161 e0/
data aip2cs( 10) / 0.0000000000 004428184 e0/
data aip2cs( 11) / -0.0000000000 000536297 e0/
data aip2cs( 12) / 0.0000000000 000069650 e0/
data aip2cs( 13) / -0.0000000000 000009620 e0/
data aip2cs( 14) / 0.0000000000 000001403 e0/
data aip2cs( 15) / -0.0000000000 000000215 e0/
c
c series for aip1 on the interval 1.25000e-01 to 1.00000e+00
c with weighted error 2.60e-17
c log weighted error 16.58
c significant figures required 14.91
c decimal places required 17.28
c
data aip1cs( 1) / 0.0358865097 808301538 e0/
data aip1cs( 2) / 0.0114668575 627764899 e0/
data aip1cs( 3) / -0.0007592073 583861400 e0/
data aip1cs( 4) / 0.0000869517 610893841 e0/
data aip1cs( 5) / -0.0000128237 294298592 e0/
data aip1cs( 6) / 0.0000022062 695681038 e0/
data aip1cs( 7) / -0.0000004222 295185921 e0/
data aip1cs( 8) / 0.0000000874 686415726 e0/
data aip1cs( 9) / -0.0000000192 773588418 e0/
data aip1cs( 10) / 0.0000000044 668460054 e0/
data aip1cs( 11) / -0.0000000010 790108052 e0/
data aip1cs( 12) / 0.0000000002 700029447 e0/
data aip1cs( 13) / -0.0000000000 696480108 e0/
data aip1cs( 14) / 0.0000000000 184489907 e0/
data aip1cs( 15) / -0.0000000000 050027817 e0/
data aip1cs( 16) / 0.0000000000 013852243 e0/
data aip1cs( 17) / -0.0000000000 003908218 e0/
data aip1cs( 18) / 0.0000000000 001121536 e0/
data aip1cs( 19) / -0.0000000000 000326862 e0/
data aip1cs( 20) / 0.0000000000 000096619 e0/
data aip1cs( 21) / -0.0000000000 000028935 e0/
data aip1cs( 22) / 0.0000000000 000008770 e0/
data aip1cs( 23) / -0.0000000000 000002688 e0/
data aip1cs( 24) / 0.0000000000 000000832 e0/
data aip1cs( 25) / -0.0000000000 000000260 e0/
c
data naif, naig, naip1, naip2 / 4*0 /
data x2sml, x3sml, x32sml, xbig / 4*0.0 /
c
if (naif.ne.0) go to 10
eta = 0.1*r1mach(3)
naif = inits (aifcs, 8, eta)
naig = inits (aigcs, 9, eta)
naip1 = inits (aip1cs, 25, eta)
naip2 = inits (aip2cs, 15, eta)
c
x2sml = sqrt (eta)
x3sml = eta**0.3333
x32sml = 1.3104*x3sml**2
xbig = r1mach(2)**0.6666
c
10 if (x.ge.(-1.0)) go to 20
call r9admp (x, xn, phi)
aide = xn * cos(phi)
return
c
20 if (x.gt.1.0) go to 30
x3 = 0.0
if (abs(x).gt.x3sml) x3 = x**3
x2 = 0.0
if (abs(x).gt.x2sml) x2 = x*x
aide = (x2*(0.125 + csevl (x3, aifcs, naif)) -
1 csevl (x3, aigcs, naig)) - 0.25
if (x.gt.x32sml) aide = aide * exp (2.0*x*sqrt(x)/3.0)
return
c
30 if (x.gt.4.0) go to 40
sqrtx = sqrt(x)
z = (16.0/(x*sqrtx) - 9.0)/7.0
aide = (-0.28125 - csevl (z, aip1cs, naip1)) * sqrt(sqrtx)
return
c
40 sqrtx = sqrt(x)
z = -1.0
if (x.lt.xbig) z = 16.0/(x*sqrtx) - 1.0
aide = (-0.28125 - csevl (z, aip2cs, naip2)) * sqrt(sqrtx)
return
c
end
|