diff options
Diffstat (limited to 'gsl-1.9/specfunc/coulomb.c')
| -rw-r--r-- | gsl-1.9/specfunc/coulomb.c | 1417 |
1 files changed, 1417 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/coulomb.c b/gsl-1.9/specfunc/coulomb.c new file mode 100644 index 0000000..7c68076 --- /dev/null +++ b/gsl-1.9/specfunc/coulomb.c @@ -0,0 +1,1417 @@ +/* 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., 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<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); + + 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 +static +int +coulomb_AS_xlt2eta(const double lam, const double eta, const double x, + gsl_sf_result * f_AS, gsl_sf_result * g_AS, + double * exponent) +{ + /* no time to do this now... */ +} +#endif /* 0 */ + + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int +gsl_sf_coulomb_wave_FG_e(const double eta, const double x, + const double lam_F, + const int k_lam_G, /* lam_G = lam_F - k_lam_G */ + gsl_sf_result * F, gsl_sf_result * Fp, + gsl_sf_result * G, gsl_sf_result * Gp, + double * exp_F, double * exp_G) +{ + const double lam_G = lam_F - k_lam_G; + + if(x < 0.0 || lam_F <= -0.5 || lam_G <= -0.5) { + GSL_SF_RESULT_SET(F, 0.0, 0.0); + GSL_SF_RESULT_SET(Fp, 0.0, 0.0); + GSL_SF_RESULT_SET(G, 0.0, 0.0); + GSL_SF_RESULT_SET(Gp, 0.0, 0.0); + *exp_F = 0.0; + *exp_G = 0.0; + GSL_ERROR ("domain error", GSL_EDOM); + } + else if(x == 0.0) { + gsl_sf_result C0; + CLeta(0.0, eta, &C0); + GSL_SF_RESULT_SET(F, 0.0, 0.0); + GSL_SF_RESULT_SET(Fp, 0.0, 0.0); + GSL_SF_RESULT_SET(G, 0.0, 0.0); /* FIXME: should be Inf */ + GSL_SF_RESULT_SET(Gp, 0.0, 0.0); /* FIXME: should be Inf */ + *exp_F = 0.0; + *exp_G = 0.0; + if(lam_F == 0.0){ + GSL_SF_RESULT_SET(Fp, C0.val, C0.err); + } + if(lam_G == 0.0) { + GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val)); + } + GSL_ERROR ("domain error", GSL_EDOM); + /* After all, since we are asking for G, this is a domain error... */ + } + else if(x < 1.2 && 2.0*M_PI*eta < 0.9*(-GSL_LOG_DBL_MIN) && fabs(eta*x) < 10.0) { + /* Reduce to a small lambda value and use the series + * representations for F and G. We cannot allow eta to + * be large and positive because the connection formula + * for G_lam is badly behaved due to an underflow in sin(phi_lam) + * [see coulomb_FG_series() and coulomb_connection() above]. + * Note that large negative eta is ok however. + */ + const double SMALL = GSL_SQRT_DBL_EPSILON; + const int N = (int)(lam_F + 0.5); + const int span = GSL_MAX(k_lam_G, N); + const double lam_min = lam_F - N; /* -1/2 <= lam_min < 1/2 */ + double F_lam_F, Fp_lam_F; + double G_lam_G, Gp_lam_G; + double F_lam_F_err, Fp_lam_F_err; + double Fp_over_F_lam_F; + double F_sign_lam_F; + double F_lam_min_unnorm, Fp_lam_min_unnorm; + double Fp_over_F_lam_min; + gsl_sf_result F_lam_min; + gsl_sf_result G_lam_min, Gp_lam_min; + double F_scale; + double Gerr_frac; + double F_scale_frac_err; + double F_unnorm_frac_err; + + /* Determine F'/F at lam_F. */ + int CF1_count; + int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count); + + int stat_ser; + int stat_Fr; + int stat_Gr; + + /* Recurse down with unnormalized F,F' values. */ + F_lam_F = SMALL; + Fp_lam_F = Fp_over_F_lam_F * F_lam_F; + if(span != 0) { + stat_Fr = coulomb_F_recur(lam_min, span, eta, x, + F_lam_F, Fp_lam_F, + &F_lam_min_unnorm, &Fp_lam_min_unnorm + ); + } + else { + F_lam_min_unnorm = F_lam_F; + Fp_lam_min_unnorm = Fp_lam_F; + stat_Fr = GSL_SUCCESS; + } + + /* Determine F and G at lam_min. */ + if(lam_min == -0.5) { + stat_ser = coulomb_FGmhalf_series(eta, x, &F_lam_min, &G_lam_min); + } + else if(lam_min == 0.0) { + stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min); + } + else if(lam_min == 0.5) { + /* This cannot happen. */ + F->val = F_lam_F; + F->err = 2.0 * GSL_DBL_EPSILON * fabs(F->val); + Fp->val = Fp_lam_F; + Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val); + G->val = G_lam_G; + G->err = 2.0 * GSL_DBL_EPSILON * fabs(G->val); + Gp->val = Gp_lam_G; + Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val); + *exp_F = 0.0; + *exp_G = 0.0; + GSL_ERROR ("error", GSL_ESANITY); + } + else { + stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min); + } + + /* Determine remaining quantities. */ + Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm; + Gp_lam_min.val = Fp_over_F_lam_min*G_lam_min.val - 1.0/F_lam_min.val; + Gp_lam_min.err = fabs(Fp_over_F_lam_min)*G_lam_min.err; + Gp_lam_min.err += fabs(1.0/F_lam_min.val) * fabs(F_lam_min.err/F_lam_min.val); + F_scale = F_lam_min.val / F_lam_min_unnorm; + + /* Apply scale to the original F,F' values. */ + F_scale_frac_err = fabs(F_lam_min.err/F_lam_min.val); + F_unnorm_frac_err = 2.0*GSL_DBL_EPSILON*(CF1_count+span+1); + F_lam_F *= F_scale; + F_lam_F_err = fabs(F_lam_F) * (F_unnorm_frac_err + F_scale_frac_err); + Fp_lam_F *= F_scale; + Fp_lam_F_err = fabs(Fp_lam_F) * (F_unnorm_frac_err + F_scale_frac_err); + + /* Recurse up to get the required G,G' values. */ + stat_Gr = coulomb_G_recur(lam_min, GSL_MAX(N-k_lam_G,0), eta, x, + G_lam_min.val, Gp_lam_min.val, + &G_lam_G, &Gp_lam_G + ); + + F->val = F_lam_F; + F->err = F_lam_F_err; + F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F); + + Fp->val = Fp_lam_F; + Fp->err = Fp_lam_F_err; + Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F); + + Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val); + + G->val = G_lam_G; + G->err = Gerr_frac * fabs(G_lam_G); + G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val); + + Gp->val = Gp_lam_G; + Gp->err = Gerr_frac * fabs(Gp->val); + Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val); + + *exp_F = 0.0; + *exp_G = 0.0; + + return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr); + } + else if(x < 2.0*eta) { + /* Use WKB approximation to obtain F and G at the two + * lambda values, and use the Wronskian and the + * continued fractions for F'/F to obtain F' and G'. + */ + gsl_sf_result F_lam_F, G_lam_F; + gsl_sf_result F_lam_G, G_lam_G; + double exp_lam_F, exp_lam_G; + int stat_lam_F; + int stat_lam_G; + int stat_CF1_lam_F; + int stat_CF1_lam_G; + int CF1_count; + double Fp_over_F_lam_F; + double Fp_over_F_lam_G; + double F_sign_lam_F; + double F_sign_lam_G; + + stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F); + if(k_lam_G == 0) { + stat_lam_G = stat_lam_F; + F_lam_G = F_lam_F; + G_lam_G = G_lam_F; + exp_lam_G = exp_lam_F; + } + else { + stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G); + } + + stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count); + if(k_lam_G == 0) { + stat_CF1_lam_G = stat_CF1_lam_F; + F_sign_lam_G = F_sign_lam_F; + Fp_over_F_lam_G = Fp_over_F_lam_F; + } + else { + stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count); + } + + F->val = F_lam_F.val; + F->err = F_lam_F.err; + + G->val = G_lam_G.val; + G->err = G_lam_G.err; + + Fp->val = Fp_over_F_lam_F * F_lam_F.val; + Fp->err = fabs(Fp_over_F_lam_F) * F_lam_F.err; + Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val); + + Gp->val = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val; + Gp->err = fabs(Fp_over_F_lam_G) * G_lam_G.err; + Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val); + + *exp_F = exp_lam_F; + *exp_G = exp_lam_G; + + if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) { + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + else { + return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G); + } + } + else { + /* x > 2 eta, so we know that we can find a lambda value such + * that x is above the turning point. We do this, evaluate + * using Steed's method at that oscillatory point, then + * use recursion on F and G to obtain the required values. + * + * lam_0 = a value of lambda such that x is below the turning point + * lam_min = minimum of lam_0 and the requested lam_G, since + * we must go at least as low as lam_G + */ + const double SMALL = GSL_SQRT_DBL_EPSILON; + const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta)); + const int N = ceil(lam_F - C + 0.5); + const double lam_0 = lam_F - GSL_MAX(N, 0); + const double lam_min = GSL_MIN(lam_0, lam_G); + double F_lam_F, Fp_lam_F; + double G_lam_G, Gp_lam_G; + double F_lam_min_unnorm, Fp_lam_min_unnorm; + double F_lam_min, Fp_lam_min; + double G_lam_min, Gp_lam_min; + double Fp_over_F_lam_F; + double Fp_over_F_lam_min; + double F_sign_lam_F, F_sign_lam_min; + double P_lam_min, Q_lam_min; + double alpha; + double gamma; + double F_scale; + + int CF1_count; + int CF2_count; + int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count); + int stat_CF2; + int stat_Fr; + int stat_Gr; + + int F_recur_count; + int G_recur_count; + + double err_amplify; + + F_lam_F = F_sign_lam_F * SMALL; /* unnormalized */ + Fp_lam_F = Fp_over_F_lam_F * F_lam_F; + + /* Backward recurrence to get F,Fp at lam_min */ + F_recur_count = GSL_MAX(k_lam_G, N); + stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x, + F_lam_F, Fp_lam_F, + &F_lam_min_unnorm, &Fp_lam_min_unnorm + ); + Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm; + + /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */ + stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count); + alpha = Fp_over_F_lam_min - P_lam_min; + gamma = alpha/Q_lam_min; + + F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ; + + F_lam_min = F_sign_lam_min / sqrt(alpha*alpha/Q_lam_min + Q_lam_min); + Fp_lam_min = Fp_over_F_lam_min * F_lam_min; + G_lam_min = gamma * F_lam_min; + Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min; + + /* Apply scale to values of F,Fp at lam_F (the top). */ + F_scale = F_lam_min / F_lam_min_unnorm; + F_lam_F *= F_scale; + Fp_lam_F *= F_scale; + + /* Forward recurrence to get G,Gp at lam_G (the top). */ + G_recur_count = GSL_MAX(N-k_lam_G,0); + stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x, + G_lam_min, Gp_lam_min, + &G_lam_G, &Gp_lam_G + ); + + err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1; + + F->val = F_lam_F; + F->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val); + + Fp->val = Fp_lam_F; + Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val); + + G->val = G_lam_G; + G->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val); + + Gp->val = Gp_lam_G; + Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val); + + *exp_F = 0.0; + *exp_G = 0.0; + + return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr); + } +} + + +int +gsl_sf_coulomb_wave_F_array(double lam_min, int kmax, + double eta, double x, + double * fc_array, + double * F_exp) +{ + if(x == 0.0) { + int k; + *F_exp = 0.0; + for(k=0; k<=kmax; k++) { + fc_array[k] = 0.0; + } + if(lam_min == 0.0){ + gsl_sf_result f_0; + CLeta(0.0, eta, &f_0); + fc_array[0] = f_0.val; + } + return GSL_SUCCESS; + } + else { + const double x_inv = 1.0/x; + const double lam_max = lam_min + kmax; + gsl_sf_result F, Fp; + gsl_sf_result G, Gp; + double G_exp; + + int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0, + &F, &Fp, &G, &Gp, F_exp, &G_exp); + + double fcl = F.val; + double fpl = Fp.val; + double lam = lam_max; + int k; + + fc_array[kmax] = F.val; + + 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 = (fcl*sl + fpl)/rl; + fc_array[k] = fc_lm1; + fpl = fc_lm1*sl - fcl*rl; + fcl = fc_lm1; + lam -= 1.0; + } + + return stat_FG; + } +} + + +int +gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax, + double eta, double x, + double * fc_array, double * gc_array, + double * F_exp, double * G_exp) +{ + const double x_inv = 1.0/x; + const double lam_max = lam_min + kmax; + gsl_sf_result F, Fp; + gsl_sf_result G, Gp; + + int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax, + &F, &Fp, &G, &Gp, F_exp, G_exp); + + double fcl = F.val; + double fpl = Fp.val; + double lam = lam_max; + int k; + + double gcl, gpl; + + fc_array[kmax] = F.val; + + 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; + fc_array[k] = fc_lm1; + fpl = fc_lm1*sl - fcl*rl; + fcl = fc_lm1; + lam -= 1.0; + } + + gcl = G.val; + gpl = Gp.val; + lam = lam_min + 1.0; + + gc_array[0] = G.val; + + 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; + gc_array[k] = gcl1; + gpl = rl*gcl - sl*gcl1; + gcl = gcl1; + lam += 1.0; + } + + return stat_FG; +} + + +int +gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax, + double eta, double x, + double * fc_array, double * fcp_array, + double * gc_array, double * gcp_array, + double * F_exp, double * G_exp) + +{ + const double x_inv = 1.0/x; + const double lam_max = lam_min + kmax; + gsl_sf_result F, Fp; + gsl_sf_result G, Gp; + + int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax, + &F, &Fp, &G, &Gp, F_exp, G_exp); + + double fcl = F.val; + double fpl = Fp.val; + double lam = lam_max; + int k; + + double gcl, gpl; + + fc_array[kmax] = F.val; + fcp_array[kmax] = Fp.val; + + 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; + fc_array[k] = fc_lm1; + fpl = fc_lm1*sl - fcl*rl; + fcp_array[k] = fpl; + fcl = fc_lm1; + lam -= 1.0; + } + + gcl = G.val; + gpl = Gp.val; + lam = lam_min + 1.0; + + gc_array[0] = G.val; + gcp_array[0] = Gp.val; + + 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; + gc_array[k] = gcl1; + gpl = rl*gcl - sl*gcl1; + gcp_array[k] = gpl; + gcl = gcl1; + lam += 1.0; + } + + return stat_FG; +} + + +int +gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax, + double eta, double x, + double * fc_array, + double * F_exp) +{ + if(x < 0.0 || lam_min < -0.5) { + GSL_ERROR ("domain error", GSL_EDOM); + } + else if(x < 10.0/GSL_DBL_MAX) { + int k; + for(k=0; k<=kmax; k++) { + fc_array[k] = 0.0; + } + if(lam_min == 0.0) { + fc_array[0] = sqrt(C0sq(eta)); + } + *F_exp = 0.0; + if(x == 0.0) + return GSL_SUCCESS; + else + GSL_ERROR ("underflow", GSL_EUNDRFLW); + } + else { + int k; + int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax, + eta, x, + fc_array, + F_exp); + + for(k=0; k<=kmax; k++) { + fc_array[k] = fc_array[k] / x; + } + return stat_F; + } +} + + |