| /* specfunc/elljac.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman |
| * |
| * 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; 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. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| /* Author: G. Jungman */ |
| |
| #include <config.h> |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_errno.h> |
| #include <gsl/gsl_sf_pow_int.h> |
| #include <gsl/gsl_sf_elljac.h> |
| |
| |
| /* GJ: See [Thompson, Atlas for Computing Mathematical Functions] */ |
| |
| /* BJG 2005-07: New algorithm based on Algorithm 5 from Numerische |
| Mathematik 7, 78-90 (1965) "Numerical Calculation of Elliptic |
| Integrals and Elliptic Functions" R. Bulirsch. |
| |
| Minor tweak is to avoid division by zero when sin(x u_l) = 0 by |
| computing reflected values sn(K-u) cn(K-u) dn(K-u) and using |
| transformation from Abramowitz & Stegun table 16.8 column "K-u"*/ |
| |
| int |
| gsl_sf_elljac_e(double u, double m, double * sn, double * cn, double * dn) |
| { |
| if(fabs(m) > 1.0) { |
| *sn = 0.0; |
| *cn = 0.0; |
| *dn = 0.0; |
| GSL_ERROR ("|m| > 1.0", GSL_EDOM); |
| } |
| else if(fabs(m) < 2.0*GSL_DBL_EPSILON) { |
| *sn = sin(u); |
| *cn = cos(u); |
| *dn = 1.0; |
| return GSL_SUCCESS; |
| } |
| else if(fabs(m - 1.0) < 2.0*GSL_DBL_EPSILON) { |
| *sn = tanh(u); |
| *cn = 1.0/cosh(u); |
| *dn = *cn; |
| return GSL_SUCCESS; |
| } |
| else { |
| int status = GSL_SUCCESS; |
| const int N = 16; |
| double mu[16]; |
| double nu[16]; |
| double c[16]; |
| double d[16]; |
| double sin_umu, cos_umu, t, r; |
| int n = 0; |
| |
| mu[0] = 1.0; |
| nu[0] = sqrt(1.0 - m); |
| |
| while( fabs(mu[n] - nu[n]) > 4.0 * GSL_DBL_EPSILON * fabs(mu[n]+nu[n])) { |
| mu[n+1] = 0.5 * (mu[n] + nu[n]); |
| nu[n+1] = sqrt(mu[n] * nu[n]); |
| ++n; |
| if(n >= N - 1) { |
| status = GSL_EMAXITER; |
| break; |
| } |
| } |
| |
| sin_umu = sin(u * mu[n]); |
| cos_umu = cos(u * mu[n]); |
| |
| /* Since sin(u*mu(n)) can be zero we switch to computing sn(K-u), |
| cn(K-u), dn(K-u) when |sin| < |cos| */ |
| |
| if (fabs(sin_umu) < fabs(cos_umu)) |
| { |
| t = sin_umu / cos_umu; |
| |
| c[n] = mu[n] * t; |
| d[n] = 1.0; |
| |
| while(n > 0) { |
| n--; |
| c[n] = d[n+1] * c[n+1]; |
| r = (c[n+1] * c[n+1]) / mu[n+1]; |
| d[n] = (r + nu[n]) / (r + mu[n]); |
| } |
| |
| *dn = sqrt(1.0-m) / d[n]; |
| *cn = (*dn) * GSL_SIGN(cos_umu) / gsl_hypot(1.0, c[n]); |
| *sn = (*cn) * c[n] /sqrt(1.0-m); |
| } |
| else |
| { |
| t = cos_umu / sin_umu; |
| |
| c[n] = mu[n] * t; |
| d[n] = 1.0; |
| |
| while(n > 0) { |
| --n; |
| c[n] = d[n+1] * c[n+1]; |
| r = (c[n+1] * c[n+1]) / mu[n+1]; |
| d[n] = (r + nu[n]) / (r + mu[n]); |
| } |
| |
| *dn = d[n]; |
| *sn = GSL_SIGN(sin_umu) / gsl_hypot(1.0, c[n]); |
| *cn = c[n] * (*sn); |
| } |
| |
| return status; |
| } |
| } |