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

   1 /* specfunc/coulomb.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 /* Evaluation of Coulomb wave functions F_L(eta, x), G_L(eta, x),
  23  * and their derivatives. A combination of Steed's method, asymptotic
  24  * results, and power series.
  25  *
  26  * Steed's method:
  27  *  [Barnett, CPC 21, 297 (1981)]
  28  * Power series and other methods:
  29  *  [Biedenharn et al., PR 97, 542 (1954)]
  30  *  [Bardin et al., CPC 3, 73 (1972)]
  31  *  [Abad+Sesma, CPC 71, 110 (1992)]
  32  */
  33 #include <config.h>
  34 #include <gsl/gsl_math.h>
  35 #include <gsl/gsl_errno.h>
  36 #include <gsl/gsl_sf_exp.h>
  37 #include <gsl/gsl_sf_psi.h>
  38 #include <gsl/gsl_sf_airy.h>
  39 #include <gsl/gsl_sf_pow_int.h>
  40 #include <gsl/gsl_sf_gamma.h>
  41 #include <gsl/gsl_sf_coulomb.h>
  42
  43 #include "error.h"
  44
  45 /* the L=0 normalization constant
  46  * [Abramowitz+Stegun 14.1.8]
  47  */
  48 static
  49 double
  50 C0sq(double eta)
  51 {
  52   double twopieta = 2.0*M_PI*eta;
  53
  54   if(fabs(eta) < GSL_DBL_EPSILON) {
  55     return 1.0;
  56   }
  57   else if(twopieta > GSL_LOG_DBL_MAX) {
  58     return 0.0;
  59   }
  60   else {
  61     gsl_sf_result scale;
  62     gsl_sf_expm1_e(twopieta, &scale);
  63     return twopieta/scale.val;
  64   }
  65 }
  66
  67
  68 /* the full definition of C_L(eta) for any valid L and eta
  69  * [Abramowitz and Stegun 14.1.7]
  70  * This depends on the complex gamma function. For large
  71  * arguments the phase of the complex gamma function is not
  72  * very accurately determined. However the modulus is, and that
  73  * is all that we need to calculate C_L.
  74  *
  75  * This is not valid for L <= -3/2  or  L = -1.
  76  */
  77 static
  78 int
  79 CLeta(double L, double eta, gsl_sf_result * result)
  80 {
  81   gsl_sf_result ln1; /* log of numerator Gamma function */
  82   gsl_sf_result ln2; /* log of denominator Gamma function */
  83   double sgn = 1.0;
  84   double arg_val, arg_err;
  85
  86   if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
  87     gsl_sf_lngamma_e(L+1.0, &ln1);
  88   }
  89   else {
  90     gsl_sf_result p1;                 /* phase of numerator Gamma -- not used */
  91     gsl_sf_lngamma_complex_e(L+1.0, eta, &ln1, &p1); /* should be ok */
  92   }
  93
  94   gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
  95   if(L < -1.0) sgn = -sgn;
  96
  97   arg_val  = L*M_LN2 - 0.5*eta*M_PI + ln1.val - ln2.val;
  98   arg_err  = ln1.err + ln2.err;
  99   arg_err += GSL_DBL_EPSILON * (fabs(L*M_LN2) + fabs(0.5*eta*M_PI));
 100   return gsl_sf_exp_err_e(arg_val, arg_err, result);
 101 }
 102
 103
 104 int
 105 gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
 106 {
 107   /* CHECK_POINTER(result) */
 108
 109   if(lam <= -1.0) {
 110     DOMAIN_ERROR(result);
 111   }
 112   else if(fabs(lam) < GSL_DBL_EPSILON) {
 113     /* saves a calculation of complex_lngamma(), otherwise not necessary */
 114     result->val = sqrt(C0sq(eta));
 115     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
 116     return GSL_SUCCESS;
 117   }
 118   else {
 119     return CLeta(lam, eta, result);
 120   }
 121 }
 122
 123
 124 /* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
 125  */
 126 int
 127 gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
 128 {
 129   int k;
 130   gsl_sf_result cl_0;
 131   gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
 132   cl[0] = cl_0.val;
 133
 134   for(k=1; k<=kmax; k++) {
 135     double L = lam_min + k;
 136     cl[k] = cl[k-1] * hypot(L, eta)/(L*(2.0*L+1.0));
 137   }
 138
 139   return GSL_SUCCESS;
 140 }
 141
 142
 143 /* Evaluate the series for Phi_L(eta,x) and Phi_L*(eta,x)
 144  * [Abramowitz+Stegun 14.1.5]
 145  * [Abramowitz+Stegun 14.1.13]
 146  *
 147  * The sequence of coefficients A_k^L is
 148  * manifestly well-controlled for L >= -1/2
 149  * and eta < 10.
 150  *
 151  * This makes sense since this is the region
 152  * away from threshold, and you expect
 153  * the evaluation to become easier as you
 154  * get farther from threshold.
 155  *
 156  * Empirically, this is quite well-behaved for
 157  *   L >= -1/2
 158  *   eta < 10
 159  *   x   < 10
 160  */
 161 #if 0
 162 static
 163 int
 164 coulomb_Phi_series(const double lam, const double eta, const double x,
 165                    double * result, double * result_star)
 166 {
 167   int kmin =   5;
 168   int kmax = 200;
 169   int k;
 170   double Akm2 = 1.0;
 171   double Akm1 = eta/(lam+1.0);
 172   double Ak;
 173
 174   double xpow = x;
 175   double sum  = Akm2 + Akm1*x;
 176   double sump = (lam+1.0)*Akm2 + (lam+2.0)*Akm1*x;
 177   double prev_abs_del   = fabs(Akm1*x);
 178   double prev_abs_del_p = (lam+2.0) * prev_abs_del;
 179
 180   for(k=2; k<kmax; k++) {
 181     double del;
 182     double del_p;
 183     double abs_del;
 184     double abs_del_p;
 185
 186     Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));
 187
 188     xpow *= x;
 189     del   = Ak*xpow;
 190     del_p = (k+lam+1.0)*del;
 191     sum  += del;
 192     sump += del_p;
 193
 194     abs_del   = fabs(del);
 195     abs_del_p = fabs(del_p);
 196
 197     if(          abs_del/(fabs(sum)+abs_del)          < GSL_DBL_EPSILON
 198        &&   prev_abs_del/(fabs(sum)+prev_abs_del)     < GSL_DBL_EPSILON
 199        &&      abs_del_p/(fabs(sump)+abs_del_p)       < GSL_DBL_EPSILON
 200        && prev_abs_del_p/(fabs(sump)+prev_abs_del_p)  < GSL_DBL_EPSILON
 201        && k > kmin
 202        ) break;
 203
 204     /* We need to keep track of the previous delta because when
 205      * eta is near zero the odd terms of the sum are very small
 206      * and this could lead to premature termination.
 207      */
 208     prev_abs_del   = abs_del;
 209     prev_abs_del_p = abs_del_p;
 210
 211     Akm2 = Akm1;
 212     Akm1 = Ak;
 213   }
 214
 215   *result      = sum;
 216   *result_star = sump;
 217
 218   if(k==kmax) {
 219     GSL_ERROR ("error", GSL_EMAXITER);
 220   }
 221   else {
 222     return GSL_SUCCESS;
 223   }
 224 }
 225 #endif /* 0 */
 226
 227
 228 /* Determine the connection phase, phi_lambda.
 229  * See coulomb_FG_series() below. We have
 230  * to be careful about sin(phi)->0. Note that
 231  * there is an underflow condition for large
 232  * positive eta in any case.
 233  */
 234 static
 235 int
 236 coulomb_connection(const double lam, const double eta,
 237                    double * cos_phi, double * sin_phi)
 238 {
 239   if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
 240     *cos_phi = 1.0;
 241     *sin_phi = 0.0;
 242     GSL_ERROR ("error", GSL_EUNDRFLW);
 243   }
 244   else if(eta > -GSL_LOG_DBL_EPSILON/(4.0*M_PI)) {
 245     const double eps = 2.0 * exp(-2.0*M_PI*eta);
 246     const double tpl = tan(M_PI * lam);
 247     const double dth = eps * tpl / (tpl*tpl + 1.0);
 248     *cos_phi = -1.0 + 0.5 * dth*dth;
 249     *sin_phi = -dth;
 250     return GSL_SUCCESS;
 251   }
 252   else {
 253     double X   = tanh(M_PI * eta) / tan(M_PI * lam);
 254     double phi = -atan(X) - (lam + 0.5) * M_PI;
 255     *cos_phi = cos(phi);
 256     *sin_phi = sin(phi);
 257     return GSL_SUCCESS;
 258   }
 259 }
 260
 261
 262 /* Evaluate the Frobenius series for F_lam(eta,x) and G_lam(eta,x).
 263  * Homegrown algebra. Evaluates the series for F_{lam} and
 264  * F_{-lam-1}, then uses
 265  *    G_{lam} = (F_{lam} cos(phi) - F_{-lam-1}) / sin(phi)
 266  * where
 267  *    phi = Arg[Gamma[1+lam+I eta]] - Arg[Gamma[-lam + I eta]] - (lam+1/2)Pi
 268  *        = Arg[Sin[Pi(-lam+I eta)] - (lam+1/2)Pi
 269  *        = atan2(-cos(lam Pi)sinh(eta Pi), -sin(lam Pi)cosh(eta Pi)) - (lam+1/2)Pi
 270  *
 271  *        = -atan(X) - (lam+1/2) Pi,  X = tanh(eta Pi)/tan(lam Pi)
 272  *
 273  * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
 274  */
 275 static
 276 int
 277 coulomb_FG_series(const double lam, const double eta, const double x,
 278                   gsl_sf_result * F, gsl_sf_result * G)
 279 {
 280   const int max_iter = 800;
 281   gsl_sf_result ClamA;
 282   gsl_sf_result ClamB;
 283   int stat_A = CLeta(lam, eta, &ClamA);
 284   int stat_B = CLeta(-lam-1.0, eta, &ClamB);
 285   const double tlp1 = 2.0*lam + 1.0;
 286   const double pow_x = pow(x, lam);
 287   double cos_phi_lam;
 288   double sin_phi_lam;
 289
 290   double uA_mm2 = 1.0;                  /* uA sum is for F_{lam} */
 291   double uA_mm1 = x*eta/(lam+1.0);
 292   double uA_m;
 293   double uB_mm2 = 1.0;                  /* uB sum is for F_{-lam-1} */
 294   double uB_mm1 = -x*eta/lam;
 295   double uB_m;
 296   double A_sum = uA_mm2 + uA_mm1;
 297   double B_sum = uB_mm2 + uB_mm1;
 298   double A_abs_del_prev = fabs(A_sum);
 299   double B_abs_del_prev = fabs(B_sum);
 300   gsl_sf_result FA, FB;
 301   int m = 2;
 302
 303   int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);
 304
 305   if(stat_conn == GSL_EUNDRFLW) {
 306     F->val = 0.0;  /* FIXME: should this be set to Inf too like G? */
 307     F->err = 0.0;
 308     OVERFLOW_ERROR(G);
 309   }
 310
 311   while(m < max_iter) {
 312     double abs_dA;
 313     double abs_dB;
 314     uA_m = x*(2.0*eta*uA_mm1 - x*uA_mm2)/(m*(m+tlp1));
 315     uB_m = x*(2.0*eta*uB_mm1 - x*uB_mm2)/(m*(m-tlp1));
 316     A_sum += uA_m;
 317     B_sum += uB_m;
 318     abs_dA = fabs(uA_m);
 319     abs_dB = fabs(uB_m);
 320     if(m > 15) {
 321       /* Don't bother checking until we have gone out a little ways;
 322        * a minor optimization. Also make sure to check both the
 323        * current and the previous increment because the odd and even
 324        * terms of the sum can have very different behaviour, depending
 325        * on the value of eta.
 326        */
 327       double max_abs_dA = GSL_MAX(abs_dA, A_abs_del_prev);
 328       double max_abs_dB = GSL_MAX(abs_dB, B_abs_del_prev);
 329       double abs_A = fabs(A_sum);
 330       double abs_B = fabs(B_sum);
 331       if(   max_abs_dA/(max_abs_dA + abs_A) < 4.0*GSL_DBL_EPSILON
 332          && max_abs_dB/(max_abs_dB + abs_B) < 4.0*GSL_DBL_EPSILON
 333          ) break;
 334     }
 335     A_abs_del_prev = abs_dA;
 336     B_abs_del_prev = abs_dB;
 337     uA_mm2 = uA_mm1;
 338     uA_mm1 = uA_m;
 339     uB_mm2 = uB_mm1;
 340     uB_mm1 = uB_m;
 341     m++;
 342   }
 343
 344   FA.val = A_sum * ClamA.val * pow_x * x;
 345   FA.err = fabs(A_sum) * ClamA.err * pow_x * x + 2.0*GSL_DBL_EPSILON*fabs(FA.val);
 346   FB.val = B_sum * ClamB.val / pow_x;
 347   FB.err = fabs(B_sum) * ClamB.err / pow_x + 2.0*GSL_DBL_EPSILON*fabs(FB.val);
 348
 349   F->val = FA.val;
 350   F->err = FA.err;
 351
 352   G->val = (FA.val * cos_phi_lam - FB.val)/sin_phi_lam;
 353   G->err = (FA.err * fabs(cos_phi_lam) + FB.err)/fabs(sin_phi_lam);
 354
 355   if(m >= max_iter)
 356     GSL_ERROR ("error", GSL_EMAXITER);
 357   else
 358     return GSL_ERROR_SELECT_2(stat_A, stat_B);
 359 }
 360
 361
 362 /* Evaluate the Frobenius series for F_0(eta,x) and G_0(eta,x).
 363  * See [Bardin et al., CPC 3, 73 (1972), (14)-(17)];
 364  * note the misprint in (17): nu_0=1 is correct, not nu_0=0.
 365  */
 366 static
 367 int
 368 coulomb_FG0_series(const double eta, const double x,
 369                    gsl_sf_result * F, gsl_sf_result * G)
 370 {
 371   const int max_iter = 800;
 372   const double x2  = x*x;
 373   const double tex = 2.0*eta*x;
 374   gsl_sf_result C0;
 375   int stat_CL = CLeta(0.0, eta, &C0);
 376   gsl_sf_result r1pie;
 377   int psi_stat = gsl_sf_psi_1piy_e(eta, &r1pie);
 378   double u_mm2 = 0.0;  /* u_0 */
 379   double u_mm1 = x;    /* u_1 */
 380   double u_m;
 381   double v_mm2 = 1.0;                               /* nu_0 */
 382   double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val);   /* nu_1 */
 383   double v_m;
 384   double u_sum = u_mm2 + u_mm1;
 385   double v_sum = v_mm2 + v_mm1;
 386   double u_abs_del_prev = fabs(u_sum);
 387   double v_abs_del_prev = fabs(v_sum);
 388   int m = 2;
 389   double u_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(u_sum);
 390   double v_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(v_sum);
 391   double ln2x = log(2.0*x);
 392
 393   while(m < max_iter) {
 394     double abs_du;
 395     double abs_dv;
 396     double m_mm1 = m*(m-1.0);
 397     u_m = (tex*u_mm1 - x2*u_mm2)/m_mm1;
 398     v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*eta*(2*m-1)*u_m)/m_mm1;
 399     u_sum += u_m;
 400     v_sum += v_m;
 401     abs_du = fabs(u_m);
 402     abs_dv = fabs(v_m);
 403     u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
 404     v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
 405     if(m > 15) {
 406       /* Don't bother checking until we have gone out a little ways;
 407        * a minor optimization. Also make sure to check both the
 408        * current and the previous increment because the odd and even
 409        * terms of the sum can have very different behaviour, depending
 410        * on the value of eta.
 411        */
 412       double max_abs_du = GSL_MAX(abs_du, u_abs_del_prev);
 413       double max_abs_dv = GSL_MAX(abs_dv, v_abs_del_prev);
 414       double abs_u = fabs(u_sum);
 415       double abs_v = fabs(v_sum);
 416       if(   max_abs_du/(max_abs_du + abs_u) < 40.0*GSL_DBL_EPSILON
 417          && max_abs_dv/(max_abs_dv + abs_v) < 40.0*GSL_DBL_EPSILON
 418          ) break;
 419     }
 420     u_abs_del_prev = abs_du;
 421     v_abs_del_prev = abs_dv;
 422     u_mm2 = u_mm1;
 423     u_mm1 = u_m;
 424     v_mm2 = v_mm1;
 425     v_mm1 = v_m;
 426     m++;
 427   }
 428
 429   F->val  = C0.val * u_sum;
 430   F->err  = C0.err * fabs(u_sum);
 431   F->err += fabs(C0.val) * u_sum_err;
 432   F->err += 2.0 * GSL_DBL_EPSILON * fabs(F->val);
 433
 434   G->val  = (v_sum + 2.0*eta*u_sum * ln2x) / C0.val;
 435   G->err  = (fabs(v_sum) + fabs(2.0*eta*u_sum * ln2x)) / fabs(C0.val) * fabs(C0.err/C0.val);
 436   G->err += (v_sum_err + fabs(2.0*eta*u_sum_err*ln2x)) / fabs(C0.val);
 437   G->err += 2.0 * GSL_DBL_EPSILON * fabs(G->val);
 438
 439   if(m == max_iter)
 440     GSL_ERROR ("error", GSL_EMAXITER);
 441   else
 442     return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
 443 }
 444
 445
 446 /* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
 447  * Homegrown algebra.
 448  */
 449 static
 450 int
 451 coulomb_FGmhalf_series(const double eta, const double x,
 452                        gsl_sf_result * F, gsl_sf_result * G)
 453 {
 454   const int max_iter = 800;
 455   const double rx  = sqrt(x);
 456   const double x2  = x*x;
 457   const double tex = 2.0*eta*x;
 458   gsl_sf_result Cmhalf;
 459   int stat_CL = CLeta(-0.5, eta, &Cmhalf);
 460   double u_mm2 = 1.0;                      /* u_0 */
 461   double u_mm1 = tex * u_mm2;              /* u_1 */
 462   double u_m;
 463   double v_mm2, v_mm1, v_m;
 464   double f_sum, g_sum;
 465   double tmp1;
 466   gsl_sf_result rpsi_1pe;
 467   gsl_sf_result rpsi_1p2e;
 468   int m = 2;
 469
 470   gsl_sf_psi_1piy_e(eta,     &rpsi_1pe);
 471   gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);
 472
 473   v_mm2 = 2.0*M_EULER - M_LN2 - rpsi_1pe.val + 2.0*rpsi_1p2e.val;
 474   v_mm1 = tex*(v_mm2 - 2.0*u_mm2);
 475
 476   f_sum = u_mm2 + u_mm1;
 477   g_sum = v_mm2 + v_mm1;
 478
 479   while(m < max_iter) {
 480     double m2 = m*m;
 481     u_m = (tex*u_mm1 - x2*u_mm2)/m2;
 482     v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*m*u_m)/m2;
 483     f_sum += u_m;
 484     g_sum += v_m;
 485     if(   f_sum != 0.0
 486        && g_sum != 0.0
 487        && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;
 488     u_mm2 = u_mm1;
 489     u_mm1 = u_m;
 490     v_mm2 = v_mm1;
 491     v_mm1 = v_m;
 492     m++;
 493   }
 494
 495   F->val = Cmhalf.val * rx * f_sum;
 496   F->err = Cmhalf.err * fabs(rx * f_sum) + 2.0*GSL_DBL_EPSILON*fabs(F->val);
 497
 498   tmp1 = f_sum*log(x);
 499   G->val = -rx*(tmp1 + g_sum)/Cmhalf.val;
 500   G->err = fabs(rx)*(fabs(tmp1) + fabs(g_sum))/fabs(Cmhalf.val) * fabs(Cmhalf.err/Cmhalf.val);
 501
 502   if(m == max_iter)
 503     GSL_ERROR ("error", GSL_EMAXITER);
 504   else
 505     return stat_CL;
 506 }
 507
 508
 509 /* Evolve the backwards recurrence for F,F'.
 510  *
 511  *    F_{lam-1}  = (S_lam F_lam + F_lam') / R_lam
 512  *    F_{lam-1}' = (S_lam F_{lam-1} - R_lam F_lam)
 513  * where
 514  *    R_lam = sqrt(1 + (eta/lam)^2)
 515  *    S_lam = lam/x + eta/lam
 516  *
 517  */
 518 static
 519 int
 520 coulomb_F_recur(double lam_min, int kmax,
 521                 double eta, double x,
 522                 double F_lam_max, double Fp_lam_max,
 523                 double * F_lam_min, double * Fp_lam_min
 524                 )
 525 {
 526   double x_inv = 1.0/x;
 527   double fcl = F_lam_max;
 528   double fpl = Fp_lam_max;
 529   double lam_max = lam_min + kmax;
 530   double lam = lam_max;
 531   int k;
 532
 533   for(k=kmax-1; k>=0; k--) {
 534     double el = eta/lam;
 535     double rl = hypot(1.0, el);
 536     double sl = el  + lam*x_inv;
 537     double fc_lm1;
 538     fc_lm1 = (fcl*sl + fpl)/rl;
 539     fpl    =  fc_lm1*sl - fcl*rl;
 540     fcl    =  fc_lm1;
 541     lam -= 1.0;
 542   }
 543
 544   *F_lam_min  = fcl;
 545   *Fp_lam_min = fpl;
 546   return GSL_SUCCESS;
 547 }
 548
 549
 550 /* Evolve the forward recurrence for G,G'.
 551  *
 552  *   G_{lam+1}  = (S_lam G_lam - G_lam')/R_lam
 553  *   G_{lam+1}' = R_{lam+1} G_lam - S_lam G_{lam+1}
 554  *
 555  * where S_lam and R_lam are as above in the F recursion.
 556  */
 557 static
 558 int
 559 coulomb_G_recur(const double lam_min, const int kmax,
 560                 const double eta, const double x,
 561                 const double G_lam_min, const double Gp_lam_min,
 562                 double * G_lam_max, double * Gp_lam_max
 563                 )
 564 {
 565   double x_inv = 1.0/x;
 566   double gcl = G_lam_min;
 567   double gpl = Gp_lam_min;
 568   double lam = lam_min + 1.0;
 569   int k;
 570
 571   for(k=1; k<=kmax; k++) {
 572     double el = eta/lam;
 573     double rl = hypot(1.0, el);
 574     double sl = el + lam*x_inv;
 575     double gcl1 = (sl*gcl - gpl)/rl;
 576     gpl   = rl*gcl - sl*gcl1;
 577     gcl   = gcl1;
 578     lam += 1.0;
 579   }
 580
 581   *G_lam_max  = gcl;
 582   *Gp_lam_max = gpl;
 583   return GSL_SUCCESS;
 584 }
 585
 586
 587 /* Evaluate the first continued fraction, giving
 588  * the ratio F'/F at the upper lambda value.
 589  * We also determine the sign of F at that point,
 590  * since it is the sign of the last denominator
 591  * in the continued fraction.
 592  */
 593 static
 594 int
 595 coulomb_CF1(double lambda,
 596             double eta, double x,
 597             double * fcl_sign,
 598             double * result,
 599             int * count
 600             )
 601 {
 602   const double CF1_small = 1.e-30;
 603   const double CF1_abort = 1.0e+05;
 604   const double CF1_acc   = 2.0*GSL_DBL_EPSILON;
 605   const double x_inv     = 1.0/x;
 606   const double px        = lambda + 1.0 + CF1_abort;
 607
 608   double pk = lambda + 1.0;
 609   double F  = eta/pk + pk*x_inv;
 610   double D, C;
 611   double df;
 612
 613   *fcl_sign = 1.0;
 614   *count = 0;
 615
 616   if(fabs(F) < CF1_small) F = CF1_small;
 617   D = 0.0;
 618   C = F;
 619
 620   do {
 621     double pk1 = pk + 1.0;
 622     double ek  = eta / pk;
 623     double rk2 = 1.0 + ek*ek;
 624     double tk  = (pk + pk1)*(x_inv + ek/pk1);
 625     D   =  tk - rk2 * D;
 626     C   =  tk - rk2 / C;
 627     if(fabs(C) < CF1_small) C = CF1_small;
 628     if(fabs(D) < CF1_small) D = CF1_small;
 629     D = 1.0/D;
 630     df = D * C;
 631     F  = F * df;
 632     if(D < 0.0) {
 633       /* sign of result depends on sign of denominator */
 634       *fcl_sign = - *fcl_sign;
 635     }
 636     pk = pk1;
 637     if( pk > px ) {
 638       *result = F;
 639       GSL_ERROR ("error", GSL_ERUNAWAY);
 640     }
 641     ++(*count);
 642   }
 643   while(fabs(df-1.0) > CF1_acc);
 644
 645   *result = F;
 646   return GSL_SUCCESS;
 647 }
 648
 649
 650 #if 0
 651 static
 652 int
 653 old_coulomb_CF1(const double lambda,
 654                 double eta, double x,
 655                 double * fcl_sign,
 656                 double * result
 657                 )
 658 {
 659   const double CF1_abort = 1.e5;
 660   const double CF1_acc   = 10.0*GSL_DBL_EPSILON;
 661   const double x_inv     = 1.0/x;
 662   const double px        = lambda + 1.0 + CF1_abort;
 663
 664   double pk = lambda + 1.0;
 665
 666   double D;
 667   double df;
 668
 669   double F;
 670   double p;
 671   double pk1;
 672   double ek;
 673
 674   double fcl = 1.0;
 675
 676   double tk;
 677
 678   while(1) {
 679     ek = eta/pk;
 680     F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;
 681     pk1 = pk + 1.0;
 682     if(fabs(eta*x + pk*pk1) > CF1_acc) break;
 683     fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));
 684     pk = 2.0 + pk;
 685   }
 686
 687   D  = 1.0/((pk + pk1)*(x_inv + ek/pk1));
 688   df = -fcl*(1.0 + ek*ek)*D;
 689
 690   if(fcl != 1.0) fcl = -1.0;
 691   if(D    < 0.0) fcl = -fcl;
 692
 693   F = F + df;
 694
 695   p = 1.0;
 696   do {
 697     pk = pk1;
 698     pk1 = pk + 1.0;
 699     ek  = eta / pk;
 700     tk  = (pk + pk1)*(x_inv + ek/pk1);
 701     D   =  tk - D*(1.0+ek*ek);
 702     if(fabs(D) < sqrt(CF1_acc)) {
 703       p += 1.0;
 704       if(p > 2.0) {
 705         printf("HELP............\n");
 706       }
 707     }
 708     D = 1.0/D;
 709     if(D < 0.0) {
 710       /* sign of result depends on sign of denominator */
 711       fcl = -fcl;
 712     }
 713     df = df*(D*tk - 1.0);
 714     F  = F + df;
 715     if( pk > px ) {
 716       GSL_ERROR ("error", GSL_ERUNAWAY);
 717     }
 718   }
 719   while(fabs(df) > fabs(F)*CF1_acc);
 720
 721   *fcl_sign = fcl;
 722   *result = F;
 723   return GSL_SUCCESS;
 724 }
 725 #endif /* 0 */
 726
 727
 728 /* Evaluate the second continued fraction to
 729  * obtain the ratio
 730  *    (G' + i F')/(G + i F) := P + i Q
 731  * at the specified lambda value.
 732  */
 733 static
 734 int
 735 coulomb_CF2(const double lambda, const double eta, const double x,
 736             double * result_P, double * result_Q, int * count
 737             )
 738 {
 739   int status = GSL_SUCCESS;
 740
 741   const double CF2_acc   = 4.0*GSL_DBL_EPSILON;
 742   const double CF2_abort = 2.0e+05;
 743
 744   const double wi    = 2.0*eta;
 745   const double x_inv = 1.0/x;
 746   const double e2mm1 = eta*eta + lambda*(lambda + 1.0);
 747
 748   double ar = -e2mm1;
 749   double ai =  eta;
 750
 751   double br =  2.0*(x - eta);
 752   double bi =  2.0;
 753
 754   double dr =  br/(br*br + bi*bi);
 755   double di = -bi/(br*br + bi*bi);
 756
 757   double dp = -x_inv*(ar*di + ai*dr);
 758   double dq =  x_inv*(ar*dr - ai*di);
 759
 760   double A, B, C, D;
 761
 762   double pk =  0.0;
 763   double P  =  0.0;
 764   double Q  =  1.0 - eta*x_inv;
 765
 766   *count = 0;
 767
 768   do {
 769     P += dp;
 770     Q += dq;
 771     pk += 2.0;
 772     ar += pk;
 773     ai += wi;
 774     bi += 2.0;
 775     D  = ar*dr - ai*di + br;
 776     di = ai*dr + ar*di + bi;
 777     C  = 1.0/(D*D + di*di);
 778     dr =  C*D;
 779     di = -C*di;
 780     A  = br*dr - bi*di - 1.;
 781     B  = bi*dr + br*di;
 782     C  = dp*A  - dq*B;
 783     dq = dp*B  + dq*A;
 784     dp = C;
 785     if(pk > CF2_abort) {
 786       status = GSL_ERUNAWAY;
 787       break;
 788     }
 789     ++(*count);
 790   }
 791   while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);
 792
 793   if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {
 794     status = GSL_ELOSS;
 795   }
 796
 797   *result_P = P;
 798   *result_Q = Q;
 799   return status;
 800 }
 801
 802
 803 /* WKB evaluation of F, G. Assumes  0 < x < turning point.
 804  * Overflows are trapped, GSL_EOVRFLW is signalled,
 805  * and an exponent is returned such that:
 806  *
 807  *   result_F = fjwkb * exp(-exponent)
 808  *   result_G = gjwkb * exp( exponent)
 809  *
 810  * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]
 811  *
 812  * Unfortunately, this is not very accurate in general. The
 813  * test cases typically have 3-4 digits of precision. One could
 814  * argue that this is ok for general use because, for instance,
 815  * F is exponentially small in this region and so the absolute
 816  * accuracy is still roughly acceptable. But it would be better
 817  * to have a systematic method for improving the precision. See
 818  * the Abad+Sesma method discussion below.
 819  */
 820 static
 821 int
 822 coulomb_jwkb(const double lam, const double eta, const double x,
 823              gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,
 824              double * exponent)
 825 {
 826   const double llp1      = lam*(lam+1.0) + 6.0/35.0;
 827   const double llp1_eff  = GSL_MAX(llp1, 0.0);
 828   const double rho_ghalf = sqrt(x*(2.0*eta - x) + llp1_eff);
 829   const double sinh_arg  = sqrt(llp1_eff/(eta*eta+llp1_eff)) * rho_ghalf / x;
 830   const double sinh_inv  = log(sinh_arg + hypot(1.0,sinh_arg));
 831
 832   const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);
 833
 834   const double zeta_half = pow(3.0*phi/2.0, 1.0/3.0);
 835   const double prefactor = sqrt(M_PI*phi*x/(6.0 * rho_ghalf));
 836
 837   double F = prefactor * 3.0/zeta_half;
 838   double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */
 839   double F_exp;
 840   double G_exp;
 841
 842   const double airy_scale_exp = phi;
 843   gsl_sf_result ai;
 844   gsl_sf_result bi;
 845   gsl_sf_airy_Ai_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &ai);
 846   gsl_sf_airy_Bi_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &bi);
 847   F *= ai.val;
 848   G *= bi.val;
 849   F_exp = log(F) - airy_scale_exp;
 850   G_exp = log(G) + airy_scale_exp;
 851
 852   if(G_exp >= GSL_LOG_DBL_MAX) {
 853     fjwkb->val = F;
 854     gjwkb->val = G;
 855     fjwkb->err = 1.0e-3 * fabs(F); /* FIXME: real error here ... could be smaller */
 856     gjwkb->err = 1.0e-3 * fabs(G);
 857     *exponent = airy_scale_exp;
 858     GSL_ERROR ("error", GSL_EOVRFLW);
 859   }
 860   else {
 861     fjwkb->val = exp(F_exp);
 862     gjwkb->val = exp(G_exp);
 863     fjwkb->err = 1.0e-3 * fabs(fjwkb->val);
 864     gjwkb->err = 1.0e-3 * fabs(gjwkb->val);
 865     *exponent = 0.0;
 866     return GSL_SUCCESS;
 867   }
 868 }
 869
 870
 871 /* Asymptotic evaluation of F and G below the minimal turning point.
 872  *
 873  * This is meant to be a drop-in replacement for coulomb_jwkb().
 874  * It uses the expressions in [Abad+Sesma]. This requires some
 875  * work because I am not sure where it is valid. They mumble
 876  * something about |x| < |lam|^(-1/2) or 8|eta x| > lam when |x| < 1.
 877  * This seems true, but I thought the result was based on a uniform
 878  * expansion and could be controlled by simply using more terms.
 879  */
 880 #if 0
 881 static
 882 int
 883 coulomb_AS_xlt2eta(const double lam, const double eta, const double x,
 884                    gsl_sf_result * f_AS, gsl_sf_result * g_AS,
 885                    double * exponent)
 886 {
 887   /* no time to do this now... */
 888 }
 889 #endif /* 0 */
 890
 891
 892
 893 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 894
 895 int
 896 gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
 897                             const double lam_F,
 898                             const int  k_lam_G,      /* lam_G = lam_F - k_lam_G */
 899                             gsl_sf_result * F, gsl_sf_result * Fp,
 900                             gsl_sf_result * G, gsl_sf_result * Gp,
 901                             double * exp_F, double * exp_G)
 902 {
 903   const double lam_G = lam_F - k_lam_G;
 904
 905   if(x < 0.0 || lam_F <= -0.5 || lam_G <= -0.5) {
 906     GSL_SF_RESULT_SET(F,  0.0, 0.0);
 907     GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
 908     GSL_SF_RESULT_SET(G,  0.0, 0.0);
 909     GSL_SF_RESULT_SET(Gp, 0.0, 0.0);
 910     *exp_F = 0.0;
 911     *exp_G = 0.0;
 912     GSL_ERROR ("domain error", GSL_EDOM);
 913   }
 914   else if(x == 0.0) {
 915     gsl_sf_result C0;
 916     CLeta(0.0, eta, &C0);
 917     GSL_SF_RESULT_SET(F,  0.0, 0.0);
 918     GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
 919     GSL_SF_RESULT_SET(G,  0.0, 0.0); /* FIXME: should be Inf */
 920     GSL_SF_RESULT_SET(Gp, 0.0, 0.0); /* FIXME: should be Inf */
 921     *exp_F = 0.0;
 922     *exp_G = 0.0;
 923     if(lam_F == 0.0){
 924       GSL_SF_RESULT_SET(Fp, C0.val, C0.err);
 925     }
 926     if(lam_G == 0.0) {
 927       GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val));
 928     }
 929     GSL_ERROR ("domain error", GSL_EDOM);
 930     /* After all, since we are asking for G, this is a domain error... */
 931   }
 932   else if(x < 1.2 && 2.0*M_PI*eta < 0.9*(-GSL_LOG_DBL_MIN) && fabs(eta*x) < 10.0) {
 933     /* Reduce to a small lambda value and use the series
 934      * representations for F and G. We cannot allow eta to
 935      * be large and positive because the connection formula
 936      * for G_lam is badly behaved due to an underflow in sin(phi_lam)
 937      * [see coulomb_FG_series() and coulomb_connection() above].
 938      * Note that large negative eta is ok however.
 939      */
 940     const double SMALL = GSL_SQRT_DBL_EPSILON;
 941     const int N    = (int)(lam_F + 0.5);
 942     const int span = GSL_MAX(k_lam_G, N);
 943     const double lam_min = lam_F - N;    /* -1/2 <= lam_min < 1/2 */
 944     double F_lam_F, Fp_lam_F;
 945     double G_lam_G, Gp_lam_G;
 946     double F_lam_F_err, Fp_lam_F_err;
 947     double Fp_over_F_lam_F;
 948     double F_sign_lam_F;
 949     double F_lam_min_unnorm, Fp_lam_min_unnorm;
 950     double Fp_over_F_lam_min;
 951     gsl_sf_result F_lam_min;
 952     gsl_sf_result G_lam_min, Gp_lam_min;
 953     double F_scale;
 954     double Gerr_frac;
 955     double F_scale_frac_err;
 956     double F_unnorm_frac_err;
 957
 958     /* Determine F'/F at lam_F. */
 959     int CF1_count;
 960     int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
 961
 962     int stat_ser;
 963     int stat_Fr;
 964     int stat_Gr;
 965
 966     /* Recurse down with unnormalized F,F' values. */
 967     F_lam_F  = SMALL;
 968     Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
 969     if(span != 0) {
 970       stat_Fr = coulomb_F_recur(lam_min, span, eta, x,
 971                                 F_lam_F, Fp_lam_F,
 972                                 &F_lam_min_unnorm, &Fp_lam_min_unnorm
 973                                 );
 974     }
 975     else {
 976       F_lam_min_unnorm  =  F_lam_F;
 977       Fp_lam_min_unnorm = Fp_lam_F;
 978       stat_Fr = GSL_SUCCESS;
 979     }
 980
 981     /* Determine F and G at lam_min. */
 982     if(lam_min == -0.5) {
 983       stat_ser = coulomb_FGmhalf_series(eta, x, &F_lam_min, &G_lam_min);
 984     }
 985     else if(lam_min == 0.0) {
 986       stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min);
 987     }
 988     else if(lam_min == 0.5) {
 989       /* This cannot happen. */
 990       F->val  = F_lam_F;
 991       F->err  = 2.0 * GSL_DBL_EPSILON * fabs(F->val);
 992       Fp->val = Fp_lam_F;
 993       Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val);
 994       G->val  = G_lam_G;
 995       G->err  = 2.0 * GSL_DBL_EPSILON * fabs(G->val);
 996       Gp->val = Gp_lam_G;
 997       Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val);
 998       *exp_F = 0.0;
 999       *exp_G = 0.0;
1000       GSL_ERROR ("error", GSL_ESANITY);
1001     }
1002     else {
1003       stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min);
1004     }
1005
1006     /* Determine remaining quantities. */
1007     Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
1008     Gp_lam_min.val  = Fp_over_F_lam_min*G_lam_min.val - 1.0/F_lam_min.val;
1009     Gp_lam_min.err  = fabs(Fp_over_F_lam_min)*G_lam_min.err;
1010     Gp_lam_min.err += fabs(1.0/F_lam_min.val) * fabs(F_lam_min.err/F_lam_min.val);
1011     F_scale     = F_lam_min.val / F_lam_min_unnorm;
1012
1013     /* Apply scale to the original F,F' values. */
1014     F_scale_frac_err  = fabs(F_lam_min.err/F_lam_min.val);
1015     F_unnorm_frac_err = 2.0*GSL_DBL_EPSILON*(CF1_count+span+1);
1016     F_lam_F     *= F_scale;
1017     F_lam_F_err  = fabs(F_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
1018     Fp_lam_F    *= F_scale;
1019     Fp_lam_F_err = fabs(Fp_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
1020
1021     /* Recurse up to get the required G,G' values. */
1022     stat_Gr = coulomb_G_recur(lam_min, GSL_MAX(N-k_lam_G,0), eta, x,
1023                               G_lam_min.val, Gp_lam_min.val,
1024                               &G_lam_G, &Gp_lam_G
1025                               );
1026
1027     F->val  = F_lam_F;
1028     F->err  = F_lam_F_err;
1029     F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F);
1030
1031     Fp->val  = Fp_lam_F;
1032     Fp->err  = Fp_lam_F_err;
1033     Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F);
1034
1035     Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val);
1036
1037     G->val  = G_lam_G;
1038     G->err  = Gerr_frac * fabs(G_lam_G);
1039     G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val);
1040
1041     Gp->val  = Gp_lam_G;
1042     Gp->err  = Gerr_frac * fabs(Gp->val);
1043     Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val);
1044
1045     *exp_F = 0.0;
1046     *exp_G = 0.0;
1047
1048     return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr);
1049   }
1050   else if(x < 2.0*eta) {
1051     /* Use WKB approximation to obtain F and G at the two
1052      * lambda values, and use the Wronskian and the
1053      * continued fractions for F'/F to obtain F' and G'.
1054      */
1055     gsl_sf_result F_lam_F, G_lam_F;
1056     gsl_sf_result F_lam_G, G_lam_G;
1057     double exp_lam_F, exp_lam_G;
1058     int stat_lam_F;
1059     int stat_lam_G;
1060     int stat_CF1_lam_F;
1061     int stat_CF1_lam_G;
1062     int CF1_count;
1063     double Fp_over_F_lam_F;
1064     double Fp_over_F_lam_G;
1065     double F_sign_lam_F;
1066     double F_sign_lam_G;
1067
1068     stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F);
1069     if(k_lam_G == 0) {
1070       stat_lam_G = stat_lam_F;
1071       F_lam_G = F_lam_F;
1072       G_lam_G = G_lam_F;
1073       exp_lam_G = exp_lam_F;
1074     }
1075     else {
1076       stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G);
1077     }
1078
1079     stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
1080     if(k_lam_G == 0) {
1081       stat_CF1_lam_G  = stat_CF1_lam_F;
1082       F_sign_lam_G    = F_sign_lam_F;
1083       Fp_over_F_lam_G = Fp_over_F_lam_F;
1084     }
1085     else {
1086       stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count);
1087     }
1088
1089     F->val = F_lam_F.val;
1090     F->err = F_lam_F.err;
1091
1092     G->val = G_lam_G.val;
1093     G->err = G_lam_G.err;
1094
1095     Fp->val  = Fp_over_F_lam_F * F_lam_F.val;
1096     Fp->err  = fabs(Fp_over_F_lam_F) * F_lam_F.err;
1097     Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val);
1098
1099     Gp->val  = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val;
1100     Gp->err  = fabs(Fp_over_F_lam_G) * G_lam_G.err;
1101     Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val);
1102
1103     *exp_F = exp_lam_F;
1104     *exp_G = exp_lam_G;
1105
1106     if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) {
1107       GSL_ERROR ("overflow", GSL_EOVRFLW);
1108     }
1109     else {
1110       return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G);
1111     }
1112   }
1113   else {
1114     /* x > 2 eta, so we know that we can find a lambda value such
1115      * that x is above the turning point. We do this, evaluate
1116      * using Steed's method at that oscillatory point, then
1117      * use recursion on F and G to obtain the required values.
1118      *
1119      * lam_0   = a value of lambda such that x is below the turning point
1120      * lam_min = minimum of lam_0 and the requested lam_G, since
1121      *           we must go at least as low as lam_G
1122      */
1123     const double SMALL = GSL_SQRT_DBL_EPSILON;
1124     const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta));
1125     const int N = ceil(lam_F - C + 0.5);
1126     const double lam_0   = lam_F - GSL_MAX(N, 0);
1127     const double lam_min = GSL_MIN(lam_0, lam_G);
1128     double F_lam_F, Fp_lam_F;
1129     double G_lam_G, Gp_lam_G;
1130     double F_lam_min_unnorm, Fp_lam_min_unnorm;
1131     double F_lam_min, Fp_lam_min;
1132     double G_lam_min, Gp_lam_min;
1133     double Fp_over_F_lam_F;
1134     double Fp_over_F_lam_min;
1135     double F_sign_lam_F, F_sign_lam_min;
1136     double P_lam_min, Q_lam_min;
1137     double alpha;
1138     double gamma;
1139     double F_scale;
1140
1141     int CF1_count;
1142     int CF2_count;
1143     int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
1144     int stat_CF2;
1145     int stat_Fr;
1146     int stat_Gr;
1147
1148     int F_recur_count;
1149     int G_recur_count;
1150
1151     double err_amplify;
1152
1153     F_lam_F  = F_sign_lam_F * SMALL;  /* unnormalized */
1154     Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
1155
1156     /* Backward recurrence to get F,Fp at lam_min */
1157     F_recur_count = GSL_MAX(k_lam_G, N);
1158     stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x,
1159                               F_lam_F, Fp_lam_F,
1160                               &F_lam_min_unnorm, &Fp_lam_min_unnorm
1161                               );
1162     Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
1163
1164     /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */
1165     stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count);
1166     alpha = Fp_over_F_lam_min - P_lam_min;
1167     gamma = alpha/Q_lam_min;
1168
1169     F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ;
1170
1171     F_lam_min  = F_sign_lam_min / sqrt(alpha*alpha/Q_lam_min + Q_lam_min);
1172     Fp_lam_min = Fp_over_F_lam_min * F_lam_min;
1173     G_lam_min  = gamma * F_lam_min;
1174     Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min;
1175
1176     /* Apply scale to values of F,Fp at lam_F (the top). */
1177     F_scale = F_lam_min / F_lam_min_unnorm;
1178     F_lam_F  *= F_scale;
1179     Fp_lam_F *= F_scale;
1180
1181     /* Forward recurrence to get G,Gp at lam_G (the top). */
1182     G_recur_count = GSL_MAX(N-k_lam_G,0);
1183     stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x,
1184                               G_lam_min, Gp_lam_min,
1185                               &G_lam_G, &Gp_lam_G
1186                               );
1187
1188     err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1;
1189
1190     F->val  = F_lam_F;
1191     F->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val);
1192
1193     Fp->val = Fp_lam_F;
1194     Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val);
1195
1196     G->val  = G_lam_G;
1197     G->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val);
1198
1199     Gp->val = Gp_lam_G;
1200     Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val);
1201
1202     *exp_F = 0.0;
1203     *exp_G = 0.0;
1204
1205     return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr);
1206   }
1207 }
1208
1209
1210 int
1211 gsl_sf_coulomb_wave_F_array(double lam_min, int kmax,
1212                                  double eta, double x,
1213                                  double * fc_array,
1214                                  double * F_exp)
1215 {
1216   if(x == 0.0) {
1217     int k;
1218     *F_exp = 0.0;
1219     for(k=0; k<=kmax; k++) {
1220       fc_array[k] = 0.0;
1221     }
1222     if(lam_min == 0.0){
1223       gsl_sf_result f_0;
1224       CLeta(0.0, eta, &f_0);
1225       fc_array[0] = f_0.val;
1226     }
1227     return GSL_SUCCESS;
1228   }
1229   else {
1230     const double x_inv = 1.0/x;
1231     const double lam_max = lam_min + kmax;
1232     gsl_sf_result F, Fp;
1233     gsl_sf_result G, Gp;
1234     double G_exp;
1235
1236     int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0,
1237                                               &F, &Fp, &G, &Gp, F_exp, &G_exp);
1238
1239     double fcl  = F.val;
1240     double fpl = Fp.val;
1241     double lam = lam_max;
1242     int k;
1243
1244     fc_array[kmax] = F.val;
1245
1246     for(k=kmax-1; k>=0; k--) {
1247       double el = eta/lam;
1248       double rl = hypot(1.0, el);
1249       double sl = el  + lam*x_inv;
1250       double fc_lm1 = (fcl*sl + fpl)/rl;
1251       fc_array[k]   = fc_lm1;
1252       fpl           =  fc_lm1*sl - fcl*rl;
1253       fcl           =  fc_lm1;
1254       lam -= 1.0;
1255     }
1256
1257     return stat_FG;
1258   }
1259 }
1260
1261
1262 int
1263 gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax,
1264                                   double eta, double x,
1265                                   double * fc_array, double * gc_array,
1266                                   double * F_exp, double * G_exp)
1267 {
1268   const double x_inv = 1.0/x;
1269   const double lam_max = lam_min + kmax;
1270   gsl_sf_result F, Fp;
1271   gsl_sf_result G, Gp;
1272
1273   int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
1274                                             &F, &Fp, &G, &Gp, F_exp, G_exp);
1275
1276   double fcl  = F.val;
1277   double fpl = Fp.val;
1278   double lam = lam_max;
1279   int k;
1280
1281   double gcl, gpl;
1282
1283   fc_array[kmax] = F.val;
1284
1285   for(k=kmax-1; k>=0; k--) {
1286     double el = eta/lam;
1287     double rl = hypot(1.0, el);
1288     double sl = el  + lam*x_inv;
1289     double fc_lm1;
1290     fc_lm1 = (fcl*sl + fpl)/rl;
1291     fc_array[k] = fc_lm1;
1292     fpl         =  fc_lm1*sl - fcl*rl;
1293     fcl         =  fc_lm1;
1294     lam -= 1.0;
1295   }
1296
1297   gcl = G.val;
1298   gpl = Gp.val;
1299   lam = lam_min + 1.0;
1300
1301   gc_array[0] = G.val;
1302
1303   for(k=1; k<=kmax; k++) {
1304     double el = eta/lam;
1305     double rl = hypot(1.0, el);
1306     double sl = el + lam*x_inv;
1307     double gcl1 = (sl*gcl - gpl)/rl;
1308     gc_array[k] = gcl1;
1309     gpl         = rl*gcl - sl*gcl1;
1310     gcl         = gcl1;
1311     lam += 1.0;
1312   }
1313
1314   return stat_FG;
1315 }
1316
1317
1318 int
1319 gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax,
1320                                    double eta, double x,
1321                                    double * fc_array, double * fcp_array,
1322                                    double * gc_array, double * gcp_array,
1323                                    double * F_exp, double * G_exp)
1324
1325 {
1326   const double x_inv = 1.0/x;
1327   const double lam_max = lam_min + kmax;
1328   gsl_sf_result F, Fp;
1329   gsl_sf_result G, Gp;
1330
1331   int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
1332                                             &F, &Fp, &G, &Gp, F_exp, G_exp);
1333
1334   double fcl  = F.val;
1335   double fpl = Fp.val;
1336   double lam = lam_max;
1337   int k;
1338
1339   double gcl, gpl;
1340
1341   fc_array[kmax]  = F.val;
1342   fcp_array[kmax] = Fp.val;
1343
1344   for(k=kmax-1; k>=0; k--) {
1345     double el = eta/lam;
1346     double rl = hypot(1.0, el);
1347     double sl = el  + lam*x_inv;
1348     double fc_lm1;
1349     fc_lm1 = (fcl*sl + fpl)/rl;
1350     fc_array[k]  = fc_lm1;
1351     fpl          = fc_lm1*sl - fcl*rl;
1352     fcp_array[k] = fpl;
1353     fcl          =  fc_lm1;
1354     lam -= 1.0;
1355   }
1356
1357   gcl = G.val;
1358   gpl = Gp.val;
1359   lam = lam_min + 1.0;
1360
1361   gc_array[0]  = G.val;
1362   gcp_array[0] = Gp.val;
1363
1364   for(k=1; k<=kmax; k++) {
1365     double el = eta/lam;
1366     double rl = hypot(1.0, el);
1367     double sl = el + lam*x_inv;
1368     double gcl1 = (sl*gcl - gpl)/rl;
1369     gc_array[k]  = gcl1;
1370     gpl          = rl*gcl - sl*gcl1;
1371     gcp_array[k] = gpl;
1372     gcl          = gcl1;
1373     lam += 1.0;
1374   }
1375
1376   return stat_FG;
1377 }
1378
1379
1380 int
1381 gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
1382                                     double eta, double x,
1383                                     double * fc_array,
1384                                     double * F_exp)
1385 {
1386   if(x < 0.0 || lam_min < -0.5) {
1387     GSL_ERROR ("domain error", GSL_EDOM);
1388   }
1389   else if(x < 10.0/GSL_DBL_MAX) {
1390     int k;
1391     for(k=0; k<=kmax; k++) {
1392       fc_array[k] = 0.0;
1393     }
1394     if(lam_min == 0.0) {
1395       fc_array[0] = sqrt(C0sq(eta));
1396     }
1397     *F_exp = 0.0;
1398     if(x == 0.0)
1399       return GSL_SUCCESS;
1400     else
1401       GSL_ERROR ("underflow", GSL_EUNDRFLW);
1402   }
1403   else {
1404     int k;
1405     int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax,
1406                                                   eta, x,
1407                                                   fc_array,
1408                                                   F_exp);
1409
1410     for(k=0; k<=kmax; k++) {
1411       fc_array[k] = fc_array[k] / x;
1412     }
1413     return stat_F;
1414   }
1415 }
1416
1417