ellpj.c
Go to the documentation of this file.
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  * IEEE phi 10000 9.2e-16* 1.4e-16*
44  * IEEE sn 50000 4.1e-15 4.6e-16
45  * IEEE cn 40000 3.6e-15 4.4e-16
46  * IEEE dn 10000 1.3e-12 1.8e-14
47  *
48  * Peak error observed in consistency check using addition
49  * theorem for sn(u+v) was 4e-16 (absolute). Also tested by
50  * the above relation to the incomplete elliptic integral.
51  * Accuracy deteriorates when u is large.
52  *
53  */
54 
55 /* ellpj.c */
56 
57 
58 /*
59  * Cephes Math Library Release 2.0: April, 1987
60  * Copyright 1984, 1987 by Stephen L. Moshier
61  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
62  */
63 
64 /* Scipy changes:
65  * - 07-18-2016: improve evaluation of dn near quarter periods
66  */
67 
68 #include "mconf.h"
69 extern double MACHEP;
70 
71 int ellpj(double u, double m, double *sn, double *cn, double *dn, double *ph)
72 {
73  double ai, b, phi, t, twon, dnfac;
74  double a[9], c[9];
75  int i;
76 
77  /* Check for special cases */
78  if (m < 0.0 || m > 1.0 || cephes_isnan(m)) {
79  sf_error("ellpj", SF_ERROR_DOMAIN, NULL);
80  *sn = NAN;
81  *cn = NAN;
82  *ph = NAN;
83  *dn = NAN;
84  return (-1);
85  }
86  if (m < 1.0e-9) {
87  t = sin(u);
88  b = cos(u);
89  ai = 0.25 * m * (u - t * b);
90  *sn = t - ai * b;
91  *cn = b + ai * t;
92  *ph = u - ai;
93  *dn = 1.0 - 0.5 * m * t * t;
94  return (0);
95  }
96  if (m >= 0.9999999999) {
97  ai = 0.25 * (1.0 - m);
98  b = cosh(u);
99  t = tanh(u);
100  phi = 1.0 / b;
101  twon = b * sinh(u);
102  *sn = t + ai * (twon - u) / (b * b);
103  *ph = 2.0 * atan(exp(u)) - M_PI_2 + ai * (twon - u) / b;
104  ai *= t * phi;
105  *cn = phi - ai * (twon - u);
106  *dn = phi + ai * (twon + u);
107  return (0);
108  }
109 
110  /* A. G. M. scale. See DLMF 22.20(ii) */
111  a[0] = 1.0;
112  b = sqrt(1.0 - m);
113  c[0] = sqrt(m);
114  twon = 1.0;
115  i = 0;
116 
117  while (fabs(c[i] / a[i]) > MACHEP) {
118  if (i > 7) {
119  sf_error("ellpj", SF_ERROR_OVERFLOW, NULL);
120  goto done;
121  }
122  ai = a[i];
123  ++i;
124  c[i] = (ai - b) / 2.0;
125  t = sqrt(ai * b);
126  a[i] = (ai + b) / 2.0;
127  b = t;
128  twon *= 2.0;
129  }
130 
131  done:
132  /* backward recurrence */
133  phi = twon * a[i] * u;
134  do {
135  t = c[i] * sin(phi) / a[i];
136  b = phi;
137  phi = (asin(t) + phi) / 2.0;
138  }
139  while (--i);
140 
141  *sn = sin(phi);
142  t = cos(phi);
143  *cn = t;
144  dnfac = cos(phi - b);
145  /* See discussion after DLMF 22.20.5 */
146  if (fabs(dnfac) < 0.1) {
147  *dn = sqrt(1 - m*(*sn)*(*sn));
148  }
149  else {
150  *dn = t / dnfac;
151  }
152  *ph = phi;
153  return (0);
154 }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
atan
const EIGEN_DEVICE_FUNC AtanReturnType atan() const
Definition: ArrayCwiseUnaryOps.h:283
cephes_isnan
#define cephes_isnan(x)
Definition: mconf.h:99
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
asin
const EIGEN_DEVICE_FUNC AsinReturnType asin() const
Definition: ArrayCwiseUnaryOps.h:311
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
ellpj
int ellpj(double u, double m, double *sn, double *cn, double *dn, double *ph)
Definition: ellpj.c:71
ceres::cos
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
cosh
const EIGEN_DEVICE_FUNC CoshReturnType cosh() const
Definition: ArrayCwiseUnaryOps.h:353
cn
static double cn[6]
Definition: fresnl.c:83
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
MACHEP
double MACHEP
Definition: const.c:54
M_PI_2
#define M_PI_2
Definition: mconf.h:118
sn
static double sn[6]
Definition: fresnl.c:63
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
tanh
const EIGEN_DEVICE_FUNC TanhReturnType tanh() const
Definition: ArrayCwiseUnaryOps.h:325
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
NULL
#define NULL
Definition: ccolamd.c:609
sinh
const EIGEN_DEVICE_FUNC SinhReturnType sinh() const
Definition: ArrayCwiseUnaryOps.h:339
align_3::t
Point2 t(10, 10)
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:29