Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/jacobianelliptic.cs @ 2575

Last change on this file since 2575 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 5.4 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright 1984, 1987, 2000 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 jacobianelliptic
33    {
34        /*************************************************************************
35        Jacobian Elliptic Functions
36
37        Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
38        and dn(u|m) of parameter m between 0 and 1, and real
39        argument u.
40
41        These functions are periodic, with quarter-period on the
42        real axis equal to the complete elliptic integral
43        ellpk(1.0-m).
44
45        Relation to incomplete elliptic integral:
46        If u = ellik(phi,m), then sn(u|m) = sin(phi),
47        and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
48
49        Computation is by means of the arithmetic-geometric mean
50        algorithm, except when m is within 1e-9 of 0 or 1.  In the
51        latter case with m close to 1, the approximation applies
52        only for phi < pi/2.
53
54        ACCURACY:
55
56        Tested at random points with u between 0 and 10, m between
57        0 and 1.
58
59                   Absolute error (* = relative error):
60        arithmetic   function   # trials      peak         rms
61           IEEE      phi         10000       9.2e-16*    1.4e-16*
62           IEEE      sn          50000       4.1e-15     4.6e-16
63           IEEE      cn          40000       3.6e-15     4.4e-16
64           IEEE      dn          10000       1.3e-12     1.8e-14
65
66         Peak error observed in consistency check using addition
67        theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
68        the above relation to the incomplete elliptic integral.
69        Accuracy deteriorates when u is large.
70
71        Cephes Math Library Release 2.8:  June, 2000
72        Copyright 1984, 1987, 2000 by Stephen L. Moshier
73        *************************************************************************/
74        public static void jacobianellipticfunctions(double u,
75            double m,
76            ref double sn,
77            ref double cn,
78            ref double dn,
79            ref double ph)
80        {
81            double ai = 0;
82            double b = 0;
83            double phi = 0;
84            double t = 0;
85            double twon = 0;
86            double[] a = new double[0];
87            double[] c = new double[0];
88            int i = 0;
89
90            System.Diagnostics.Debug.Assert((double)(m)>=(double)(0) & (double)(m)<=(double)(1), "Domain error in JacobianEllipticFunctions: m<0 or m>1");
91            a = new double[8+1];
92            c = new double[8+1];
93            if( (double)(m)<(double)(1.0e-9) )
94            {
95                t = Math.Sin(u);
96                b = Math.Cos(u);
97                ai = 0.25*m*(u-t*b);
98                sn = t-ai*b;
99                cn = b+ai*t;
100                ph = u-ai;
101                dn = 1.0-0.5*m*t*t;
102                return;
103            }
104            if( (double)(m)>=(double)(0.9999999999) )
105            {
106                ai = 0.25*(1.0-m);
107                b = Math.Cosh(u);
108                t = Math.Tanh(u);
109                phi = 1.0/b;
110                twon = b*Math.Sinh(u);
111                sn = t+ai*(twon-u)/(b*b);
112                ph = 2.0*Math.Atan(Math.Exp(u))-1.57079632679489661923+ai*(twon-u)/b;
113                ai = ai*t*phi;
114                cn = phi-ai*(twon-u);
115                dn = phi+ai*(twon+u);
116                return;
117            }
118            a[0] = 1.0;
119            b = Math.Sqrt(1.0-m);
120            c[0] = Math.Sqrt(m);
121            twon = 1.0;
122            i = 0;
123            while( (double)(Math.Abs(c[i]/a[i]))>(double)(AP.Math.MachineEpsilon) )
124            {
125                if( i>7 )
126                {
127                    System.Diagnostics.Debug.Assert(false, "Overflow in JacobianEllipticFunctions");
128                    break;
129                }
130                ai = a[i];
131                i = i+1;
132                c[i] = 0.5*(ai-b);
133                t = Math.Sqrt(ai*b);
134                a[i] = 0.5*(ai+b);
135                b = t;
136                twon = twon*2.0;
137            }
138            phi = twon*a[i]*u;
139            do
140            {
141                t = c[i]*Math.Sin(phi)/a[i];
142                b = phi;
143                phi = (Math.Asin(t)+phi)/2.0;
144                i = i-1;
145            }
146            while( i!=0 );
147            sn = Math.Sin(phi);
148            t = Math.Cos(phi);
149            cn = t;
150            dn = t/Math.Cos(phi-b);
151            ph = phi;
152        }
153    }
154}
Note: See TracBrowser for help on using the repository browser.