bie.f 9.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206
  1. *DECK BIE
  2. FUNCTION BIE (X)
  3. C***BEGIN PROLOGUE BIE
  4. C***PURPOSE Calculate the Bairy function for a negative argument and an
  5. C exponentially scaled Bairy function for a non-negative
  6. C argument.
  7. C***LIBRARY SLATEC (FNLIB)
  8. C***CATEGORY C10D
  9. C***TYPE SINGLE PRECISION (BIE-S, DBIE-D)
  10. C***KEYWORDS BAIRY FUNCTION, EXPONENTIALLY SCALED, FNLIB,
  11. C SPECIAL FUNCTIONS
  12. C***AUTHOR Fullerton, W., (LANL)
  13. C***DESCRIPTION
  14. C
  15. C Evaluate BI(X) for X .LE. 0 and BI(X)*EXP(ZETA) where
  16. C ZETA = 2/3 * X**(3/2) for X .GE. 0.0
  17. C
  18. C Series for BIF on the interval -1.00000D+00 to 1.00000D+00
  19. C with weighted error 1.88E-19
  20. C log weighted error 18.72
  21. C significant figures required 17.74
  22. C decimal places required 19.20
  23. C
  24. C Series for BIG on the interval -1.00000D+00 to 1.00000D+00
  25. C with weighted error 2.61E-17
  26. C log weighted error 16.58
  27. C significant figures required 15.17
  28. C decimal places required 17.03
  29. C
  30. C Series for BIF2 on the interval 1.00000D+00 to 8.00000D+00
  31. C with weighted error 1.11E-17
  32. C log weighted error 16.95
  33. C approx significant figures required 16.5
  34. C decimal places required 17.45
  35. C
  36. C Series for BIG2 on the interval 1.00000D+00 to 8.00000D+00
  37. C with weighted error 1.19E-18
  38. C log weighted error 17.92
  39. C approx significant figures required 17.2
  40. C decimal places required 18.42
  41. C
  42. C Series for BIP on the interval 1.25000D-01 to 3.53553D-01
  43. C with weighted error 1.91E-17
  44. C log weighted error 16.72
  45. C significant figures required 15.35
  46. C decimal places required 17.41
  47. C
  48. C Series for BIP2 on the interval 0. to 1.25000D-01
  49. C with weighted error 1.05E-18
  50. C log weighted error 17.98
  51. C significant figures required 16.74
  52. C decimal places required 18.71
  53. C
  54. C***REFERENCES (NONE)
  55. C***ROUTINES CALLED CSEVL, INITS, R1MACH, R9AIMP
  56. C***REVISION HISTORY (YYMMDD)
  57. C 770701 DATE WRITTEN
  58. C 890206 REVISION DATE from Version 3.2
  59. C 891214 Prologue converted to Version 4.0 format. (BAB)
  60. C***END PROLOGUE BIE
  61. LOGICAL FIRST
  62. DIMENSION BIFCS(9), BIGCS(8), BIF2CS(10), BIG2CS(10), BIPCS(24),
  63. 1 BIP2CS(29)
  64. SAVE BIFCS, BIGCS, BIF2CS, BIG2CS, BIPCS, BIP2CS, ATR, BTR,
  65. 1 NBIF, NBIG, NBIF2, NBIG2, NBIP, NBIP2, X3SML, X32SML, XBIG, FIRST
  66. DATA BIFCS( 1) / -.0167302164 7198664948E0 /
  67. DATA BIFCS( 2) / .1025233583 424944561E0 /
  68. DATA BIFCS( 3) / .0017083092 5073815165E0 /
  69. DATA BIFCS( 4) / .0000118625 4546774468E0 /
  70. DATA BIFCS( 5) / .0000000449 3290701779E0 /
  71. DATA BIFCS( 6) / .0000000001 0698207143E0 /
  72. DATA BIFCS( 7) / .0000000000 0017480643E0 /
  73. DATA BIFCS( 8) / .0000000000 0000020810E0 /
  74. DATA BIFCS( 9) / .0000000000 0000000018E0 /
  75. DATA BIGCS( 1) / .0224662232 4857452E0 /
  76. DATA BIGCS( 2) / .0373647754 5301955E0 /
  77. DATA BIGCS( 3) / .0004447621 8957212E0 /
  78. DATA BIGCS( 4) / .0000024708 0756363E0 /
  79. DATA BIGCS( 5) / .0000000079 1913533E0 /
  80. DATA BIGCS( 6) / .0000000000 1649807E0 /
  81. DATA BIGCS( 7) / .0000000000 0002411E0 /
  82. DATA BIGCS( 8) / .0000000000 0000002E0 /
  83. DATA BIF2CS( 1) / 0.0998457269 3816041E0 /
  84. DATA BIF2CS( 2) / .4786249778 63005538E0 /
  85. DATA BIF2CS( 3) / .0251552119 604330118E0 /
  86. DATA BIF2CS( 4) / .0005820693 885232645E0 /
  87. DATA BIF2CS( 5) / .0000074997 659644377E0 /
  88. DATA BIF2CS( 6) / .0000000613 460287034E0 /
  89. DATA BIF2CS( 7) / .0000000003 462753885E0 /
  90. DATA BIF2CS( 8) / .0000000000 014288910E0 /
  91. DATA BIF2CS( 9) / .0000000000 000044962E0 /
  92. DATA BIF2CS(10) / .0000000000 000000111E0 /
  93. DATA BIG2CS( 1) / .0333056621 45514340E0 /
  94. DATA BIG2CS( 2) / .1613092151 23197068E0 /
  95. DATA BIG2CS( 3) / .0063190073 096134286E0 /
  96. DATA BIG2CS( 4) / .0001187904 568162517E0 /
  97. DATA BIG2CS( 5) / .0000013045 345886200E0 /
  98. DATA BIG2CS( 6) / .0000000093 741259955E0 /
  99. DATA BIG2CS( 7) / .0000000000 474580188E0 /
  100. DATA BIG2CS( 8) / .0000000000 001783107E0 /
  101. DATA BIG2CS( 9) / .0000000000 000005167E0 /
  102. DATA BIG2CS(10) / .0000000000 000000011E0 /
  103. DATA BIPCS( 1) / -.0832204747 7943447E0 /
  104. DATA BIPCS( 2) / .0114611892 7371174E0 /
  105. DATA BIPCS( 3) / .0004289644 0718911E0 /
  106. DATA BIPCS( 4) / -.0001490663 9379950E0 /
  107. DATA BIPCS( 5) / -.0000130765 9726787E0 /
  108. DATA BIPCS( 6) / .0000063275 9839610E0 /
  109. DATA BIPCS( 7) / -.0000004222 6696982E0 /
  110. DATA BIPCS( 8) / -.0000001914 7186298E0 /
  111. DATA BIPCS( 9) / .0000000645 3106284E0 /
  112. DATA BIPCS(10) / -.0000000078 4485467E0 /
  113. DATA BIPCS(11) / -.0000000009 6077216E0 /
  114. DATA BIPCS(12) / .0000000007 0004713E0 /
  115. DATA BIPCS(13) / -.0000000001 7731789E0 /
  116. DATA BIPCS(14) / .0000000000 2272089E0 /
  117. DATA BIPCS(15) / .0000000000 0165404E0 /
  118. DATA BIPCS(16) / -.0000000000 0185171E0 /
  119. DATA BIPCS(17) / .0000000000 0059576E0 /
  120. DATA BIPCS(18) / -.0000000000 0012194E0 /
  121. DATA BIPCS(19) / .0000000000 0001334E0 /
  122. DATA BIPCS(20) / .0000000000 0000172E0 /
  123. DATA BIPCS(21) / -.0000000000 0000145E0 /
  124. DATA BIPCS(22) / .0000000000 0000049E0 /
  125. DATA BIPCS(23) / -.0000000000 0000011E0 /
  126. DATA BIPCS(24) / .0000000000 0000001E0 /
  127. DATA BIP2CS( 1) / -.1135967375 85988679E0 /
  128. DATA BIP2CS( 2) / .0041381473 947881595E0 /
  129. DATA BIP2CS( 3) / .0001353470 622119332E0 /
  130. DATA BIP2CS( 4) / .0000104273 166530153E0 /
  131. DATA BIP2CS( 5) / .0000013474 954767849E0 /
  132. DATA BIP2CS( 6) / .0000001696 537405438E0 /
  133. DATA BIP2CS( 7) / -.0000000100 965008656E0 /
  134. DATA BIP2CS( 8) / -.0000000167 291194937E0 /
  135. DATA BIP2CS( 9) / -.0000000045 815364485E0 /
  136. DATA BIP2CS(10) / .0000000003 736681366E0 /
  137. DATA BIP2CS(11) / .0000000005 766930320E0 /
  138. DATA BIP2CS(12) / .0000000000 621812650E0 /
  139. DATA BIP2CS(13) / -.0000000000 632941202E0 /
  140. DATA BIP2CS(14) / -.0000000000 149150479E0 /
  141. DATA BIP2CS(15) / .0000000000 078896213E0 /
  142. DATA BIP2CS(16) / .0000000000 024960513E0 /
  143. DATA BIP2CS(17) / -.0000000000 012130075E0 /
  144. DATA BIP2CS(18) / -.0000000000 003740493E0 /
  145. DATA BIP2CS(19) / .0000000000 002237727E0 /
  146. DATA BIP2CS(20) / .0000000000 000474902E0 /
  147. DATA BIP2CS(21) / -.0000000000 000452616E0 /
  148. DATA BIP2CS(22) / -.0000000000 000030172E0 /
  149. DATA BIP2CS(23) / .0000000000 000091058E0 /
  150. DATA BIP2CS(24) / -.0000000000 000009814E0 /
  151. DATA BIP2CS(25) / -.0000000000 000016429E0 /
  152. DATA BIP2CS(26) / .0000000000 000005533E0 /
  153. DATA BIP2CS(27) / .0000000000 000002175E0 /
  154. DATA BIP2CS(28) / -.0000000000 000001737E0 /
  155. DATA BIP2CS(29) / -.0000000000 000000010E0 /
  156. DATA ATR / 8.750690570 8484345 E0 /
  157. DATA BTR / -2.093836321 356054 E0 /
  158. DATA FIRST /.TRUE./
  159. C***FIRST EXECUTABLE STATEMENT BIE
  160. IF (FIRST) THEN
  161. ETA = 0.1*R1MACH(3)
  162. NBIF = INITS (BIFCS, 9, ETA)
  163. NBIG = INITS (BIGCS, 8, ETA)
  164. NBIF2 = INITS (BIF2CS, 10, ETA)
  165. NBIG2 = INITS (BIG2CS, 10, ETA)
  166. NBIP = INITS (BIPCS , 24, ETA)
  167. NBIP2 = INITS (BIP2CS, 29, ETA)
  168. C
  169. X3SML = ETA**0.3333
  170. X32SML = 1.3104*X3SML**2
  171. XBIG = R1MACH(2)**0.6666
  172. ENDIF
  173. FIRST = .FALSE.
  174. C
  175. IF (X.GE.(-1.0)) GO TO 20
  176. CALL R9AIMP (X, XM, THETA)
  177. BIE = XM * SIN(THETA)
  178. RETURN
  179. C
  180. 20 IF (X.GT.1.0) GO TO 30
  181. Z = 0.0
  182. IF (ABS(X).GT.X3SML) Z = X**3
  183. BIE = 0.625 + CSEVL (Z, BIFCS, NBIF) + X*(0.4375 +
  184. 1 CSEVL (Z, BIGCS, NBIG))
  185. IF (X.GT.X32SML) BIE = BIE * EXP(-2.0*X*SQRT(X)/3.0)
  186. RETURN
  187. C
  188. 30 IF (X.GT.2.0) GO TO 40
  189. Z = (2.0*X**3 - 9.0) / 7.0
  190. BIE = EXP(-2.0*X*SQRT(X)/3.0) * (1.125 + CSEVL (Z, BIF2CS, NBIF2)
  191. 1 + X*(0.625 + CSEVL (Z, BIG2CS, NBIG2)) )
  192. RETURN
  193. C
  194. 40 IF (X.GT.4.0) GO TO 50
  195. SQRTX = SQRT(X)
  196. Z = ATR/(X*SQRTX) + BTR
  197. BIE = (0.625 + CSEVL (Z, BIPCS, NBIP)) / SQRT(SQRTX)
  198. RETURN
  199. C
  200. 50 SQRTX = SQRT(X)
  201. Z = -1.0
  202. IF (X.LT.XBIG) Z = 16.0/(X*SQRTX) - 1.0
  203. BIE = (0.625 + CSEVL (Z, BIP2CS, NBIP2))/SQRT(SQRTX)
  204. RETURN
  205. C
  206. END