123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372 |
- /* fresnl.c
- *
- * Fresnel integral
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, S, C;
- * void fresnl();
- *
- * fresnl( x, _&S, _&C );
- *
- *
- * DESCRIPTION:
- *
- * Evaluates the Fresnel integrals
- *
- * x
- * -
- * | |
- * C(x) = | cos(pi/2 t**2) dt,
- * | |
- * -
- * 0
- *
- * x
- * -
- * | |
- * S(x) = | sin(pi/2 t**2) dt.
- * | |
- * -
- * 0
- *
- *
- * The integrals are evaluated by a power series for x < 1.
- * For x >= 1 auxiliary functions f(x) and g(x) are employed
- * such that
- *
- * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
- * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
- *
- *
- *
- * ACCURACY:
- *
- * Relative error.
- *
- * Arithmetic function domain # trials peak rms
- * IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16
- * IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16
- * DEC S(x) 0, 10 6000 2.2e-16 3.9e-17
- * DEC C(x) 0, 10 5000 2.3e-16 3.9e-17
- */
- /*
- Cephes Math Library Release 2.8: June, 2000
- Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- /* S(x) for small x */
- #ifdef UNK
- static double sn[6] = {
- -2.99181919401019853726E3, 7.08840045257738576863E5,
- -6.29741486205862506537E7, 2.54890880573376359104E9,
- -4.42979518059697779103E10, 3.18016297876567817986E11,
- };
- static double sd[6] = {
- /* 1.00000000000000000000E0,*/
- 2.81376268889994315696E2, 4.55847810806532581675E4,
- 5.17343888770096400730E6, 4.19320245898111231129E8,
- 2.24411795645340920940E10, 6.07366389490084639049E11,
- };
- #endif
- #ifdef DEC
- static unsigned short sn[24] = {
- 0143072, 0176433, 0065455, 0127034, 0045055, 0007200, 0134540, 0026661,
- 0146560, 0035061, 0023667, 0127545, 0050027, 0166503, 0002673, 0153756,
- 0151045, 0002721, 0121737, 0102066, 0051624, 0013177, 0033451, 0021271,
- };
- static unsigned short sd[24] = {
- /*0040200,0000000,0000000,0000000,*/
- 0042214, 0130051, 0112070, 0101617, 0044062, 0010307, 0172346, 0152510,
- 0045635, 0160575, 0143200, 0136642, 0047307, 0171215, 0127457, 0052361,
- 0050647, 0031447, 0032621, 0013510, 0052015, 0064733, 0117362, 0012653,
- };
- #endif
- #ifdef IBMPC
- static unsigned short sn[24] = {
- 0xb5c3, 0x6d65, 0x5fa3, 0xc0a7, 0x05b6, 0x172c, 0xa1d0, 0x4125,
- 0xf5ed, 0x24f6, 0x0746, 0xc18e, 0x7afe, 0x60b7, 0xfda8, 0x41e2,
- 0xf087, 0x347b, 0xa0ba, 0xc224, 0x2457, 0xe6e5, 0x82cf, 0x4252,
- };
- static unsigned short sd[24] = {
- /*0x0000,0x0000,0x0000,0x3ff0,*/
- 0x1072, 0x3287, 0x9605, 0x4071, 0xdaa9, 0xfe9c, 0x4218, 0x40e6,
- 0x17b4, 0xb8d0, 0xbc2f, 0x4153, 0xea9e, 0xb5e5, 0xfe51, 0x41b8,
- 0x22e9, 0xe6b2, 0xe664, 0x4214, 0x42b5, 0x73de, 0xad3b, 0x4261,
- };
- #endif
- #ifdef MIEEE
- static unsigned short sn[24] = {
- 0xc0a7, 0x5fa3, 0x6d65, 0xb5c3, 0x4125, 0xa1d0, 0x172c, 0x05b6,
- 0xc18e, 0x0746, 0x24f6, 0xf5ed, 0x41e2, 0xfda8, 0x60b7, 0x7afe,
- 0xc224, 0xa0ba, 0x347b, 0xf087, 0x4252, 0x82cf, 0xe6e5, 0x2457,
- };
- static unsigned short sd[24] = {
- /*0x3ff0,0x0000,0x0000,0x0000,*/
- 0x4071, 0x9605, 0x3287, 0x1072, 0x40e6, 0x4218, 0xfe9c, 0xdaa9,
- 0x4153, 0xbc2f, 0xb8d0, 0x17b4, 0x41b8, 0xfe51, 0xb5e5, 0xea9e,
- 0x4214, 0xe664, 0xe6b2, 0x22e9, 0x4261, 0xad3b, 0x73de, 0x42b5,
- };
- #endif
- /* C(x) for small x */
- #ifdef UNK
- static double cn[6] = {
- -4.98843114573573548651E-8, 9.50428062829859605134E-6,
- -6.45191435683965050962E-4, 1.88843319396703850064E-2,
- -2.05525900955013891793E-1, 9.99999999999999998822E-1,
- };
- static double cd[7] = {
- 3.99982968972495980367E-12, 9.15439215774657478799E-10,
- 1.25001862479598821474E-7, 1.22262789024179030997E-5,
- 8.68029542941784300606E-4, 4.12142090722199792936E-2,
- 1.00000000000000000118E0,
- };
- #endif
- #ifdef DEC
- static unsigned short cn[24] = {
- 0132126, 0040141, 0063733, 0013231, 0034037, 0072223, 0010200, 0075637,
- 0135451, 0021020, 0073264, 0036057, 0036632, 0131520, 0101316, 0060233,
- 0137522, 0072541, 0136124, 0132202, 0040200, 0000000, 0000000, 0000000,
- };
- static unsigned short cd[28] = {
- 0026614, 0135503, 0051776, 0032631, 0030573, 0121116, 0154033,
- 0126712, 0032406, 0034100, 0012442, 0106212, 0034115, 0017567,
- 0150520, 0164623, 0035543, 0106171, 0177336, 0146351, 0037050,
- 0150073, 0000607, 0171635, 0040200, 0000000, 0000000, 0000000,
- };
- #endif
- #ifdef IBMPC
- static unsigned short cn[24] = {
- 0x62d3, 0x2cfb, 0xc80c, 0xbe6a, 0x0f74, 0x6210, 0xee92, 0x3ee3,
- 0x8786, 0x0ed6, 0x2442, 0xbf45, 0xcc13, 0x1059, 0x566a, 0x3f93,
- 0x9690, 0x378a, 0x4eac, 0xbfca, 0x0000, 0x0000, 0x0000, 0x3ff0,
- };
- static unsigned short cd[28] = {
- 0xc6b3, 0x6a7f, 0x9768, 0x3d91, 0x75b9, 0xdb03, 0x7449,
- 0x3e0f, 0x5191, 0x02a4, 0xc708, 0x3e80, 0x1d32, 0xfa2a,
- 0xa3ee, 0x3ee9, 0xd99d, 0x3fdb, 0x718f, 0x3f4c, 0xfe74,
- 0x6030, 0x1a07, 0x3fa5, 0x0000, 0x0000, 0x0000, 0x3ff0,
- };
- #endif
- #ifdef MIEEE
- static unsigned short cn[24] = {
- 0xbe6a, 0xc80c, 0x2cfb, 0x62d3, 0x3ee3, 0xee92, 0x6210, 0x0f74,
- 0xbf45, 0x2442, 0x0ed6, 0x8786, 0x3f93, 0x566a, 0x1059, 0xcc13,
- 0xbfca, 0x4eac, 0x378a, 0x9690, 0x3ff0, 0x0000, 0x0000, 0x0000,
- };
- static unsigned short cd[28] = {
- 0x3d91, 0x9768, 0x6a7f, 0xc6b3, 0x3e0f, 0x7449, 0xdb03,
- 0x75b9, 0x3e80, 0xc708, 0x02a4, 0x5191, 0x3ee9, 0xa3ee,
- 0xfa2a, 0x1d32, 0x3f4c, 0x718f, 0x3fdb, 0xd99d, 0x3fa5,
- 0x1a07, 0x6030, 0xfe74, 0x3ff0, 0x0000, 0x0000, 0x0000,
- };
- #endif
- /* Auxiliary function f(x) */
- #ifdef UNK
- static double fn[10] = {
- 4.21543555043677546506E-1, 1.43407919780758885261E-1,
- 1.15220955073585758835E-2, 3.45017939782574027900E-4,
- 4.63613749287867322088E-6, 3.05568983790257605827E-8,
- 1.02304514164907233465E-10, 1.72010743268161828879E-13,
- 1.34283276233062758925E-16, 3.76329711269987889006E-20,
- };
- static double fd[10] = {
- /* 1.00000000000000000000E0,*/
- 7.51586398353378947175E-1, 1.16888925859191382142E-1,
- 6.44051526508858611005E-3, 1.55934409164153020873E-4,
- 1.84627567348930545870E-6, 1.12699224763999035261E-8,
- 3.60140029589371370404E-11, 5.88754533621578410010E-14,
- 4.52001434074129701496E-17, 1.25443237090011264384E-20,
- };
- #endif
- #ifdef DEC
- static unsigned short fn[40] = {
- 0037727, 0152216, 0106601, 0016214, 0037422, 0154606, 0112710, 0071355,
- 0036474, 0143453, 0154253, 0166545, 0035264, 0161606, 0022250, 0073743,
- 0033633, 0110036, 0024653, 0136246, 0032003, 0036652, 0041164, 0036413,
- 0027740, 0174122, 0046305, 0036726, 0025501, 0125270, 0121317, 0167667,
- 0023032, 0150555, 0076175, 0047443, 0020061, 0133570, 0070130, 0027657,
- };
- static unsigned short fd[40] = {
- /*0040200,0000000,0000000,0000000,*/
- 0040100, 0063767, 0054413, 0151452, 0037357, 0061566, 0007243, 0065754,
- 0036323, 0005365, 0033552, 0133625, 0035043, 0101123, 0000275, 0165402,
- 0033367, 0146614, 0110623, 0023647, 0031501, 0116644, 0125222, 0144263,
- 0027436, 0062051, 0117235, 0001411, 0025204, 0111543, 0056370, 0036201,
- 0022520, 0071351, 0015227, 0122144, 0017554, 0172240, 0112713, 0005006,
- };
- #endif
- #ifdef IBMPC
- static unsigned short fn[40] = {
- 0x2391, 0xd1b0, 0xfa91, 0x3fda, 0x0e5e, 0xd2b9, 0x5b30, 0x3fc2,
- 0x7dad, 0x7b15, 0x98e5, 0x3f87, 0x0efc, 0xc495, 0x9c70, 0x3f36,
- 0x7795, 0xc535, 0x7203, 0x3ed3, 0x87a1, 0x484e, 0x67b5, 0x3e60,
- 0xa7bb, 0x4998, 0x1f0a, 0x3ddc, 0xfdf7, 0x1459, 0x3557, 0x3d48,
- 0xa9e4, 0xaf8f, 0x5a2d, 0x3ca3, 0x05f6, 0x0e0b, 0x36ef, 0x3be6,
- };
- static unsigned short fd[40] = {
- /*0x0000,0x0000,0x0000,0x3ff0,*/
- 0x7a65, 0xeb21, 0x0cfe, 0x3fe8, 0x6d7d, 0xc1d4, 0xec6e, 0x3fbd,
- 0x56f3, 0xa6ed, 0x615e, 0x3f7a, 0xbd60, 0x6017, 0x704a, 0x3f24,
- 0x64f5, 0x9232, 0xf9b1, 0x3ebe, 0x5916, 0x9552, 0x33b4, 0x3e48,
- 0xa061, 0x33d3, 0xcc85, 0x3dc3, 0x0790, 0x6b9f, 0x926c, 0x3d30,
- 0xf48d, 0x2352, 0x0e5d, 0x3c8a, 0x6141, 0x12b9, 0x9e94, 0x3bcd,
- };
- #endif
- #ifdef MIEEE
- static unsigned short fn[40] = {
- 0x3fda, 0xfa91, 0xd1b0, 0x2391, 0x3fc2, 0x5b30, 0xd2b9, 0x0e5e,
- 0x3f87, 0x98e5, 0x7b15, 0x7dad, 0x3f36, 0x9c70, 0xc495, 0x0efc,
- 0x3ed3, 0x7203, 0xc535, 0x7795, 0x3e60, 0x67b5, 0x484e, 0x87a1,
- 0x3ddc, 0x1f0a, 0x4998, 0xa7bb, 0x3d48, 0x3557, 0x1459, 0xfdf7,
- 0x3ca3, 0x5a2d, 0xaf8f, 0xa9e4, 0x3be6, 0x36ef, 0x0e0b, 0x05f6,
- };
- static unsigned short fd[40] = {
- /*0x3ff0,0x0000,0x0000,0x0000,*/
- 0x3fe8, 0x0cfe, 0xeb21, 0x7a65, 0x3fbd, 0xec6e, 0xc1d4, 0x6d7d,
- 0x3f7a, 0x615e, 0xa6ed, 0x56f3, 0x3f24, 0x704a, 0x6017, 0xbd60,
- 0x3ebe, 0xf9b1, 0x9232, 0x64f5, 0x3e48, 0x33b4, 0x9552, 0x5916,
- 0x3dc3, 0xcc85, 0x33d3, 0xa061, 0x3d30, 0x926c, 0x6b9f, 0x0790,
- 0x3c8a, 0x0e5d, 0x2352, 0xf48d, 0x3bcd, 0x9e94, 0x12b9, 0x6141,
- };
- #endif
- /* Auxiliary function g(x) */
- #ifdef UNK
- static double gn[11] = {
- 5.04442073643383265887E-1, 1.97102833525523411709E-1,
- 1.87648584092575249293E-2, 6.84079380915393090172E-4,
- 1.15138826111884280931E-5, 9.82852443688422223854E-8,
- 4.45344415861750144738E-10, 1.08268041139020870318E-12,
- 1.37555460633261799868E-15, 8.36354435630677421531E-19,
- 1.86958710162783235106E-22,
- };
- static double gd[11] = {
- /* 1.00000000000000000000E0,*/
- 1.47495759925128324529E0, 3.37748989120019970451E-1,
- 2.53603741420338795122E-2, 8.14679107184306179049E-4,
- 1.27545075667729118702E-5, 1.04314589657571990585E-7,
- 4.60680728146520428211E-10, 1.10273215066240270757E-12,
- 1.38796531259578871258E-15, 8.39158816283118707363E-19,
- 1.86958710162783236342E-22,
- };
- #endif
- #ifdef DEC
- static unsigned short gn[44] = {
- 0040001, 0021435, 0120406, 0053123, 0037511, 0152523, 0037703, 0122011,
- 0036631, 0134302, 0122721, 0110235, 0035463, 0051712, 0043215, 0114732,
- 0034101, 0025677, 0147725, 0057630, 0032323, 0010342, 0067523, 0002206,
- 0030364, 0152247, 0110007, 0054107, 0026230, 0057654, 0035464, 0047124,
- 0023706, 0036401, 0167705, 0045440, 0021166, 0154447, 0105632, 0142461,
- 0016142, 0002353, 0011175, 0170530,
- };
- static unsigned short gd[44] = {
- /*0040200,0000000,0000000,0000000,*/
- 0040274, 0145551, 0016742, 0127005, 0037654, 0166557, 0076416, 0015165,
- 0036717, 0140217, 0030675, 0050111, 0035525, 0110060, 0076405, 0070502,
- 0034125, 0176061, 0060120, 0031730, 0032340, 0001615, 0054343, 0120501,
- 0030375, 0041414, 0070747, 0107060, 0026233, 0031034, 0160757, 0074526,
- 0023710, 0003341, 0137100, 0144664, 0021167, 0126414, 0023774, 0015435,
- 0016142, 0002353, 0011175, 0170530,
- };
- #endif
- #ifdef IBMPC
- static unsigned short gn[44] = {
- 0xcaca, 0xb420, 0x2463, 0x3fe0, 0x7481, 0x67f8, 0x3aaa, 0x3fc9, 0x3214,
- 0x54ba, 0x3718, 0x3f93, 0xb33b, 0x48d1, 0x6a79, 0x3f46, 0xabf3, 0xf9fa,
- 0x2577, 0x3ee8, 0x6091, 0x4dea, 0x621c, 0x3e7a, 0xeb09, 0xf200, 0x9a94,
- 0x3dfe, 0x89cb, 0x8766, 0x0bf5, 0x3d73, 0xa964, 0x3df8, 0xc7a0, 0x3cd8,
- 0x58a6, 0xf173, 0xdb24, 0x3c2e, 0xbe2b, 0x624f, 0x409d, 0x3b6c,
- };
- static unsigned short gd[44] = {
- /*0x0000,0x0000,0x0000,0x3ff0,*/
- 0x55c1, 0x23bc, 0x996d, 0x3ff7, 0xc34f, 0xefa1, 0x9dad, 0x3fd5, 0xaa09,
- 0xe637, 0xf811, 0x3f99, 0xae28, 0x0fa0, 0xb206, 0x3f4a, 0x067b, 0x2c0a,
- 0xbf86, 0x3eea, 0x7428, 0xab1c, 0x0071, 0x3e7c, 0xf1c6, 0x8e3c, 0xa861,
- 0x3dff, 0xef2b, 0x9c3d, 0x6643, 0x3d73, 0x1936, 0x37c8, 0x00dc, 0x3cd9,
- 0x8364, 0x84ff, 0xf5a1, 0x3c2e, 0xbe2b, 0x624f, 0x409d, 0x3b6c,
- };
- #endif
- #ifdef MIEEE
- static unsigned short gn[44] = {
- 0x3fe0, 0x2463, 0xb420, 0xcaca, 0x3fc9, 0x3aaa, 0x67f8, 0x7481, 0x3f93,
- 0x3718, 0x54ba, 0x3214, 0x3f46, 0x6a79, 0x48d1, 0xb33b, 0x3ee8, 0x2577,
- 0xf9fa, 0xabf3, 0x3e7a, 0x621c, 0x4dea, 0x6091, 0x3dfe, 0x9a94, 0xf200,
- 0xeb09, 0x3d73, 0x0bf5, 0x8766, 0x89cb, 0x3cd8, 0xc7a0, 0x3df8, 0xa964,
- 0x3c2e, 0xdb24, 0xf173, 0x58a6, 0x3b6c, 0x409d, 0x624f, 0xbe2b,
- };
- static unsigned short gd[44] = {
- /*0x3ff0,0x0000,0x0000,0x0000,*/
- 0x3ff7, 0x996d, 0x23bc, 0x55c1, 0x3fd5, 0x9dad, 0xefa1, 0xc34f, 0x3f99,
- 0xf811, 0xe637, 0xaa09, 0x3f4a, 0xb206, 0x0fa0, 0xae28, 0x3eea, 0xbf86,
- 0x2c0a, 0x067b, 0x3e7c, 0x0071, 0xab1c, 0x7428, 0x3dff, 0xa861, 0x8e3c,
- 0xf1c6, 0x3d73, 0x6643, 0x9c3d, 0xef2b, 0x3cd9, 0x00dc, 0x37c8, 0x1936,
- 0x3c2e, 0xf5a1, 0x84ff, 0x8364, 0x3b6c, 0x409d, 0x624f, 0xbe2b,
- };
- #endif
- #ifdef ANSIPROT
- extern double fabs(double);
- extern double cos(double);
- extern double sin(double);
- extern double polevl(double, void *, int);
- extern double p1evl(double, void *, int);
- #else
- double fabs(), cos(), sin(), polevl(), p1evl();
- #endif
- extern double PI, PIO2, MACHEP;
- int fresnl(xxa, ssa, cca) double xxa, *ssa, *cca;
- {
- double f, g, cc, ss, c, s, t, u;
- double x, x2;
- x = fabs(xxa);
- x2 = x * x;
- if (x2 < 2.5625) {
- t = x2 * x2;
- ss = x * x2 * polevl(t, sn, 5) / p1evl(t, sd, 6);
- cc = x * polevl(t, cn, 5) / polevl(t, cd, 6);
- goto done;
- }
- if (x > 36974.0) {
- cc = 0.5;
- ss = 0.5;
- goto done;
- }
- /* Asymptotic power series auxiliary functions
- * for large argument
- */
- x2 = x * x;
- t = PI * x2;
- u = 1.0 / (t * t);
- t = 1.0 / t;
- f = 1.0 - u * polevl(u, fn, 9) / p1evl(u, fd, 10);
- g = t * polevl(u, gn, 10) / p1evl(u, gd, 11);
- t = PIO2 * x2;
- c = cos(t);
- s = sin(t);
- t = PI * x;
- cc = 0.5 + (f * s - g * c) / t;
- ss = 0.5 - (f * c + g * s) / t;
- done:
- if (xxa < 0.0) {
- cc = -cc;
- ss = -ss;
- }
- *cca = cc;
- *ssa = ss;
- return (0);
- }
|