mathieu_charv.c

来自「math library from gnu」· C语言 代码 · 共 850 行 · 第 1/2 页

/* specfunc/mathieu_charv.c
 * 
 * Copyright (C) 2002 Lowell Johnson
 * 
 * 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 3 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_mathieu.h>


/* prototypes */
static double solve_cubic(double c2, double c1, double c0);


static double ceer(int order, double qq, double aa, int nterms)
{

  double term, term1;
  int ii, n1;
  

  if (order == 0)
      term = 0.0;
  else
  {      
      term = 2.0*qq*qq/aa;

      if (order != 2)
      {
          n1 = order/2 - 1;

          for (ii=0; ii<n1; ii++)
              term = qq*qq/(aa - 4.0*(ii+1)*(ii+1) - term);
      }
  }
  
  term += order*order;

  term1 = 0.0;

  for (ii=0; ii<nterms; ii++)
      term1 = qq*qq/
        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  if (order == 0)
      term1 *= 2.0;
  
  return (term + term1 - aa);
}


static double ceor(int order, double qq, double aa, int nterms)
{
  double term, term1;
  int ii, n1;

  term = qq;
  n1 = (int)((float)order/2.0 - 0.5);

  for (ii=0; ii<n1; ii++)
      term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
  term += order*order;

  term1 = 0.0;
  for (ii=0; ii<nterms; ii++)
      term1 = qq*qq/
        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);
}


static double seer(int order, double qq, double aa, int nterms)
{
  double term, term1;
  int ii, n1;

  term = 0.0;
  n1 = order/2 - 1;

  for (ii=0; ii<n1; ii++)
      term = qq*qq/(aa - 4.0*(ii + 1)*(ii + 1) - term);
  term += order*order;

  term1 = 0.0;
  for (ii=0; ii<nterms; ii++)
      term1 = qq*qq/
        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);
}


static double seor(int order, double qq, double aa, int nterms)
{
  double term, term1;
  int ii, n1;


  term = -1.0*qq;
  n1 = (int)((float)order/2.0 - 0.5);
  for (ii=0; ii<n1; ii++)
      term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
  term += order*order;

  term1 = 0.0;
  for (ii=0; ii<nterms; ii++)
      term1 = qq*qq/
        (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);

  return (term + term1 - aa);
}


/*----------------------------------------------------------------------------
 * Asymptotic and approximation routines for the characteristic value.
 *
 * Adapted from F.A. Alhargan's paper,
 * "Algorithms for the Computation of All Mathieu Functions of Integer
 * Orders," ACM Transactions on Mathematical Software, Vol. 26, No. 3,
 * September 2000, pp. 390-407.
 *--------------------------------------------------------------------------*/
static double asymptotic(int order, double qq)
{
  double asymp;
  double nn, n2, n4, n6;
  double hh, ah, ah2, ah3, ah4, ah5;


  /* Set up temporary variables to simplify the readability. */
  nn = 2*order + 1;
  n2 = nn*nn;
  n4 = n2*n2;
  n6 = n4*n2;
  
  hh = 2*sqrt(qq);
  ah = 16*hh;
  ah2 = ah*ah;
  ah3 = ah2*ah;
  ah4 = ah3*ah;
  ah5 = ah4*ah;

  /* Equation 38, p. 397 of Alhargan's paper. */
  asymp = -2*qq + nn*hh - 0.125*(n2 + 1);
  asymp -= 0.25*nn*(                          n2 +     3)/ah;
  asymp -= 0.25*   (             5*n4 +    34*n2 +     9)/ah2;
  asymp -= 0.25*nn*(            33*n4 +   410*n2 +   405)/ah3;
  asymp -=         ( 63*n6 +  1260*n4 +  2943*n2 +   486)/ah4;
  asymp -=      nn*(527*n6 + 15617*n4 + 69001*n2 + 41607)/ah5;

  return asymp;
}


/* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
static double solve_cubic(double c2, double c1, double c0)
{
  double qq, rr, ww, ss, tt;

  
  qq = (3*c1 - c2*c2)/9;
  rr = (9*c2*c1 - 27*c0 - 2*c2*c2*c2)/54;
  ww = qq*qq*qq + rr*rr;
  
  if (ww >= 0)
  {
      double t1 = rr + sqrt(ww);
      ss = fabs(t1)/t1*pow(fabs(t1), 1/3.);
      t1 = rr - sqrt(ww);
      tt = fabs(t1)/t1*pow(fabs(t1), 1/3.);
  }
  else
  {
      double theta = acos(rr/sqrt(-qq*qq*qq));
      ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.);
      tt = 0.0;
  }
  
  return (ss + tt - c2/3);
}


/* Compute an initial approximation for the characteristic value. */
static double approx_c(int order, double qq)
{
  double approx;
  double c0, c1, c2;


  if (order < 0)
  {
    GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
  }
  
  switch (order)
  {
      case 0:
          if (qq <= 4)
              return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */
          else
              return asymptotic(order, qq);
          break;

      case 1:
          if (qq <= 4)
              return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */
          else
              return asymptotic(order, qq);
          break;

      case 2:
          if (qq <= 3)
          {
              c2 = -8.0;  /* Eqn. 33 */
              c1 = -48 - 3*qq*qq;
              c0 = 20*qq*qq;
          }
          else
              return asymptotic(order, qq);
          break;

      case 3:
          if (qq <= 6.25)
          {
              c2 = -qq - 8;  /* Eqn. 34 */
              c1 = 16*qq - 128 - 2*qq*qq;
              c0 = qq*qq*(qq + 8);
          }
          else
              return asymptotic(order, qq);
          break;

      default:
          if (order < 70)
          {
              if (1.7*order > 2*sqrt(qq))
              {
                  /* Eqn. 30 */
                  double n2 = (double)(order*order);
                  double n22 = (double)((n2 - 1)*(n2 - 1));
                  double q2 = qq*qq;
                  double q4 = q2*q2;
                  approx = n2 + 0.5*q2/(n2 - 1);
                  approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
                  approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
                      (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
                  if (1.4*order < 2*sqrt(qq))
                  {
                      approx += asymptotic(order, qq);
                      approx *= 0.5;
                  }
              }
              else
                  approx = asymptotic(order, qq);

              return approx;
          }
          else
              return order*order;
  }

  /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
  approx = solve_cubic(c2, c1, c0);
      
  if ( approx < 0 && sqrt(qq) > 0.1*order )
      return asymptotic(order-1, qq);
  else
      return (order*order + fabs(approx));
}

  
static double approx_s(int order, double qq)
{
  double approx;
  double c0, c1, c2;

  
  if (order < 1)
  {
    GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
  }
  
  switch (order)
  {
      case 1:
          if (qq <= 4)
              return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */
          else
              return asymptotic(order-1, qq);
          break;

      case 2:
          if (qq <= 5)
              return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */
          else
              return asymptotic(order-1, qq);
          break;

      case 3:
          if (qq <= 6.25)
          {
              c2 = qq - 8; /* Eqn. 37 */
              c1 = -128 - 16*qq - 2*qq*qq;
              c0 = qq*qq*(8 - qq);
          }
          else
              return asymptotic(order-1, qq);
          break;

      default:
          if (order < 70)
          {
              if (1.7*order > 2*sqrt(qq))
              {
                  /* Eqn. 30 */
                  double n2 = (double)(order*order);
                  double n22 = (double)((n2 - 1)*(n2 - 1));
                  double q2 = qq*qq;
                  double q4 = q2*q2;
                  approx = n2 + 0.5*q2/(n2 - 1);
                  approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
                  approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
                      (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
                  if (1.4*order < 2*sqrt(qq))
                  {
                      approx += asymptotic(order-1, qq);
                      approx *= 0.5;
                  }
              }
              else
                  approx = asymptotic(order-1, qq);

              return approx;
          }
          else
              return order*order;
  }

  /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
  approx = solve_cubic(c2, c1, c0);
      
  if ( approx < 0 && sqrt(qq) > 0.1*order )
      return asymptotic(order-1, qq);
  else
      return (order*order + fabs(approx));
}


int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result *result)
{
  int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
  double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;


  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* If the argument is 0, then the coefficient is simply the square of
     the order. */
  if (qq == 0)
  {
      result->val = order*order;
      result->err = 0.0;
      return GSL_SUCCESS;
  }

  /* Use symmetry characteristics of the functions to handle cases with
     negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */
  if (order < 0)
      order *= -1;
  if (qq < 0.0)
  {
      if (even_odd == 0)
          return gsl_sf_mathieu_a(order, -qq, result);
      else
          return gsl_sf_mathieu_b(order, -qq, result);
  }
  
  /* Compute an initial approximation for the characteristic value. */
  aa = approx_c(order, qq);

  /* Save the original approximation for later comparison. */
  aa_orig = aa;
  
  /* Loop as long as the final value is not near the approximate value
     (with a max limit to avoid potential infinite loop). */
  while (counter < maxcount)
  {
      a1 = aa + 0.001;
      ii = 0;
      if (even_odd == 0)
          fa1 = ceer(order, qq, a1, nterms);
      else
          fa1 = ceor(order, qq, a1, nterms);

      for (;;)
      {
          if (even_odd == 0)
              fa = ceer(order, qq, aa, nterms);
          else
              fa = ceor(order, qq, aa, nterms);
      
          a2 = a1;