/************************************************************************* Cephes Math Library Release 2.8: June, 2000 Copyright by Stephen L. Moshier Contributors: * Sergey Bochkanov (ALGLIB project). Translation from C to pseudocode. See subroutines comments for additional copyrights. >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; namespace alglib { public class fresnel { /************************************************************************* Fresnel integral 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 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ public static void fresnelintegral(double x, ref double c, ref double s) { double xxa = 0; double f = 0; double g = 0; double cc = 0; double ss = 0; double t = 0; double u = 0; double x2 = 0; double sn = 0; double sd = 0; double cn = 0; double cd = 0; double fn = 0; double fd = 0; double gn = 0; double gd = 0; double mpi = 0; double mpio2 = 0; mpi = 3.14159265358979323846; mpio2 = 1.57079632679489661923; xxa = x; x = Math.Abs(xxa); x2 = x*x; if( (double)(x2)<(double)(2.5625) ) { t = x2*x2; sn = -2.99181919401019853726E3; sn = sn*t+7.08840045257738576863E5; sn = sn*t-6.29741486205862506537E7; sn = sn*t+2.54890880573376359104E9; sn = sn*t-4.42979518059697779103E10; sn = sn*t+3.18016297876567817986E11; sd = 1.00000000000000000000E0; sd = sd*t+2.81376268889994315696E2; sd = sd*t+4.55847810806532581675E4; sd = sd*t+5.17343888770096400730E6; sd = sd*t+4.19320245898111231129E8; sd = sd*t+2.24411795645340920940E10; sd = sd*t+6.07366389490084639049E11; cn = -4.98843114573573548651E-8; cn = cn*t+9.50428062829859605134E-6; cn = cn*t-6.45191435683965050962E-4; cn = cn*t+1.88843319396703850064E-2; cn = cn*t-2.05525900955013891793E-1; cn = cn*t+9.99999999999999998822E-1; cd = 3.99982968972495980367E-12; cd = cd*t+9.15439215774657478799E-10; cd = cd*t+1.25001862479598821474E-7; cd = cd*t+1.22262789024179030997E-5; cd = cd*t+8.68029542941784300606E-4; cd = cd*t+4.12142090722199792936E-2; cd = cd*t+1.00000000000000000118E0; s = Math.Sign(xxa)*x*x2*sn/sd; c = Math.Sign(xxa)*x*cn/cd; return; } if( (double)(x)>(double)(36974.0) ) { c = Math.Sign(xxa)*0.5; s = Math.Sign(xxa)*0.5; return; } x2 = x*x; t = mpi*x2; u = 1/(t*t); t = 1/t; fn = 4.21543555043677546506E-1; fn = fn*u+1.43407919780758885261E-1; fn = fn*u+1.15220955073585758835E-2; fn = fn*u+3.45017939782574027900E-4; fn = fn*u+4.63613749287867322088E-6; fn = fn*u+3.05568983790257605827E-8; fn = fn*u+1.02304514164907233465E-10; fn = fn*u+1.72010743268161828879E-13; fn = fn*u+1.34283276233062758925E-16; fn = fn*u+3.76329711269987889006E-20; fd = 1.00000000000000000000E0; fd = fd*u+7.51586398353378947175E-1; fd = fd*u+1.16888925859191382142E-1; fd = fd*u+6.44051526508858611005E-3; fd = fd*u+1.55934409164153020873E-4; fd = fd*u+1.84627567348930545870E-6; fd = fd*u+1.12699224763999035261E-8; fd = fd*u+3.60140029589371370404E-11; fd = fd*u+5.88754533621578410010E-14; fd = fd*u+4.52001434074129701496E-17; fd = fd*u+1.25443237090011264384E-20; gn = 5.04442073643383265887E-1; gn = gn*u+1.97102833525523411709E-1; gn = gn*u+1.87648584092575249293E-2; gn = gn*u+6.84079380915393090172E-4; gn = gn*u+1.15138826111884280931E-5; gn = gn*u+9.82852443688422223854E-8; gn = gn*u+4.45344415861750144738E-10; gn = gn*u+1.08268041139020870318E-12; gn = gn*u+1.37555460633261799868E-15; gn = gn*u+8.36354435630677421531E-19; gn = gn*u+1.86958710162783235106E-22; gd = 1.00000000000000000000E0; gd = gd*u+1.47495759925128324529E0; gd = gd*u+3.37748989120019970451E-1; gd = gd*u+2.53603741420338795122E-2; gd = gd*u+8.14679107184306179049E-4; gd = gd*u+1.27545075667729118702E-5; gd = gd*u+1.04314589657571990585E-7; gd = gd*u+4.60680728146520428211E-10; gd = gd*u+1.10273215066240270757E-12; gd = gd*u+1.38796531259578871258E-15; gd = gd*u+8.39158816283118707363E-19; gd = gd*u+1.86958710162783236342E-22; f = 1-u*fn/fd; g = t*gn/gd; t = mpio2*x2; cc = Math.Cos(t); ss = Math.Sin(t); t = mpi*x; c = 0.5+(f*ss-g*cc)/t; s = 0.5-(f*cc+g*ss)/t; c = c*Math.Sign(xxa); s = s*Math.Sign(xxa); } } }