htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/elljac.c

   1 /* specfunc/elljac.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_pow_int.h>
  26 #include <gsl/gsl_sf_elljac.h>
  27
  28
  29 /* GJ: See [Thompson, Atlas for Computing Mathematical Functions] */
  30
  31 /* BJG 2005-07: New algorithm based on Algorithm 5 from Numerische
  32    Mathematik 7, 78-90 (1965) "Numerical Calculation of Elliptic
  33    Integrals and Elliptic Functions" R. Bulirsch.
  34
  35    Minor tweak is to avoid division by zero when sin(x u_l) = 0 by
  36    computing reflected values sn(K-u) cn(K-u) dn(K-u) and using
  37    transformation from Abramowitz & Stegun table 16.8 column "K-u"*/
  38
  39 int
  40 gsl_sf_elljac_e(double u, double m, double * sn, double * cn, double * dn)
  41 {
  42   if(fabs(m) > 1.0) {
  43     *sn = 0.0;
  44     *cn = 0.0;
  45     *dn = 0.0;
  46     GSL_ERROR ("|m| > 1.0", GSL_EDOM);
  47   }
  48   else if(fabs(m) < 2.0*GSL_DBL_EPSILON) {
  49     *sn = sin(u);
  50     *cn = cos(u);
  51     *dn = 1.0;
  52     return GSL_SUCCESS;
  53   }
  54   else if(fabs(m - 1.0) < 2.0*GSL_DBL_EPSILON) {
  55     *sn = tanh(u);
  56     *cn = 1.0/cosh(u);
  57     *dn = *cn;
  58     return GSL_SUCCESS;
  59   }
  60   else {
  61     int status = GSL_SUCCESS;
  62     const int N = 16;
  63     double mu[16];
  64     double nu[16];
  65     double c[16];
  66     double d[16];
  67     double sin_umu, cos_umu, t, r;
  68     int n = 0;
  69
  70     mu[0] = 1.0;
  71     nu[0] = sqrt(1.0 - m);
  72
  73     while( fabs(mu[n] - nu[n]) > 4.0 * GSL_DBL_EPSILON * fabs(mu[n]+nu[n])) {
  74       mu[n+1] = 0.5 * (mu[n] + nu[n]);
  75       nu[n+1] = sqrt(mu[n] * nu[n]);
  76       ++n;
  77       if(n >= N - 1) {
  78         status = GSL_EMAXITER;
  79         break;
  80       }
  81     }
  82
  83     sin_umu = sin(u * mu[n]);
  84     cos_umu = cos(u * mu[n]);
  85
  86     /* Since sin(u*mu(n)) can be zero we switch to computing sn(K-u),
  87        cn(K-u), dn(K-u) when |sin| < |cos| */
  88
  89     if (fabs(sin_umu) < fabs(cos_umu))
  90       {
  91         t = sin_umu / cos_umu;
  92
  93         c[n] = mu[n] * t;
  94         d[n] = 1.0;
  95
  96         while(n > 0) {
  97           n--;
  98           c[n] = d[n+1] * c[n+1];
  99           r = (c[n+1] * c[n+1]) / mu[n+1];
 100           d[n] = (r + nu[n]) / (r + mu[n]);
 101           }
 102
 103         *dn = sqrt(1.0-m) / d[n];
 104         *cn = (*dn) * GSL_SIGN(cos_umu) / gsl_hypot(1.0, c[n]);
 105         *sn = (*cn) * c[n] /sqrt(1.0-m);
 106       }
 107     else
 108       {
 109         t = cos_umu / sin_umu;
 110
 111         c[n] = mu[n] * t;
 112         d[n] = 1.0;
 113
 114         while(n > 0) {
 115           --n;
 116           c[n] = d[n+1] * c[n+1];
 117           r = (c[n+1] * c[n+1]) / mu[n+1];
 118           d[n] = (r + nu[n]) / (r + mu[n]);
 119         }
 120
 121         *dn = d[n];
 122         *sn = GSL_SIGN(sin_umu) / gsl_hypot(1.0, c[n]);
 123         *cn = c[n] * (*sn);
 124       }
 125
 126     return status;
 127   }
 128 }