/* 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 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 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] * hypot(L, 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

    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);

  v_mm2 = 2.0*M_EULER - M_LN2 - rpsi_1pe.val + 2.0*rpsi_1p2e.val;

  v_mm1 = tex*(v_mm2 - 2.0*u_mm2);

  f_sum = u_mm2 + u_mm1;

  g_sum = v_mm2 + v_mm1;

  while(m < max_iter) {

    double m2 = m*m;

    u_m = (tex*u_mm1 - x2*u_mm2)/m2;

    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*m*u_m)/m2;

    f_sum += u_m;

    g_sum += v_m;

    if(   f_sum != 0.0

       && g_sum != 0.0

       && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;

    u_mm2 = u_mm1;

    u_mm1 = u_m;

    v_mm2 = v_mm1;

    v_mm1 = v_m;

    m++;

  }

  F->val = Cmhalf.val * rx * f_sum;

  F->err = Cmhalf.err * fabs(rx * f_sum) + 2.0*GSL_DBL_EPSILON*fabs(F->val);

  tmp1 = f_sum*log(x);

  G->val = -rx*(tmp1 + g_sum)/Cmhalf.val;

  G->err = fabs(rx)*(fabs(tmp1) + fabs(g_sum))/fabs(Cmhalf.val) * fabs(Cmhalf.err/Cmhalf.val);

  if(m == max_iter)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return stat_CL;

}

/* Evolve the backwards recurrence for F,F'.

 *

 *    F_{lam-1}  = (S_lam F_lam + F_lam') / R_lam

 *    F_{lam-1}' = (S_lam F_{lam-1} - R_lam F_lam)

 * where

 *    R_lam = sqrt(1 + (eta/lam)^2)

 *    S_lam = lam/x + eta/lam

 *

 */

static

int

coulomb_F_recur(double lam_min, int kmax,

                double eta, double x,

                double F_lam_max, double Fp_lam_max,

                double * F_lam_min, double * Fp_lam_min

                )

{

  double x_inv = 1.0/x;

  double fcl = F_lam_max;

  double fpl = Fp_lam_max;

  double lam_max = lam_min + kmax;

  double lam = lam_max;

  int k;

  for(k=kmax-1; k>=0; k--) {

    double el = eta/lam;

    double rl = hypot(1.0, el);

    double sl = el  + lam*x_inv;

    double fc_lm1;

    fc_lm1 = (fcl*sl + fpl)/rl;

    fpl    =  fc_lm1*sl - fcl*rl;

    fcl    =  fc_lm1;

    lam -= 1.0;

  }

  *F_lam_min  = fcl;

  *Fp_lam_min = fpl;  

  return GSL_SUCCESS;

}

/* Evolve the forward recurrence for G,G'.

 *

 *   G_{lam+1}  = (S_lam G_lam - G_lam')/R_lam

 *   G_{lam+1}' = R_{lam+1} G_lam - S_lam G_{lam+1}

 *

 * where S_lam and R_lam are as above in the F recursion.

 */

static

int

coulomb_G_recur(const double lam_min, const int kmax,

                const double eta, const double x,

                const double G_lam_min, const double Gp_lam_min,

                double * G_lam_max, double * Gp_lam_max

                )

{

  double x_inv = 1.0/x;

  double gcl = G_lam_min;

  double gpl = Gp_lam_min;

  double lam = lam_min + 1.0;

  int k;

  for(k=1; k<=kmax; k++) {

    double el = eta/lam;

    double rl = hypot(1.0, el);

    double sl = el + lam*x_inv;

    double gcl1 = (sl*gcl - gpl)/rl;

    gpl   = rl*gcl - sl*gcl1;

    gcl   = gcl1;

    lam += 1.0;

  }

  *G_lam_max  = gcl;

  *Gp_lam_max = gpl;

  return GSL_SUCCESS;

}

/* Evaluate the first continued fraction, giving

 * the ratio F'/F at the upper lambda value.

 * We also determine the sign of F at that point,

 * since it is the sign of the last denominator

 * in the continued fraction.

 */

static

int

coulomb_CF1(double lambda,

            double eta, double x,

            double * fcl_sign,

            double * result,

            int * count

            )

{

  const double CF1_small = 1.e-30;

  const double CF1_abort = 1.0e+05;

  const double CF1_acc   = 2.0*GSL_DBL_EPSILON;

  const double x_inv     = 1.0/x;

  const double px        = lambda + 1.0 + CF1_abort;

  double pk = lambda + 1.0;

  double F  = eta/pk + pk*x_inv;

  double D, C;

  double df;

  *fcl_sign = 1.0;

  *count = 0;

  if(fabs(F) < CF1_small) F = CF1_small;

  D = 0.0;

  C = F;

  do {

    double pk1 = pk + 1.0;

    double ek  = eta / pk;

    double rk2 = 1.0 + ek*ek;

    double tk  = (pk + pk1)*(x_inv + ek/pk1);

    D   =  tk - rk2 * D;

    C   =  tk - rk2 / C;

    if(fabs(C) < CF1_small) C = CF1_small;

    if(fabs(D) < CF1_small) D = CF1_small;

    D = 1.0/D;

    df = D * C;

    F  = F * df;

    if(D < 0.0) {

      /* sign of result depends on sign of denominator */

      *fcl_sign = - *fcl_sign;

    }

    pk = pk1;

    if( pk > px ) {

      *result = F;

      GSL_ERROR ("error", GSL_ERUNAWAY);

    }

    ++(*count);

  }

  while(fabs(df-1.0) > CF1_acc);

  *result = F;

  return GSL_SUCCESS;

}

#if 0

static

int

old_coulomb_CF1(const double lambda,

                double eta, double x,

                double * fcl_sign,

                double * result

                )

{

  const double CF1_abort = 1.e5;

  const double CF1_acc   = 10.0*GSL_DBL_EPSILON;

  const double x_inv     = 1.0/x;

  const double px        = lambda + 1.0 + CF1_abort;

  double pk = lambda + 1.0;

  double D;

  double df;

  double F;

  double p;

  double pk1;

  double ek;

  double fcl = 1.0;

  double tk;

  while(1) {

    ek = eta/pk;

    F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;

    pk1 = pk + 1.0;

    if(fabs(eta*x + pk*pk1) > CF1_acc) break;

    fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));

    pk = 2.0 + pk;

  }

  D  = 1.0/((pk + pk1)*(x_inv + ek/pk1));

  df = -fcl*(1.0 + ek*ek)*D;

  if(fcl != 1.0) fcl = -1.0;

  if(D    < 0.0) fcl = -fcl;

  F = F + df;

  p = 1.0;

  do {

    pk = pk1;

    pk1 = pk + 1.0;

    ek  = eta / pk;

    tk  = (pk + pk1)*(x_inv + ek/pk1);

    D   =  tk - D*(1.0+ek*ek);

    if(fabs(D) < sqrt(CF1_acc)) {

      p += 1.0;

      if(p > 2.0) {

        printf("HELP............\n");

      }

    }

    D = 1.0/D;

    if(D < 0.0) {

      /* sign of result depends on sign of denominator */

      fcl = -fcl;

    }

    df = df*(D*tk - 1.0);

    F  = F + df;

    if( pk > px ) {

      GSL_ERROR ("error", GSL_ERUNAWAY);

    }

  }

  while(fabs(df) > fabs(F)*CF1_acc);

  *fcl_sign = fcl;

  *result = F;

  return GSL_SUCCESS;

}

#endif /* 0 */

/* Evaluate the second continued fraction to 

 * obtain the ratio

 *    (G' + i F')/(G + i F) := P + i Q

 * at the specified lambda value.

 */

static

int

coulomb_CF2(const double lambda, const double eta, const double x,

            double * result_P, double * result_Q, int * count

            )

{

  int status = GSL_SUCCESS;

  const double CF2_acc   = 4.0*GSL_DBL_EPSILON;

  const double CF2_abort = 2.0e+05;

  const double wi    = 2.0*eta;

  const double x_inv = 1.0/x;

  const double e2mm1 = eta*eta + lambda*(lambda + 1.0);

  double ar = -e2mm1;

  double ai =  eta;

  double br =  2.0*(x - eta);

  double bi =  2.0;

  double dr =  br/(br*br + bi*bi);

  double di = -bi/(br*br + bi*bi);

  double dp = -x_inv*(ar*di + ai*dr);

  double dq =  x_inv*(ar*dr - ai*di);

  double A, B, C, D;

  double pk =  0.0;

  double P  =  0.0;

  double Q  =  1.0 - eta*x_inv;

  *count = 0;

  do {

    P += dp;

    Q += dq;

    pk += 2.0;

    ar += pk;

    ai += wi;

    bi += 2.0;

    D  = ar*dr - ai*di + br;

    di = ai*dr + ar*di + bi;

    C  = 1.0/(D*D + di*di);

    dr =  C*D;

    di = -C*di;

    A  = br*dr - bi*di - 1.;

    B  = bi*dr + br*di;

    C  = dp*A  - dq*B;

    dq = dp*B  + dq*A;

    dp = C;

    if(pk > CF2_abort) {

      status = GSL_ERUNAWAY;

      break;

    }

    ++(*count);

  }

  while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);

  if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {

    status = GSL_ELOSS;

  }

  *result_P = P;

  *result_Q = Q;

  return status;

}

/* WKB evaluation of F, G. Assumes  0 < x < turning point.

 * Overflows are trapped, GSL_EOVRFLW is signalled,

 * and an exponent is returned such that:

 *

 *   result_F = fjwkb * exp(-exponent)

 *   result_G = gjwkb * exp( exponent)

 *

 * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]

 *

 * Unfortunately, this is not very accurate in general. The

 * test cases typically have 3-4 digits of precision. One could

 * argue that this is ok for general use because, for instance,

 * F is exponentially small in this region and so the absolute

 * accuracy is still roughly acceptable. But it would be better

 * to have a systematic method for improving the precision. See

 * the Abad+Sesma method discussion below.

 */

static

int

coulomb_jwkb(const double lam, const double eta, const double x,

             gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,

             double * exponent)

{

  const double llp1      = lam*(lam+1.0) + 6.0/35.0;

  const double llp1_eff  = GSL_MAX(llp1, 0.0);

  const double rho_ghalf = sqrt(x*(2.0*eta - x) + llp1_eff);

  const double sinh_arg  = sqrt(llp1_eff/(eta*eta+llp1_eff)) * rho_ghalf / x;

  const double sinh_inv  = log(sinh_arg + hypot(1.0,sinh_arg));

  const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);

  const double zeta_half = pow(3.0*phi/2.0, 1.0/3.0);

  const double prefactor = sqrt(M_PI*phi*x/(6.0 * rho_ghalf));

  double F = prefactor * 3.0/zeta_half;

  double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */

  double F_exp;

  double G_exp;

  const double airy_scale_exp = phi;

  gsl_sf_result ai;

  gsl_sf_result bi;

  gsl_sf_airy_Ai_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &ai;;

  gsl_sf_airy_Bi_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &bi);

  F *= ai.val;

  G *= bi.val;

  F_exp = log(F) - airy_scale_exp;

  G_exp = log(G) + airy_scale_exp;

  if(G_exp >= GSL_LOG_DBL_MAX) {

    fjwkb->val = F;

    gjwkb->val = G;

    fjwkb->err = 1.0e-3 * fabs(F); /* FIXME: real error here ... could be smaller */

    gjwkb->err = 1.0e-3 * fabs(G);

    *exponent = airy_scale_exp;

    GSL_ERROR ("error", GSL_EOVRFLW);

  }

  else {

    fjwkb->val = exp(F_exp);

    gjwkb->val = exp(G_exp);

    fjwkb->err = 1.0e-3 * fabs(fjwkb->val);

    gjwkb->err = 1.0e-3 * fabs(gjwkb->val);

    *exponent = 0.0;

    return GSL_SUCCESS;

  }

}

/* Asymptotic evaluation of F and G below the minimal turning point.

 *

 * This is meant to be a drop-in replacement for coulomb_jwkb().

 * It uses the expressions in [Abad+Sesma]. This requires some

 * work because I am not sure where it is valid. They mumble

 * something about |x| < |lam|^(-1/2) or 8|eta x| > lam when |x| < 1.

 * This seems true, but I thought the result was based on a uniform

 * expansion and could be controlled by simply using more terms.

 */

#if 0