coulomb.c

多个常用的特殊函数的例子

/* specfunc/coulomb.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

/* Evaluation of Coulomb wave functions F_L(eta, x), G_L(eta, x),
 * and their derivatives. A combination of Steed's method, asymptotic
 * results, and power series.
 *
 * Steed's method:
 *  [Barnett, CPC 21, 297 (1981)]
 * Power series and other methods:
 *  [Biedenharn et al., PR 97, 542 (1954)]
 *  [Bardin et al., CPC 3, 73 (1972)]
 *  [Abad+Sesma, CPC 71, 110 (1992)]
 */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_coulomb.h>

#include "error.h"

/* the L=0 normalization constant
 * [Abramowitz+Stegun 14.1.8]
 */
static
double
C0sq(double eta)
{
  double twopieta = 2.0*M_PI*eta;

  if(fabs(eta) < GSL_DBL_EPSILON) {
    return 1.0;
  }
  else if(twopieta > GSL_LOG_DBL_MAX) {
    return 0.0;
  }
  else {
    gsl_sf_result scale;
    gsl_sf_expm1_e(twopieta, &scale);
    return twopieta/scale.val;
  }
}


/* the full definition of C_L(eta) for any valid L and eta
 * [Abramowitz and Stegun 14.1.7]
 * This depends on the complex gamma function. For large
 * arguments the phase of the complex gamma function is not
 * very accurately determined. However the modulus is, and that
 * is all that we need to calculate C_L.
 *
 * This is not valid for L <= -3/2  or  L = -1.
 */
static
int
CLeta(double L, double eta, gsl_sf_result * result)
{
  gsl_sf_result ln1; /* log of numerator Gamma function */
  gsl_sf_result ln2; /* log of denominator Gamma function */
  double sgn = 1.0;
  double arg_val, arg_err;

  if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
    gsl_sf_lngamma_e(L+1.0, &ln1);
  }
  else {
    gsl_sf_result p1;                 /* phase of numerator Gamma -- not used */
    gsl_sf_lngamma_complex_e(L+1.0, eta, &ln1, &p1); /* should be ok */
  }

  gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
  if(L < -1.0) sgn = -sgn;

  arg_val  = L*M_LN2 - 0.5*eta*M_PI + ln1.val - ln2.val;
  arg_err  = ln1.err + ln2.err;
  arg_err += GSL_DBL_EPSILON * (fabs(L*M_LN2) + fabs(0.5*eta*M_PI));
  return gsl_sf_exp_err_e(arg_val, arg_err, result);
}


int
gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(lam <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(lam) < GSL_DBL_EPSILON) {
    /* saves a calculation of complex_lngamma(), otherwise not necessary */
    result->val = sqrt(C0sq(eta));
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    return CLeta(lam, eta, result);
  }
}


/* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
 */
int
gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
{
  int k;
  gsl_sf_result cl_0;
  gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
  cl[0] = cl_0.val;

  for(k=1; k<=kmax; k++) {
    double L = lam_min + k;
    cl[k] = cl[k-1] * sqrt(L*L + eta*eta)/(L*(2.0*L+1.0));
  }

  return GSL_SUCCESS;
}


/* Evaluate the series for Phi_L(eta,x) and Phi_L*(eta,x)
 * [Abramowitz+Stegun 14.1.5]
 * [Abramowitz+Stegun 14.1.13]
 *
 * The sequence of coefficients A_k^L is
 * manifestly well-controlled for L >= -1/2
 * and eta < 10.
 *
 * This makes sense since this is the region
 * away from threshold, and you expect
 * the evaluation to become easier as you
 * get farther from threshold.
 *
 * Empirically, this is quite well-behaved for
 *   L >= -1/2
 *   eta < 10
 *   x   < 10
 */
#if 0
static
int
coulomb_Phi_series(const double lam, const double eta, const double x,
                   double * result, double * result_star)
{
  int kmin =   5;
  int kmax = 200;
  int k;
  double Akm2 = 1.0;
  double Akm1 = eta/(lam+1.0);
  double Ak;

  double xpow = x;
  double sum  = Akm2 + Akm1*x;
  double sump = (lam+1.0)*Akm2 + (lam+2.0)*Akm1*x;
  double prev_abs_del   = fabs(Akm1*x);
  double prev_abs_del_p = (lam+2.0) * prev_abs_del;

  for(k=2; k<kmax; k++) {
    double del;
    double del_p;
    double abs_del;
    double abs_del_p;

    Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));

    xpow *= x;
    del   = Ak*xpow;
    del_p = (k+lam+1.0)*del;
    sum  += del;
    sump += del_p;

    abs_del   = fabs(del);
    abs_del_p = fabs(del_p);

    if(          abs_del/(fabs(sum)+abs_del)          < GSL_DBL_EPSILON
       &&   prev_abs_del/(fabs(sum)+prev_abs_del)     < GSL_DBL_EPSILON
       &&      abs_del_p/(fabs(sump)+abs_del_p)       < GSL_DBL_EPSILON
       && prev_abs_del_p/(fabs(sump)+prev_abs_del_p)  < GSL_DBL_EPSILON
       && k > kmin
       ) break;

    /* We need to keep track of the previous delta because when
     * eta is near zero the odd terms of the sum are very small
     * and this could lead to premature termination.
     */
    prev_abs_del   = abs_del;
    prev_abs_del_p = abs_del_p;

    Akm2 = Akm1;
    Akm1 = Ak;
  }

  *result      = sum;
  *result_star = sump;

  if(k==kmax) {
    GSL_ERROR ("error", GSL_EMAXITER);
  }
  else {
    return GSL_SUCCESS;
  }
}
#endif /* 0 */


/* Determine the connection phase, phi_lambda.
 * See coulomb_FG_series() below. We have
 * to be careful about sin(phi)->0. Note that
 * there is an underflow condition for large 
 * positive eta in any case.
 */
static
int
coulomb_connection(const double lam, const double eta,
                   double * cos_phi, double * sin_phi)
{
  if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
    *cos_phi = 1.0;
    *sin_phi = 0.0;
    GSL_ERROR ("error", GSL_EUNDRFLW);
  }
  else if(eta > -GSL_LOG_DBL_EPSILON/(4.0*M_PI)) {
    const double eps = 2.0 * exp(-2.0*M_PI*eta);
    const double tpl = tan(M_PI * lam);
    const double dth = eps * tpl / (tpl*tpl + 1.0);
    *cos_phi = -1.0 + 0.5 * dth*dth;
    *sin_phi = -dth;
    return GSL_SUCCESS;
  }
  else {
    double X   = tanh(M_PI * eta) / tan(M_PI * lam);
    double phi = -atan(X) - (lam + 0.5) * M_PI;
    *cos_phi = cos(phi);
    *sin_phi = sin(phi);
    return GSL_SUCCESS;
  }
}


/* Evaluate the Frobenius series for F_lam(eta,x) and G_lam(eta,x).
 * Homegrown algebra. Evaluates the series for F_{lam} and
 * F_{-lam-1}, then uses
 *    G_{lam} = (F_{lam} cos(phi) - F_{-lam-1}) / sin(phi)
 * where
 *    phi = Arg[Gamma[1+lam+I eta]] - Arg[Gamma[-lam + I eta]] - (lam+1/2)Pi
 *        = Arg[Sin[Pi(-lam+I eta)] - (lam+1/2)Pi
 *        = atan2(-cos(lam Pi)sinh(eta Pi), -sin(lam Pi)cosh(eta Pi)) - (lam+1/2)Pi
 *
 *        = -atan(X) - (lam+1/2) Pi,  X = tanh(eta Pi)/tan(lam Pi)
 *
 * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
 */
static
int
coulomb_FG_series(const double lam, const double eta, const double x,
                  gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  gsl_sf_result ClamA;
  gsl_sf_result ClamB;
  int stat_A = CLeta(lam, eta, &ClamA);
  int stat_B = CLeta(-lam-1.0, eta, &ClamB);
  const double tlp1 = 2.0*lam + 1.0;
  const double pow_x = pow(x, lam);
  double cos_phi_lam;
  double sin_phi_lam;

  double uA_mm2 = 1.0;                  /* uA sum is for F_{lam} */
  double uA_mm1 = x*eta/(lam+1.0);
  double uA_m;
  double uB_mm2 = 1.0;                  /* uB sum is for F_{-lam-1} */
  double uB_mm1 = -x*eta/lam;
  double uB_m;
  double A_sum = uA_mm2 + uA_mm1;
  double B_sum = uB_mm2 + uB_mm1;
  double A_abs_del_prev = fabs(A_sum);
  double B_abs_del_prev = fabs(B_sum);
  gsl_sf_result FA, FB;
  int m = 2;

  int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);

  if(stat_conn == GSL_EUNDRFLW) {
    F->val = 0.0;  /* FIXME: should this be set to Inf too like G? */
    F->err = 0.0;
    OVERFLOW_ERROR(G);
  }

  while(m < max_iter) {
    double abs_dA;
    double abs_dB;
    uA_m = x*(2.0*eta*uA_mm1 - x*uA_mm2)/(m*(m+tlp1));
    uB_m = x*(2.0*eta*uB_mm1 - x*uB_mm2)/(m*(m-tlp1));
    A_sum += uA_m;
    B_sum += uB_m;
    abs_dA = fabs(uA_m);
    abs_dB = fabs(uB_m);
    if(m > 15) {
      /* Don't bother checking until we have gone out a little ways;
       * a minor optimization. Also make sure to check both the
       * current and the previous increment because the odd and even
       * terms of the sum can have very different behaviour, depending
       * on the value of eta.
       */
      double max_abs_dA = GSL_MAX(abs_dA, A_abs_del_prev);
      double max_abs_dB = GSL_MAX(abs_dB, B_abs_del_prev);
      double abs_A = fabs(A_sum);
      double abs_B = fabs(B_sum);
      if(   max_abs_dA/(max_abs_dA + abs_A) < 4.0*GSL_DBL_EPSILON
         && max_abs_dB/(max_abs_dB + abs_B) < 4.0*GSL_DBL_EPSILON
         ) break;
    }
    A_abs_del_prev = abs_dA;
    B_abs_del_prev = abs_dB;
    uA_mm2 = uA_mm1;
    uA_mm1 = uA_m;
    uB_mm2 = uB_mm1;
    uB_mm1 = uB_m;
    m++;
  }

  FA.val = A_sum * ClamA.val * pow_x * x;
  FA.err = fabs(A_sum) * ClamA.err * pow_x * x + 2.0*GSL_DBL_EPSILON*fabs(FA.val);
  FB.val = B_sum * ClamB.val / pow_x;
  FB.err = fabs(B_sum) * ClamB.err / pow_x + 2.0*GSL_DBL_EPSILON*fabs(FB.val);

  F->val = FA.val;
  F->err = FA.err;

  G->val = (FA.val * cos_phi_lam - FB.val)/sin_phi_lam;
  G->err = (FA.err * fabs(cos_phi_lam) + FB.err)/fabs(sin_phi_lam);

  if(m >= max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_ERROR_SELECT_2(stat_A, stat_B);
}


/* Evaluate the Frobenius series for F_0(eta,x) and G_0(eta,x).
 * See [Bardin et al., CPC 3, 73 (1972), (14)-(17)];
 * note the misprint in (17): nu_0=1 is correct, not nu_0=0.
 */
static
int
coulomb_FG0_series(const double eta, const double x,
                   gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  const double x2  = x*x;
  const double tex = 2.0*eta*x;
  gsl_sf_result C0;
  int stat_CL = CLeta(0.0, eta, &C0);
  gsl_sf_result r1pie;
  int psi_stat = gsl_sf_psi_1piy_e(eta, &r1pie);
  double u_mm2 = 0.0;  /* u_0 */
  double u_mm1 = x;    /* u_1 */
  double u_m;
  double v_mm2 = 1.0;                               /* nu_0 */
  double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val);   /* nu_1 */
  double v_m;
  double u_sum = u_mm2 + u_mm1;
  double v_sum = v_mm2 + v_mm1;
  double u_abs_del_prev = fabs(u_sum);
  double v_abs_del_prev = fabs(v_sum);
  int m = 2;
  double u_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(u_sum);
  double v_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(v_sum);
  double ln2x = log(2.0*x);

  while(m < max_iter) {
    double abs_du;
    double abs_dv;
    double m_mm1 = m*(m-1.0);
    u_m = (tex*u_mm1 - x2*u_mm2)/m_mm1;
    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*eta*(2*m-1)*u_m)/m_mm1;
    u_sum += u_m;
    v_sum += v_m;
    abs_du = fabs(u_m);
    abs_dv = fabs(v_m);
    u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
    v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
    if(m > 15) {
      /* Don't bother checking until we have gone out a little ways;
       * a minor optimization. Also make sure to check both the
       * current and the previous increment because the odd and even
       * terms of the sum can have very different behaviour, depending
       * on the value of eta.
       */
      double max_abs_du = GSL_MAX(abs_du, u_abs_del_prev);
      double max_abs_dv = GSL_MAX(abs_dv, v_abs_del_prev);
      double abs_u = fabs(u_sum);
      double abs_v = fabs(v_sum);
      if(   max_abs_du/(max_abs_du + abs_u) < 40.0*GSL_DBL_EPSILON
         && max_abs_dv/(max_abs_dv + abs_v) < 40.0*GSL_DBL_EPSILON
         ) break;
    }
    u_abs_del_prev = abs_du;
    v_abs_del_prev = abs_dv;
    u_mm2 = u_mm1;
    u_mm1 = u_m;
    v_mm2 = v_mm1;
    v_mm1 = v_m;
    m++;
  }

  F->val  = C0.val * u_sum;
  F->err  = C0.err * fabs(u_sum);
  F->err += fabs(C0.val) * u_sum_err;
  F->err += 2.0 * GSL_DBL_EPSILON * fabs(F->val);

  G->val  = (v_sum + 2.0*eta*u_sum * ln2x) / C0.val;
  G->err  = (fabs(v_sum) + fabs(2.0*eta*u_sum * ln2x)) / fabs(C0.val) * fabs(C0.err/C0.val);
  G->err += (v_sum_err + fabs(2.0*eta*u_sum_err*ln2x)) / fabs(C0.val);
  G->err += 2.0 * GSL_DBL_EPSILON * fabs(G->val);

  if(m == max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
}


/* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
 * Homegrown algebra.
 */
static
int
coulomb_FGmhalf_series(const double eta, const double x,
                       gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  const double rx  = sqrt(x);
  const double x2  = x*x;
  const double tex = 2.0*eta*x;
  gsl_sf_result Cmhalf;
  int stat_CL = CLeta(-0.5, eta, &Cmhalf);
  double u_mm2 = 1.0;                      /* u_0 */
  double u_mm1 = tex * u_mm2;              /* u_1 */
  double u_m;
  double v_mm2, v_mm1, v_m;
  double f_sum, g_sum;
  double tmp1;
  gsl_sf_result rpsi_1pe;
  gsl_sf_result rpsi_1p2e;
  int m = 2;

  gsl_sf_psi_1piy_e(eta,     &rpsi_1pe);
  gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);