← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/ellpj.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/ellpj.c
 
1
/*                                                      ellpj.c
 
2
 *
 
3
 *      Jacobian Elliptic Functions
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double u, m, sn, cn, dn, phi;
 
10
 * int ellpj();
 
11
 *
 
12
 * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
 
13
 *
 
14
 *
 
15
 *
 
16
 * DESCRIPTION:
 
17
 *
 
18
 *
 
19
 * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
 
20
 * and dn(u|m) of parameter m between 0 and 1, and real
 
21
 * argument u.
 
22
 *
 
23
 * These functions are periodic, with quarter-period on the
 
24
 * real axis equal to the complete elliptic integral
 
25
 * ellpk(m).
 
26
 *
 
27
 * Relation to incomplete elliptic integral:
 
28
 * If u = ellik(phi,m), then sn(u|m) = sin(phi),
 
29
 * and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
 
30
 *
 
31
 * Computation is by means of the arithmetic-geometric mean
 
32
 * algorithm, except when m is within 1e-9 of 0 or 1.  In the
 
33
 * latter case with m close to 1, the approximation applies
 
34
 * only for phi < pi/2.
 
35
 *
 
36
 * ACCURACY:
 
37
 *
 
38
 * Tested at random points with u between 0 and 10, m between
 
39
 * 0 and 1.
 
40
 *
 
41
 *            Absolute error (* = relative error):
 
42
 * arithmetic   function   # trials      peak         rms
 
43
 *    DEC       sn           1800       4.5e-16     8.7e-17
 
44
 *    IEEE      phi         10000       9.2e-16*    1.4e-16*
 
45
 *    IEEE      sn          50000       4.1e-15     4.6e-16
 
46
 *    IEEE      cn          40000       3.6e-15     4.4e-16
 
47
 *    IEEE      dn          10000       1.3e-12     1.8e-14
 
48
 *
 
49
 *  Peak error observed in consistency check using addition
 
50
 * theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
 
51
 * the above relation to the incomplete elliptic integral.
 
52
 * Accuracy deteriorates when u is large.
 
53
 *
 
54
 */
 
55
 
 
56
/*                                                      ellpj.c         */
 
57
 
 
58
 
 
59
/*
 
60
Cephes Math Library Release 2.0:  April, 1987
 
61
Copyright 1984, 1987 by Stephen L. Moshier
 
62
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
63
*/
 
64
 
 
65
#include "mconf.h"
 
66
#ifndef ANSIPROT
 
67
double sqrt(), fabs(), sin(), cos(), asin(), tanh();
 
68
double sinh(), cosh(), atan(), exp();
 
69
#endif
 
70
extern double PIO2, MACHEP, NAN;
 
71
 
 
72
int ellpj( u, m, sn, cn, dn, ph )
 
73
double u, m;
 
74
double *sn, *cn, *dn, *ph;
 
75
{
 
76
double ai, b, phi, t, twon;
 
77
double a[9], c[9];
 
78
int i;
 
79
 
 
80
 
 
81
/* Check for special cases */
 
82
 
 
83
if( m < 0.0 || m > 1.0 )
 
84
        {
 
85
        mtherr( "ellpj", DOMAIN );
 
86
        *sn = NAN;
 
87
        *cn = NAN;
 
88
        *ph = NAN;
 
89
        *dn = NAN;
 
90
        return(-1);
 
91
        }
 
92
if( m < 1.0e-9 )
 
93
        {
 
94
        t = sin(u);
 
95
        b = cos(u);
 
96
        ai = 0.25 * m * (u - t*b);
 
97
        *sn = t - ai*b;
 
98
        *cn = b + ai*t;
 
99
        *ph = u - ai;
 
100
        *dn = 1.0 - 0.5*m*t*t;
 
101
        return(0);
 
102
        }
 
103
 
 
104
if( m >= 0.9999999999 )
 
105
        {
 
106
        ai = 0.25 * (1.0-m);
 
107
        b = cosh(u);
 
108
        t = tanh(u);
 
109
        phi = 1.0/b;
 
110
        twon = b * sinh(u);
 
111
        *sn = t + ai * (twon - u)/(b*b);
 
112
        *ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;
 
113
        ai *= t * phi;
 
114
        *cn = phi - ai * (twon - u);
 
115
        *dn = phi + ai * (twon + u);
 
116
        return(0);
 
117
        }
 
118
 
 
119
 
 
120
/*      A. G. M. scale          */
 
121
a[0] = 1.0;
 
122
b = sqrt(1.0 - m);
 
123
c[0] = sqrt(m);
 
124
twon = 1.0;
 
125
i = 0;
 
126
 
 
127
while( fabs(c[i]/a[i]) > MACHEP )
 
128
        {
 
129
        if( i > 7 )
 
130
                {
 
131
                mtherr( "ellpj", OVERFLOW );
 
132
                goto done;
 
133
                }
 
134
        ai = a[i];
 
135
        ++i;
 
136
        c[i] = ( ai - b )/2.0;
 
137
        t = sqrt( ai * b );
 
138
        a[i] = ( ai + b )/2.0;
 
139
        b = t;
 
140
        twon *= 2.0;
 
141
        }
 
142
 
 
143
done:
 
144
 
 
145
/* backward recurrence */
 
146
phi = twon * a[i] * u;
 
147
do
 
148
        {
 
149
        t = c[i] * sin(phi) / a[i];
 
150
        b = phi;
 
151
        phi = (asin(t) + phi)/2.0;
 
152
        }
 
153
while( --i );
 
154
 
 
155
*sn = sin(phi);
 
156
t = cos(phi);
 
157
*cn = t;
 
158
*dn = t/cos(phi-b);
 
159
*ph = phi;
 
160
return(0);
 
161
}