Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/fresnel.cs @ 5229

Last change on this file since 5229 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)
File size: 7.4 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace alglib
31{
32    public class fresnel
33    {
34        /*************************************************************************
35        Fresnel integral
36
37        Evaluates the Fresnel integrals
38
39                  x
40                  -
41                 | |
42        C(x) =   |   cos(pi/2 t**2) dt,
43               | |
44                -
45                 0
46
47                  x
48                  -
49                 | |
50        S(x) =   |   sin(pi/2 t**2) dt.
51               | |
52                -
53                 0
54
55
56        The integrals are evaluated by a power series for x < 1.
57        For x >= 1 auxiliary functions f(x) and g(x) are employed
58        such that
59
60        C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
61        S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
62
63
64
65        ACCURACY:
66
67         Relative error.
68
69        Arithmetic  function   domain     # trials      peak         rms
70          IEEE       S(x)      0, 10       10000       2.0e-15     3.2e-16
71          IEEE       C(x)      0, 10       10000       1.8e-15     3.3e-16
72
73        Cephes Math Library Release 2.8:  June, 2000
74        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
75        *************************************************************************/
76        public static void fresnelintegral(double x,
77            ref double c,
78            ref double s)
79        {
80            double xxa = 0;
81            double f = 0;
82            double g = 0;
83            double cc = 0;
84            double ss = 0;
85            double t = 0;
86            double u = 0;
87            double x2 = 0;
88            double sn = 0;
89            double sd = 0;
90            double cn = 0;
91            double cd = 0;
92            double fn = 0;
93            double fd = 0;
94            double gn = 0;
95            double gd = 0;
96            double mpi = 0;
97            double mpio2 = 0;
98
99            mpi = 3.14159265358979323846;
100            mpio2 = 1.57079632679489661923;
101            xxa = x;
102            x = Math.Abs(xxa);
103            x2 = x*x;
104            if( (double)(x2)<(double)(2.5625) )
105            {
106                t = x2*x2;
107                sn = -2.99181919401019853726E3;
108                sn = sn*t+7.08840045257738576863E5;
109                sn = sn*t-6.29741486205862506537E7;
110                sn = sn*t+2.54890880573376359104E9;
111                sn = sn*t-4.42979518059697779103E10;
112                sn = sn*t+3.18016297876567817986E11;
113                sd = 1.00000000000000000000E0;
114                sd = sd*t+2.81376268889994315696E2;
115                sd = sd*t+4.55847810806532581675E4;
116                sd = sd*t+5.17343888770096400730E6;
117                sd = sd*t+4.19320245898111231129E8;
118                sd = sd*t+2.24411795645340920940E10;
119                sd = sd*t+6.07366389490084639049E11;
120                cn = -4.98843114573573548651E-8;
121                cn = cn*t+9.50428062829859605134E-6;
122                cn = cn*t-6.45191435683965050962E-4;
123                cn = cn*t+1.88843319396703850064E-2;
124                cn = cn*t-2.05525900955013891793E-1;
125                cn = cn*t+9.99999999999999998822E-1;
126                cd = 3.99982968972495980367E-12;
127                cd = cd*t+9.15439215774657478799E-10;
128                cd = cd*t+1.25001862479598821474E-7;
129                cd = cd*t+1.22262789024179030997E-5;
130                cd = cd*t+8.68029542941784300606E-4;
131                cd = cd*t+4.12142090722199792936E-2;
132                cd = cd*t+1.00000000000000000118E0;
133                s = Math.Sign(xxa)*x*x2*sn/sd;
134                c = Math.Sign(xxa)*x*cn/cd;
135                return;
136            }
137            if( (double)(x)>(double)(36974.0) )
138            {
139                c = Math.Sign(xxa)*0.5;
140                s = Math.Sign(xxa)*0.5;
141                return;
142            }
143            x2 = x*x;
144            t = mpi*x2;
145            u = 1/(t*t);
146            t = 1/t;
147            fn = 4.21543555043677546506E-1;
148            fn = fn*u+1.43407919780758885261E-1;
149            fn = fn*u+1.15220955073585758835E-2;
150            fn = fn*u+3.45017939782574027900E-4;
151            fn = fn*u+4.63613749287867322088E-6;
152            fn = fn*u+3.05568983790257605827E-8;
153            fn = fn*u+1.02304514164907233465E-10;
154            fn = fn*u+1.72010743268161828879E-13;
155            fn = fn*u+1.34283276233062758925E-16;
156            fn = fn*u+3.76329711269987889006E-20;
157            fd = 1.00000000000000000000E0;
158            fd = fd*u+7.51586398353378947175E-1;
159            fd = fd*u+1.16888925859191382142E-1;
160            fd = fd*u+6.44051526508858611005E-3;
161            fd = fd*u+1.55934409164153020873E-4;
162            fd = fd*u+1.84627567348930545870E-6;
163            fd = fd*u+1.12699224763999035261E-8;
164            fd = fd*u+3.60140029589371370404E-11;
165            fd = fd*u+5.88754533621578410010E-14;
166            fd = fd*u+4.52001434074129701496E-17;
167            fd = fd*u+1.25443237090011264384E-20;
168            gn = 5.04442073643383265887E-1;
169            gn = gn*u+1.97102833525523411709E-1;
170            gn = gn*u+1.87648584092575249293E-2;
171            gn = gn*u+6.84079380915393090172E-4;
172            gn = gn*u+1.15138826111884280931E-5;
173            gn = gn*u+9.82852443688422223854E-8;
174            gn = gn*u+4.45344415861750144738E-10;
175            gn = gn*u+1.08268041139020870318E-12;
176            gn = gn*u+1.37555460633261799868E-15;
177            gn = gn*u+8.36354435630677421531E-19;
178            gn = gn*u+1.86958710162783235106E-22;
179            gd = 1.00000000000000000000E0;
180            gd = gd*u+1.47495759925128324529E0;
181            gd = gd*u+3.37748989120019970451E-1;
182            gd = gd*u+2.53603741420338795122E-2;
183            gd = gd*u+8.14679107184306179049E-4;
184            gd = gd*u+1.27545075667729118702E-5;
185            gd = gd*u+1.04314589657571990585E-7;
186            gd = gd*u+4.60680728146520428211E-10;
187            gd = gd*u+1.10273215066240270757E-12;
188            gd = gd*u+1.38796531259578871258E-15;
189            gd = gd*u+8.39158816283118707363E-19;
190            gd = gd*u+1.86958710162783236342E-22;
191            f = 1-u*fn/fd;
192            g = t*gn/gd;
193            t = mpio2*x2;
194            cc = Math.Cos(t);
195            ss = Math.Sin(t);
196            t = mpi*x;
197            c = 0.5+(f*ss-g*cc)/t;
198            s = 0.5-(f*cc+g*ss)/t;
199            c = c*Math.Sign(xxa);
200            s = s*Math.Sign(xxa);
201        }
202    }
203}
Note: See TracBrowser for help on using the repository browser.