Blame view
fvn_fnlib/daie.f
11.4 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 |
double precision function daie (x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate ai(x) for x .le. 0.0 and ai(x)*exp(zeta) where
c zeta = 2/3 * x**(3/2) for x .ge. 0.0
c
double precision x, aifcs(13), aigcs(13), aip1cs(57), aip2cs(37),
1 sqrtx, theta, xbig, xm, x3sml, x32sml, z, d1mach, dcsevl,
2 dcos, dexp, dsqrt
external d1mach, 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
c series for aip1 on the interval 1.25000e-01 to 1.00000e+00
c with weighted error 3.69e-32
c log weighted error 31.43
c significant figures required 29.55
c decimal places required 32.31
c
data aip1cs( 1) / -.2146951858 9105384554 6086346777 8 d-1 /
data aip1cs( 2) / -.7535382535 0433011662 1972086556 5 d-2 /
data aip1cs( 3) / +.5971527949 0263808520 3538888199 4 d-3 /
data aip1cs( 4) / -.7283251254 2076106485 0236829154 8 d-4 /
data aip1cs( 5) / +.1110297130 7392996665 1738182114 0 d-4 /
data aip1cs( 6) / -.1950386152 2844057103 4693031403 3 d-5 /
data aip1cs( 7) / +.3786973885 1595151938 8531967005 7 d-6 /
data aip1cs( 8) / -.7929675297 3509782790 3907287915 4 d-7 /
data aip1cs( 9) / +.1762247638 6742560755 6842012220 2 d-7 /
data aip1cs( 10) / -.4110767539 6671950450 2989659389 3 d-8 /
data aip1cs( 11) / +.9984770057 8578922471 8341410754 4 d-9 /
data aip1cs( 12) / -.2510093251 3871222113 4986773003 4 d-9 /
data aip1cs( 13) / +.6500501929 8606954092 7203860172 5 d-10 /
data aip1cs( 14) / -.1727818405 3936165154 7887710736 6 d-10 /
data aip1cs( 15) / +.4699378842 8245125783 6229287230 7 d-11 /
data aip1cs( 16) / -.1304675656 2977439144 9124124627 2 d-11 /
data aip1cs( 17) / +.3689698478 4626788104 7394838228 2 d-12 /
data aip1cs( 18) / -.1061087206 6468061736 5035967903 5 d-12 /
data aip1cs( 19) / +.3098414384 8781874386 6021007011 0 d-13 /
data aip1cs( 20) / -.9174908079 8241393078 3342354785 1 d-14 /
data aip1cs( 21) / +.2752049140 3472108956 9357906227 1 d-14 /
data aip1cs( 22) / -.8353750115 9220465580 9139330188 0 d-15 /
data aip1cs( 23) / +.2563931129 3579349475 6863616861 2 d-15 /
data aip1cs( 24) / -.7950633762 5988549832 7374728982 2 d-16 /
data aip1cs( 25) / +.2489283634 6030699774 3728117564 4 d-16 /
data aip1cs( 26) / -.7864326933 9287355696 6462622129 6 d-17 /
data aip1cs( 27) / +.2505687311 4399756723 2447064501 9 d-17 /
data aip1cs( 28) / -.8047420364 1639095245 3795868224 1 d-18 /
data aip1cs( 29) / +.2604097118 9520539644 4340110439 2 d-18 /
data aip1cs( 30) / -.8486954164 0564122594 8248883418 4 d-19 /
data aip1cs( 31) / +.2784706882 1423378433 5942918602 7 d-19 /
data aip1cs( 32) / -.9195858953 4986129136 8722415135 4 d-20 /
data aip1cs( 33) / +.3055304318 3742387422 4766822558 3 d-20 /
data aip1cs( 34) / -.1021035455 4794778759 0217704843 9 d-20 /
data aip1cs( 35) / +.3431118190 7437578440 0055568083 6 d-21 /
data aip1cs( 36) / -.1159129341 7977495133 7692246310 9 d-21 /
data aip1cs( 37) / +.3935772844 2002556108 3626822915 4 d-22 /
data aip1cs( 38) / -.1342880980 2967176119 5671898903 8 d-22 /
data aip1cs( 39) / +.4603287883 5200027416 5919030531 4 d-23 /
data aip1cs( 40) / -.1585043927 0040642278 1077249938 7 d-23 /
data aip1cs( 41) / +.5481275667 7296759089 2552375500 8 d-24 /
data aip1cs( 42) / -.1903349371 8550472590 6401794894 5 d-24 /
data aip1cs( 43) / +.6635682302 3740087167 7761211596 8 d-25 /
data aip1cs( 44) / -.2322311650 0263143079 7520098645 3 d-25 /
data aip1cs( 45) / +.8157640113 4291793131 4274369535 9 d-26 /
data aip1cs( 46) / -.2875824240 6329004900 5748992955 7 d-26 /
data aip1cs( 47) / +.1017329450 9429014350 7971431901 8 d-26 /
data aip1cs( 48) / -.3610879108 7422164465 7570349055 9 d-27 /
data aip1cs( 49) / +.1285788540 3639934212 5664034269 8 d-27 /
data aip1cs( 50) / -.4592901037 3785474251 6069302271 9 d-28 /
data aip1cs( 51) / +.1645597033 8207137258 1210248533 3 d-28 /
data aip1cs( 52) / -.5913421299 8435018420 8792027136 0 d-29 /
data aip1cs( 53) / +.2131057006 6049933034 7936950954 6 d-29 /
data aip1cs( 54) / -.7701158157 7875982169 8276174506 6 d-30 /
data aip1cs( 55) / +.2790533307 9689304175 8178377728 0 d-30 /
data aip1cs( 56) / -.1013807715 1112840064 5224136703 9 d-30 /
data aip1cs( 57) / +.3692580158 7196240936 5828621653 3 d-31 /
c
c series for aip2 on the interval 0. to 1.25000e-01
c with weighted error 3.48e-32
c log weighted error 31.46
c significant figures required 28.74
c decimal places required 32.24
c
data aip2cs( 1) / -.1743144969 2937551339 0355844011 d-2 /
data aip2cs( 2) / -.1678938543 2554167163 2190613480 d-2 /
data aip2cs( 3) / +.3596534033 5216603588 5983858114 d-4 /
data aip2cs( 4) / -.1380818602 7392283545 7399383100 d-5 /
data aip2cs( 5) / +.7411228077 3150529884 8699095233 d-7 /
data aip2cs( 6) / -.5002382039 0013301313 0422866325 d-8 /
data aip2cs( 7) / +.4006939174 1718424067 5446866355 d-9 /
data aip2cs( 8) / -.3673312427 9590504419 9318496207 d-10 /
data aip2cs( 9) / +.3760344395 9237385243 9592002918 d-11 /
data aip2cs( 10) / -.4223213327 1874753802 6564938968 d-12 /
data aip2cs( 11) / +.5135094540 3365707091 9618754120 d-13 /
data aip2cs( 12) / -.6690958503 9047759565 1681356676 d-14 /
data aip2cs( 13) / +.9266675456 4129064823 9550724382 d-15 /
data aip2cs( 14) / -.1355143824 1607057633 3397356591 d-15 /
data aip2cs( 15) / +.2081154963 1283099529 9006549335 d-16 /
data aip2cs( 16) / -.3341164991 5917685687 1277570256 d-17 /
data aip2cs( 17) / +.5585785845 8592431686 8032946585 d-18 /
data aip2cs( 18) / -.9692190401 5236524751 8658209109 d-19 /
data aip2cs( 19) / +.1740457001 2889320646 5696557738 d-19 /
data aip2cs( 20) / -.3226409797 3113040024 7846333098 d-20 /
data aip2cs( 21) / +.6160744711 0662525853 3259618986 d-21 /
data aip2cs( 22) / -.1209363479 8249005907 6420676266 d-21 /
data aip2cs( 23) / +.2436327633 1013810826 1570095786 d-22 /
data aip2cs( 24) / -.5029142214 9745746894 3403144533 d-23 /
data aip2cs( 25) / +.1062241755 4363568949 5470626133 d-23 /
data aip2cs( 26) / -.2292842848 9598924150 9856324266 d-24 /
data aip2cs( 27) / +.5051817339 2950374498 6884778666 d-25 /
data aip2cs( 28) / -.1134981237 1441240497 9793920000 d-25 /
data aip2cs( 29) / +.2597655659 8560698069 8374144000 d-26 /
data aip2cs( 30) / -.6051246215 4293950617 2231679999 d-27 /
data aip2cs( 31) / +.1433597779 6677280072 0295253333 d-27 /
data aip2cs( 32) / -.3451477570 6089998628 0721066666 d-28 /
data aip2cs( 33) / +.8438751902 1364674042 7025066666 d-29 /
data aip2cs( 34) / -.2093961422 9818816943 4453333333 d-29 /
data aip2cs( 35) / +.5270088734 7894550318 2848000000 d-30 /
data aip2cs( 36) / -.1344574330 1455338578 9030399999 d-30 /
data aip2cs( 37) / +.3475709645 2660114734 0117333333 d-31 /
c
data naif, naig, naip1, naip2 / 4*0 /
data x3sml, x32sml, xbig / 3*0.d0 /
c
if (naif.ne.0) go to 10
eta = 0.1*sngl(d1mach(3))
naif = initds (aifcs, 13, eta)
naig = initds (aigcs, 13, eta)
naip1 = initds (aip1cs, 57, eta)
naip2 = initds (aip2cs, 37, eta)
c
x3sml = eta**0.3333e0
x32sml = 1.3104d0*x3sml**2
xbig = d1mach(2)**0.6666d0
c
10 if (x.ge.(-1.0d0)) go to 20
call d9aimp (x, xm, theta)
daie = 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
daie = 0.375d0 + (dcsevl (z, aifcs, naif) - x*(0.25d0 +
1 dcsevl (z, aigcs, naig)) )
if (x.gt.x32sml) daie = daie * dexp (2.0d0*x*dsqrt(x)/3.0d0)
return
c
30 if (x.gt.4.0d0) go to 40
sqrtx = dsqrt(x)
z = (16.d0/(x*sqrtx) - 9.d0)/7.d0
daie = (0.28125d0 + dcsevl (z, aip1cs, naip1))/dsqrt(sqrtx)
return
c
40 sqrtx = dsqrt(x)
z = -1.0d0
if (x.lt.xbig) z = 16.0d0/(x*sqrtx) - 1.0d0
daie = (0.28125d0 + dcsevl (z, aip2cs, naip2))/dsqrt(sqrtx)
return
c
end
|