Blame view
fvn_fnlib/r9aimp.f
9.45 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 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 |
subroutine r9aimp (x, ampl, theta)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate the airy modulus and phase for x .le. -1.0
c
dimension am21cs(40), ath1cs(36), am22cs(33), ath2cs(32)
external csevl, inits, r1mach
c
c series for am21 on the interval -1.25000d-01 to 0.
c with weighted error 2.89e-17
c log weighted error 16.54
c significant figures required 14.15
c decimal places required 17.34
c
data am21cs( 1) / .0065809191 761485e0 /
data am21cs( 2) / .0023675984 685722e0 /
data am21cs( 3) / .0001324741 670371e0 /
data am21cs( 4) / .0000157600 904043e0 /
data am21cs( 5) / .0000027529 702663e0 /
data am21cs( 6) / .0000006102 679017e0 /
data am21cs( 7) / .0000001595 088468e0 /
data am21cs( 8) / .0000000471 033947e0 /
data am21cs( 9) / .0000000152 933871e0 /
data am21cs(10) / .0000000053 590722e0 /
data am21cs(11) / .0000000020 000910e0 /
data am21cs(12) / .0000000007 872292e0 /
data am21cs(13) / .0000000003 243103e0 /
data am21cs(14) / .0000000001 390106e0 /
data am21cs(15) / .0000000000 617011e0 /
data am21cs(16) / .0000000000 282491e0 /
data am21cs(17) / .0000000000 132979e0 /
data am21cs(18) / .0000000000 064188e0 /
data am21cs(19) / .0000000000 031697e0 /
data am21cs(20) / .0000000000 015981e0 /
data am21cs(21) / .0000000000 008213e0 /
data am21cs(22) / .0000000000 004296e0 /
data am21cs(23) / .0000000000 002284e0 /
data am21cs(24) / .0000000000 001232e0 /
data am21cs(25) / .0000000000 000675e0 /
data am21cs(26) / .0000000000 000374e0 /
data am21cs(27) / .0000000000 000210e0 /
data am21cs(28) / .0000000000 000119e0 /
data am21cs(29) / .0000000000 000068e0 /
data am21cs(30) / .0000000000 000039e0 /
data am21cs(31) / .0000000000 000023e0 /
data am21cs(32) / .0000000000 000013e0 /
data am21cs(33) / .0000000000 000008e0 /
data am21cs(34) / .0000000000 000005e0 /
data am21cs(35) / .0000000000 000003e0 /
data am21cs(36) / .0000000000 000001e0 /
data am21cs(37) / .0000000000 000001e0 /
data am21cs(38) / .0000000000 000000e0 /
data am21cs(39) / .0000000000 000000e0 /
data am21cs(40) / .0000000000 000000e0 /
c
c series for ath1 on the interval -1.25000d-01 to 0.
c with weighted error 2.53e-17
c log weighted error 16.60
c significant figures required 15.15
c decimal places required 17.38
c
data ath1cs( 1) / -.0712583781 5669365e0 /
data ath1cs( 2) / -.0059047197 9831451e0 /
data ath1cs( 3) / -.0001211454 4069499e0 /
data ath1cs( 4) / -.0000098860 8542270e0 /
data ath1cs( 5) / -.0000013808 4097352e0 /
data ath1cs( 6) / -.0000002614 2640172e0 /
data ath1cs( 7) / -.0000000605 0432589e0 /
data ath1cs( 8) / -.0000000161 8436223e0 /
data ath1cs( 9) / -.0000000048 3464911e0 /
data ath1cs(10) / -.0000000015 7655272e0 /
data ath1cs(11) / -.0000000005 5231518e0 /
data ath1cs(12) / -.0000000002 0545441e0 /
data ath1cs(13) / -.0000000000 8043412e0 /
data ath1cs(14) / -.0000000000 3291252e0 /
data ath1cs(15) / -.0000000000 1399875e0 /
data ath1cs(16) / -.0000000000 0616151e0 /
data ath1cs(17) / -.0000000000 0279614e0 /
data ath1cs(18) / -.0000000000 0130428e0 /
data ath1cs(19) / -.0000000000 0062373e0 /
data ath1cs(20) / -.0000000000 0030512e0 /
data ath1cs(21) / -.0000000000 0015239e0 /
data ath1cs(22) / -.0000000000 0007758e0 /
data ath1cs(23) / -.0000000000 0004020e0 /
data ath1cs(24) / -.0000000000 0002117e0 /
data ath1cs(25) / -.0000000000 0001132e0 /
data ath1cs(26) / -.0000000000 0000614e0 /
data ath1cs(27) / -.0000000000 0000337e0 /
data ath1cs(28) / -.0000000000 0000188e0 /
data ath1cs(29) / -.0000000000 0000105e0 /
data ath1cs(30) / -.0000000000 0000060e0 /
data ath1cs(31) / -.0000000000 0000034e0 /
data ath1cs(32) / -.0000000000 0000020e0 /
data ath1cs(33) / -.0000000000 0000011e0 /
data ath1cs(34) / -.0000000000 0000007e0 /
data ath1cs(35) / -.0000000000 0000004e0 /
data ath1cs(36) / -.0000000000 0000002e0 /
c
c series for am22 on the interval -1.00000d+00 to -1.25000d-01
c with weighted error 2.99e-17
c log weighted error 16.52
c significant figures required 14.57
c decimal places required 17.28
c
data am22cs( 1) / -.0156284448 0625341e0 /
data am22cs( 2) / .0077833644 5239681e0 /
data am22cs( 3) / .0008670577 7047718e0 /
data am22cs( 4) / .0001569662 7315611e0 /
data am22cs( 5) / .0000356396 2571432e0 /
data am22cs( 6) / .0000092459 8335425e0 /
data am22cs( 7) / .0000026211 0161850e0 /
data am22cs( 8) / .0000007918 8221651e0 /
data am22cs( 9) / .0000002510 4152792e0 /
data am22cs(10) / .0000000826 5223206e0 /
data am22cs(11) / .0000000280 5711662e0 /
data am22cs(12) / .0000000097 6821090e0 /
data am22cs(13) / .0000000034 7407923e0 /
data am22cs(14) / .0000000012 5828132e0 /
data am22cs(15) / .0000000004 6298826e0 /
data am22cs(16) / .0000000001 7272825e0 /
data am22cs(17) / .0000000000 6523192e0 /
data am22cs(18) / .0000000000 2490471e0 /
data am22cs(19) / .0000000000 0960156e0 /
data am22cs(20) / .0000000000 0373448e0 /
data am22cs(21) / .0000000000 0146417e0 /
data am22cs(22) / .0000000000 0057826e0 /
data am22cs(23) / .0000000000 0022991e0 /
data am22cs(24) / .0000000000 0009197e0 /
data am22cs(25) / .0000000000 0003700e0 /
data am22cs(26) / .0000000000 0001496e0 /
data am22cs(27) / .0000000000 0000608e0 /
data am22cs(28) / .0000000000 0000248e0 /
data am22cs(29) / .0000000000 0000101e0 /
data am22cs(30) / .0000000000 0000041e0 /
data am22cs(31) / .0000000000 0000017e0 /
data am22cs(32) / .0000000000 0000007e0 /
data am22cs(33) / .0000000000 0000002e0 /
c
c series for ath2 on the interval -1.00000d+00 to -1.25000d-01
c with weighted error 2.57e-17
c log weighted error 16.59
c significant figures required 15.07
c decimal places required 17.34
c
data ath2cs( 1) / .0044052734 5871877e0 /
data ath2cs( 2) / -.0304291945 2318455e0 /
data ath2cs( 3) / -.0013856532 8377179e0 /
data ath2cs( 4) / -.0001804443 9089549e0 /
data ath2cs( 5) / -.0000338084 7108327e0 /
data ath2cs( 6) / -.0000076781 8353522e0 /
data ath2cs( 7) / -.0000019678 3944371e0 /
data ath2cs( 8) / -.0000005483 7271158e0 /
data ath2cs( 9) / -.0000001625 4615505e0 /
data ath2cs(10) / -.0000000505 3049981e0 /
data ath2cs(11) / -.0000000163 1580701e0 /
data ath2cs(12) / -.0000000054 3420411e0 /
data ath2cs(13) / -.0000000018 5739855e0 /
data ath2cs(14) / -.0000000006 4895120e0 /
data ath2cs(15) / -.0000000002 3105948e0 /
data ath2cs(16) / -.0000000000 8363282e0 /
data ath2cs(17) / -.0000000000 3071196e0 /
data ath2cs(18) / -.0000000000 1142367e0 /
data ath2cs(19) / -.0000000000 0429811e0 /
data ath2cs(20) / -.0000000000 0163389e0 /
data ath2cs(21) / -.0000000000 0062693e0 /
data ath2cs(22) / -.0000000000 0024260e0 /
data ath2cs(23) / -.0000000000 0009461e0 /
data ath2cs(24) / -.0000000000 0003716e0 /
data ath2cs(25) / -.0000000000 0001469e0 /
data ath2cs(26) / -.0000000000 0000584e0 /
data ath2cs(27) / -.0000000000 0000233e0 /
data ath2cs(28) / -.0000000000 0000093e0 /
data ath2cs(29) / -.0000000000 0000037e0 /
data ath2cs(30) / -.0000000000 0000015e0 /
data ath2cs(31) / -.0000000000 0000006e0 /
data ath2cs(32) / -.0000000000 0000002e0 /
c
data pi4 / 0.7853981633 9744831 e0 /
data nam21, nath1, nam22, nath2, xsml / 4*0, 0.0 /
c
if (nam21.ne.0) go to 10
eta = 0.1*r1mach(3)
nam21 = inits (am21cs, 40, eta)
nath1 = inits (ath1cs, 36, eta)
nam22 = inits (am22cs, 33, eta)
nath2 = inits (ath2cs, 32, eta)
c
xsml = -(16.0/r1mach(3))**0.3333
c
10 if (x.ge.(-2.0)) go to 20
z = 1.0
if (x.gt.xsml) z = 16.0/x**3 + 1.0
ampl = 0.3125 + csevl(z, am21cs, nam21)
theta = -0.625 + csevl (z, ath1cs, nath1)
go to 30
c
20 if (x.gt.(-1.0)) call seteru (25hr9aimp x must be le -1.0, 25, 1,
1 2)
c
z = (16.0/x**3 + 9.0)/7.0
ampl = 0.3125 + csevl (z, am22cs, nam22)
theta = -0.625 + csevl (z, ath2cs, nath2)
c
30 sqrtx = sqrt(-x)
ampl = sqrt (ampl/sqrtx)
theta = pi4 - x*sqrtx * theta
c
return
end
|