fresnl.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372
  1. /* fresnl.c
  2. *
  3. * Fresnel integral
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, S, C;
  10. * void fresnl();
  11. *
  12. * fresnl( x, _&S, _&C );
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Evaluates the Fresnel integrals
  18. *
  19. * x
  20. * -
  21. * | |
  22. * C(x) = | cos(pi/2 t**2) dt,
  23. * | |
  24. * -
  25. * 0
  26. *
  27. * x
  28. * -
  29. * | |
  30. * S(x) = | sin(pi/2 t**2) dt.
  31. * | |
  32. * -
  33. * 0
  34. *
  35. *
  36. * The integrals are evaluated by a power series for x < 1.
  37. * For x >= 1 auxiliary functions f(x) and g(x) are employed
  38. * such that
  39. *
  40. * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
  41. * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
  42. *
  43. *
  44. *
  45. * ACCURACY:
  46. *
  47. * Relative error.
  48. *
  49. * Arithmetic function domain # trials peak rms
  50. * IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16
  51. * IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16
  52. * DEC S(x) 0, 10 6000 2.2e-16 3.9e-17
  53. * DEC C(x) 0, 10 5000 2.3e-16 3.9e-17
  54. */
  55. /*
  56. Cephes Math Library Release 2.8: June, 2000
  57. Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
  58. */
  59. #include "mconf.h"
  60. /* S(x) for small x */
  61. #ifdef UNK
  62. static double sn[6] = {
  63. -2.99181919401019853726E3, 7.08840045257738576863E5,
  64. -6.29741486205862506537E7, 2.54890880573376359104E9,
  65. -4.42979518059697779103E10, 3.18016297876567817986E11,
  66. };
  67. static double sd[6] = {
  68. /* 1.00000000000000000000E0,*/
  69. 2.81376268889994315696E2, 4.55847810806532581675E4,
  70. 5.17343888770096400730E6, 4.19320245898111231129E8,
  71. 2.24411795645340920940E10, 6.07366389490084639049E11,
  72. };
  73. #endif
  74. #ifdef DEC
  75. static unsigned short sn[24] = {
  76. 0143072, 0176433, 0065455, 0127034, 0045055, 0007200, 0134540, 0026661,
  77. 0146560, 0035061, 0023667, 0127545, 0050027, 0166503, 0002673, 0153756,
  78. 0151045, 0002721, 0121737, 0102066, 0051624, 0013177, 0033451, 0021271,
  79. };
  80. static unsigned short sd[24] = {
  81. /*0040200,0000000,0000000,0000000,*/
  82. 0042214, 0130051, 0112070, 0101617, 0044062, 0010307, 0172346, 0152510,
  83. 0045635, 0160575, 0143200, 0136642, 0047307, 0171215, 0127457, 0052361,
  84. 0050647, 0031447, 0032621, 0013510, 0052015, 0064733, 0117362, 0012653,
  85. };
  86. #endif
  87. #ifdef IBMPC
  88. static unsigned short sn[24] = {
  89. 0xb5c3, 0x6d65, 0x5fa3, 0xc0a7, 0x05b6, 0x172c, 0xa1d0, 0x4125,
  90. 0xf5ed, 0x24f6, 0x0746, 0xc18e, 0x7afe, 0x60b7, 0xfda8, 0x41e2,
  91. 0xf087, 0x347b, 0xa0ba, 0xc224, 0x2457, 0xe6e5, 0x82cf, 0x4252,
  92. };
  93. static unsigned short sd[24] = {
  94. /*0x0000,0x0000,0x0000,0x3ff0,*/
  95. 0x1072, 0x3287, 0x9605, 0x4071, 0xdaa9, 0xfe9c, 0x4218, 0x40e6,
  96. 0x17b4, 0xb8d0, 0xbc2f, 0x4153, 0xea9e, 0xb5e5, 0xfe51, 0x41b8,
  97. 0x22e9, 0xe6b2, 0xe664, 0x4214, 0x42b5, 0x73de, 0xad3b, 0x4261,
  98. };
  99. #endif
  100. #ifdef MIEEE
  101. static unsigned short sn[24] = {
  102. 0xc0a7, 0x5fa3, 0x6d65, 0xb5c3, 0x4125, 0xa1d0, 0x172c, 0x05b6,
  103. 0xc18e, 0x0746, 0x24f6, 0xf5ed, 0x41e2, 0xfda8, 0x60b7, 0x7afe,
  104. 0xc224, 0xa0ba, 0x347b, 0xf087, 0x4252, 0x82cf, 0xe6e5, 0x2457,
  105. };
  106. static unsigned short sd[24] = {
  107. /*0x3ff0,0x0000,0x0000,0x0000,*/
  108. 0x4071, 0x9605, 0x3287, 0x1072, 0x40e6, 0x4218, 0xfe9c, 0xdaa9,
  109. 0x4153, 0xbc2f, 0xb8d0, 0x17b4, 0x41b8, 0xfe51, 0xb5e5, 0xea9e,
  110. 0x4214, 0xe664, 0xe6b2, 0x22e9, 0x4261, 0xad3b, 0x73de, 0x42b5,
  111. };
  112. #endif
  113. /* C(x) for small x */
  114. #ifdef UNK
  115. static double cn[6] = {
  116. -4.98843114573573548651E-8, 9.50428062829859605134E-6,
  117. -6.45191435683965050962E-4, 1.88843319396703850064E-2,
  118. -2.05525900955013891793E-1, 9.99999999999999998822E-1,
  119. };
  120. static double cd[7] = {
  121. 3.99982968972495980367E-12, 9.15439215774657478799E-10,
  122. 1.25001862479598821474E-7, 1.22262789024179030997E-5,
  123. 8.68029542941784300606E-4, 4.12142090722199792936E-2,
  124. 1.00000000000000000118E0,
  125. };
  126. #endif
  127. #ifdef DEC
  128. static unsigned short cn[24] = {
  129. 0132126, 0040141, 0063733, 0013231, 0034037, 0072223, 0010200, 0075637,
  130. 0135451, 0021020, 0073264, 0036057, 0036632, 0131520, 0101316, 0060233,
  131. 0137522, 0072541, 0136124, 0132202, 0040200, 0000000, 0000000, 0000000,
  132. };
  133. static unsigned short cd[28] = {
  134. 0026614, 0135503, 0051776, 0032631, 0030573, 0121116, 0154033,
  135. 0126712, 0032406, 0034100, 0012442, 0106212, 0034115, 0017567,
  136. 0150520, 0164623, 0035543, 0106171, 0177336, 0146351, 0037050,
  137. 0150073, 0000607, 0171635, 0040200, 0000000, 0000000, 0000000,
  138. };
  139. #endif
  140. #ifdef IBMPC
  141. static unsigned short cn[24] = {
  142. 0x62d3, 0x2cfb, 0xc80c, 0xbe6a, 0x0f74, 0x6210, 0xee92, 0x3ee3,
  143. 0x8786, 0x0ed6, 0x2442, 0xbf45, 0xcc13, 0x1059, 0x566a, 0x3f93,
  144. 0x9690, 0x378a, 0x4eac, 0xbfca, 0x0000, 0x0000, 0x0000, 0x3ff0,
  145. };
  146. static unsigned short cd[28] = {
  147. 0xc6b3, 0x6a7f, 0x9768, 0x3d91, 0x75b9, 0xdb03, 0x7449,
  148. 0x3e0f, 0x5191, 0x02a4, 0xc708, 0x3e80, 0x1d32, 0xfa2a,
  149. 0xa3ee, 0x3ee9, 0xd99d, 0x3fdb, 0x718f, 0x3f4c, 0xfe74,
  150. 0x6030, 0x1a07, 0x3fa5, 0x0000, 0x0000, 0x0000, 0x3ff0,
  151. };
  152. #endif
  153. #ifdef MIEEE
  154. static unsigned short cn[24] = {
  155. 0xbe6a, 0xc80c, 0x2cfb, 0x62d3, 0x3ee3, 0xee92, 0x6210, 0x0f74,
  156. 0xbf45, 0x2442, 0x0ed6, 0x8786, 0x3f93, 0x566a, 0x1059, 0xcc13,
  157. 0xbfca, 0x4eac, 0x378a, 0x9690, 0x3ff0, 0x0000, 0x0000, 0x0000,
  158. };
  159. static unsigned short cd[28] = {
  160. 0x3d91, 0x9768, 0x6a7f, 0xc6b3, 0x3e0f, 0x7449, 0xdb03,
  161. 0x75b9, 0x3e80, 0xc708, 0x02a4, 0x5191, 0x3ee9, 0xa3ee,
  162. 0xfa2a, 0x1d32, 0x3f4c, 0x718f, 0x3fdb, 0xd99d, 0x3fa5,
  163. 0x1a07, 0x6030, 0xfe74, 0x3ff0, 0x0000, 0x0000, 0x0000,
  164. };
  165. #endif
  166. /* Auxiliary function f(x) */
  167. #ifdef UNK
  168. static double fn[10] = {
  169. 4.21543555043677546506E-1, 1.43407919780758885261E-1,
  170. 1.15220955073585758835E-2, 3.45017939782574027900E-4,
  171. 4.63613749287867322088E-6, 3.05568983790257605827E-8,
  172. 1.02304514164907233465E-10, 1.72010743268161828879E-13,
  173. 1.34283276233062758925E-16, 3.76329711269987889006E-20,
  174. };
  175. static double fd[10] = {
  176. /* 1.00000000000000000000E0,*/
  177. 7.51586398353378947175E-1, 1.16888925859191382142E-1,
  178. 6.44051526508858611005E-3, 1.55934409164153020873E-4,
  179. 1.84627567348930545870E-6, 1.12699224763999035261E-8,
  180. 3.60140029589371370404E-11, 5.88754533621578410010E-14,
  181. 4.52001434074129701496E-17, 1.25443237090011264384E-20,
  182. };
  183. #endif
  184. #ifdef DEC
  185. static unsigned short fn[40] = {
  186. 0037727, 0152216, 0106601, 0016214, 0037422, 0154606, 0112710, 0071355,
  187. 0036474, 0143453, 0154253, 0166545, 0035264, 0161606, 0022250, 0073743,
  188. 0033633, 0110036, 0024653, 0136246, 0032003, 0036652, 0041164, 0036413,
  189. 0027740, 0174122, 0046305, 0036726, 0025501, 0125270, 0121317, 0167667,
  190. 0023032, 0150555, 0076175, 0047443, 0020061, 0133570, 0070130, 0027657,
  191. };
  192. static unsigned short fd[40] = {
  193. /*0040200,0000000,0000000,0000000,*/
  194. 0040100, 0063767, 0054413, 0151452, 0037357, 0061566, 0007243, 0065754,
  195. 0036323, 0005365, 0033552, 0133625, 0035043, 0101123, 0000275, 0165402,
  196. 0033367, 0146614, 0110623, 0023647, 0031501, 0116644, 0125222, 0144263,
  197. 0027436, 0062051, 0117235, 0001411, 0025204, 0111543, 0056370, 0036201,
  198. 0022520, 0071351, 0015227, 0122144, 0017554, 0172240, 0112713, 0005006,
  199. };
  200. #endif
  201. #ifdef IBMPC
  202. static unsigned short fn[40] = {
  203. 0x2391, 0xd1b0, 0xfa91, 0x3fda, 0x0e5e, 0xd2b9, 0x5b30, 0x3fc2,
  204. 0x7dad, 0x7b15, 0x98e5, 0x3f87, 0x0efc, 0xc495, 0x9c70, 0x3f36,
  205. 0x7795, 0xc535, 0x7203, 0x3ed3, 0x87a1, 0x484e, 0x67b5, 0x3e60,
  206. 0xa7bb, 0x4998, 0x1f0a, 0x3ddc, 0xfdf7, 0x1459, 0x3557, 0x3d48,
  207. 0xa9e4, 0xaf8f, 0x5a2d, 0x3ca3, 0x05f6, 0x0e0b, 0x36ef, 0x3be6,
  208. };
  209. static unsigned short fd[40] = {
  210. /*0x0000,0x0000,0x0000,0x3ff0,*/
  211. 0x7a65, 0xeb21, 0x0cfe, 0x3fe8, 0x6d7d, 0xc1d4, 0xec6e, 0x3fbd,
  212. 0x56f3, 0xa6ed, 0x615e, 0x3f7a, 0xbd60, 0x6017, 0x704a, 0x3f24,
  213. 0x64f5, 0x9232, 0xf9b1, 0x3ebe, 0x5916, 0x9552, 0x33b4, 0x3e48,
  214. 0xa061, 0x33d3, 0xcc85, 0x3dc3, 0x0790, 0x6b9f, 0x926c, 0x3d30,
  215. 0xf48d, 0x2352, 0x0e5d, 0x3c8a, 0x6141, 0x12b9, 0x9e94, 0x3bcd,
  216. };
  217. #endif
  218. #ifdef MIEEE
  219. static unsigned short fn[40] = {
  220. 0x3fda, 0xfa91, 0xd1b0, 0x2391, 0x3fc2, 0x5b30, 0xd2b9, 0x0e5e,
  221. 0x3f87, 0x98e5, 0x7b15, 0x7dad, 0x3f36, 0x9c70, 0xc495, 0x0efc,
  222. 0x3ed3, 0x7203, 0xc535, 0x7795, 0x3e60, 0x67b5, 0x484e, 0x87a1,
  223. 0x3ddc, 0x1f0a, 0x4998, 0xa7bb, 0x3d48, 0x3557, 0x1459, 0xfdf7,
  224. 0x3ca3, 0x5a2d, 0xaf8f, 0xa9e4, 0x3be6, 0x36ef, 0x0e0b, 0x05f6,
  225. };
  226. static unsigned short fd[40] = {
  227. /*0x3ff0,0x0000,0x0000,0x0000,*/
  228. 0x3fe8, 0x0cfe, 0xeb21, 0x7a65, 0x3fbd, 0xec6e, 0xc1d4, 0x6d7d,
  229. 0x3f7a, 0x615e, 0xa6ed, 0x56f3, 0x3f24, 0x704a, 0x6017, 0xbd60,
  230. 0x3ebe, 0xf9b1, 0x9232, 0x64f5, 0x3e48, 0x33b4, 0x9552, 0x5916,
  231. 0x3dc3, 0xcc85, 0x33d3, 0xa061, 0x3d30, 0x926c, 0x6b9f, 0x0790,
  232. 0x3c8a, 0x0e5d, 0x2352, 0xf48d, 0x3bcd, 0x9e94, 0x12b9, 0x6141,
  233. };
  234. #endif
  235. /* Auxiliary function g(x) */
  236. #ifdef UNK
  237. static double gn[11] = {
  238. 5.04442073643383265887E-1, 1.97102833525523411709E-1,
  239. 1.87648584092575249293E-2, 6.84079380915393090172E-4,
  240. 1.15138826111884280931E-5, 9.82852443688422223854E-8,
  241. 4.45344415861750144738E-10, 1.08268041139020870318E-12,
  242. 1.37555460633261799868E-15, 8.36354435630677421531E-19,
  243. 1.86958710162783235106E-22,
  244. };
  245. static double gd[11] = {
  246. /* 1.00000000000000000000E0,*/
  247. 1.47495759925128324529E0, 3.37748989120019970451E-1,
  248. 2.53603741420338795122E-2, 8.14679107184306179049E-4,
  249. 1.27545075667729118702E-5, 1.04314589657571990585E-7,
  250. 4.60680728146520428211E-10, 1.10273215066240270757E-12,
  251. 1.38796531259578871258E-15, 8.39158816283118707363E-19,
  252. 1.86958710162783236342E-22,
  253. };
  254. #endif
  255. #ifdef DEC
  256. static unsigned short gn[44] = {
  257. 0040001, 0021435, 0120406, 0053123, 0037511, 0152523, 0037703, 0122011,
  258. 0036631, 0134302, 0122721, 0110235, 0035463, 0051712, 0043215, 0114732,
  259. 0034101, 0025677, 0147725, 0057630, 0032323, 0010342, 0067523, 0002206,
  260. 0030364, 0152247, 0110007, 0054107, 0026230, 0057654, 0035464, 0047124,
  261. 0023706, 0036401, 0167705, 0045440, 0021166, 0154447, 0105632, 0142461,
  262. 0016142, 0002353, 0011175, 0170530,
  263. };
  264. static unsigned short gd[44] = {
  265. /*0040200,0000000,0000000,0000000,*/
  266. 0040274, 0145551, 0016742, 0127005, 0037654, 0166557, 0076416, 0015165,
  267. 0036717, 0140217, 0030675, 0050111, 0035525, 0110060, 0076405, 0070502,
  268. 0034125, 0176061, 0060120, 0031730, 0032340, 0001615, 0054343, 0120501,
  269. 0030375, 0041414, 0070747, 0107060, 0026233, 0031034, 0160757, 0074526,
  270. 0023710, 0003341, 0137100, 0144664, 0021167, 0126414, 0023774, 0015435,
  271. 0016142, 0002353, 0011175, 0170530,
  272. };
  273. #endif
  274. #ifdef IBMPC
  275. static unsigned short gn[44] = {
  276. 0xcaca, 0xb420, 0x2463, 0x3fe0, 0x7481, 0x67f8, 0x3aaa, 0x3fc9, 0x3214,
  277. 0x54ba, 0x3718, 0x3f93, 0xb33b, 0x48d1, 0x6a79, 0x3f46, 0xabf3, 0xf9fa,
  278. 0x2577, 0x3ee8, 0x6091, 0x4dea, 0x621c, 0x3e7a, 0xeb09, 0xf200, 0x9a94,
  279. 0x3dfe, 0x89cb, 0x8766, 0x0bf5, 0x3d73, 0xa964, 0x3df8, 0xc7a0, 0x3cd8,
  280. 0x58a6, 0xf173, 0xdb24, 0x3c2e, 0xbe2b, 0x624f, 0x409d, 0x3b6c,
  281. };
  282. static unsigned short gd[44] = {
  283. /*0x0000,0x0000,0x0000,0x3ff0,*/
  284. 0x55c1, 0x23bc, 0x996d, 0x3ff7, 0xc34f, 0xefa1, 0x9dad, 0x3fd5, 0xaa09,
  285. 0xe637, 0xf811, 0x3f99, 0xae28, 0x0fa0, 0xb206, 0x3f4a, 0x067b, 0x2c0a,
  286. 0xbf86, 0x3eea, 0x7428, 0xab1c, 0x0071, 0x3e7c, 0xf1c6, 0x8e3c, 0xa861,
  287. 0x3dff, 0xef2b, 0x9c3d, 0x6643, 0x3d73, 0x1936, 0x37c8, 0x00dc, 0x3cd9,
  288. 0x8364, 0x84ff, 0xf5a1, 0x3c2e, 0xbe2b, 0x624f, 0x409d, 0x3b6c,
  289. };
  290. #endif
  291. #ifdef MIEEE
  292. static unsigned short gn[44] = {
  293. 0x3fe0, 0x2463, 0xb420, 0xcaca, 0x3fc9, 0x3aaa, 0x67f8, 0x7481, 0x3f93,
  294. 0x3718, 0x54ba, 0x3214, 0x3f46, 0x6a79, 0x48d1, 0xb33b, 0x3ee8, 0x2577,
  295. 0xf9fa, 0xabf3, 0x3e7a, 0x621c, 0x4dea, 0x6091, 0x3dfe, 0x9a94, 0xf200,
  296. 0xeb09, 0x3d73, 0x0bf5, 0x8766, 0x89cb, 0x3cd8, 0xc7a0, 0x3df8, 0xa964,
  297. 0x3c2e, 0xdb24, 0xf173, 0x58a6, 0x3b6c, 0x409d, 0x624f, 0xbe2b,
  298. };
  299. static unsigned short gd[44] = {
  300. /*0x3ff0,0x0000,0x0000,0x0000,*/
  301. 0x3ff7, 0x996d, 0x23bc, 0x55c1, 0x3fd5, 0x9dad, 0xefa1, 0xc34f, 0x3f99,
  302. 0xf811, 0xe637, 0xaa09, 0x3f4a, 0xb206, 0x0fa0, 0xae28, 0x3eea, 0xbf86,
  303. 0x2c0a, 0x067b, 0x3e7c, 0x0071, 0xab1c, 0x7428, 0x3dff, 0xa861, 0x8e3c,
  304. 0xf1c6, 0x3d73, 0x6643, 0x9c3d, 0xef2b, 0x3cd9, 0x00dc, 0x37c8, 0x1936,
  305. 0x3c2e, 0xf5a1, 0x84ff, 0x8364, 0x3b6c, 0x409d, 0x624f, 0xbe2b,
  306. };
  307. #endif
  308. #ifdef ANSIPROT
  309. extern double fabs(double);
  310. extern double cos(double);
  311. extern double sin(double);
  312. extern double polevl(double, void *, int);
  313. extern double p1evl(double, void *, int);
  314. #else
  315. double fabs(), cos(), sin(), polevl(), p1evl();
  316. #endif
  317. extern double PI, PIO2, MACHEP;
  318. int fresnl(xxa, ssa, cca) double xxa, *ssa, *cca;
  319. {
  320. double f, g, cc, ss, c, s, t, u;
  321. double x, x2;
  322. x = fabs(xxa);
  323. x2 = x * x;
  324. if (x2 < 2.5625) {
  325. t = x2 * x2;
  326. ss = x * x2 * polevl(t, sn, 5) / p1evl(t, sd, 6);
  327. cc = x * polevl(t, cn, 5) / polevl(t, cd, 6);
  328. goto done;
  329. }
  330. if (x > 36974.0) {
  331. cc = 0.5;
  332. ss = 0.5;
  333. goto done;
  334. }
  335. /* Asymptotic power series auxiliary functions
  336. * for large argument
  337. */
  338. x2 = x * x;
  339. t = PI * x2;
  340. u = 1.0 / (t * t);
  341. t = 1.0 / t;
  342. f = 1.0 - u * polevl(u, fn, 9) / p1evl(u, fd, 10);
  343. g = t * polevl(u, gn, 10) / p1evl(u, gd, 11);
  344. t = PIO2 * x2;
  345. c = cos(t);
  346. s = sin(t);
  347. t = PI * x;
  348. cc = 0.5 + (f * s - g * c) / t;
  349. ss = 0.5 - (f * c + g * s) / t;
  350. done:
  351. if (xxa < 0.0) {
  352. cc = -cc;
  353. ss = -ss;
  354. }
  355. *cca = cc;
  356. *ssa = ss;
  357. return (0);
  358. }