3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
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.
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.
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.
20 /* Author: G. Jungman */
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.
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)]
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>
45 /* the L=0 normalization constant
46 * [Abramowitz+Stegun 14.1.8]
52 double twopieta = 2.0*M_PI*eta;
54 if(fabs(eta) < GSL_DBL_EPSILON) {
57 else if(twopieta > GSL_LOG_DBL_MAX) {
62 gsl_sf_expm1_e(twopieta, &scale);
63 return twopieta/scale.val;
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.
75 * This is not valid for L <= -3/2 or L = -1.
79 CLeta(double L, double eta, gsl_sf_result * result)
81 gsl_sf_result ln1; /* log of numerator Gamma function */
82 gsl_sf_result ln2; /* log of denominator Gamma function */
84 double arg_val, arg_err;
86 if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
87 gsl_sf_lngamma_e(L+1.0, &ln1);
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 */
94 gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
95 if(L < -1.0) sgn = -sgn;
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);
105 gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
107 /* CHECK_POINTER(result) */
110 DOMAIN_ERROR(result);
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;
119 return CLeta(lam, eta, result);
124 /* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
127 gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
131 gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
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));
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]
147 * The sequence of coefficients A_k^L is
148 * manifestly well-controlled for L >= -1/2
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.
156 * Empirically, this is quite well-behaved for
164 coulomb_Phi_series(const double lam, const double eta, const double x,
165 double * result, double * result_star)
171 double Akm1 = eta/(lam+1.0);
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;
180 for(k=2; k<kmax; k++) {
186 Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));
190 del_p = (k+lam+1.0)*del;
195 abs_del_p = fabs(del_p);
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
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.
208 prev_abs_del = abs_del;
209 prev_abs_del_p = abs_del_p;
219 GSL_ERROR ("error", GSL_EMAXITER);
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.
236 coulomb_connection(const double lam, const double eta,
237 double * cos_phi, double * sin_phi)
239 if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
242 GSL_ERROR ("error", GSL_EUNDRFLW);
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;
253 double X = tanh(M_PI * eta) / tan(M_PI * lam);
254 double phi = -atan(X) - (lam + 0.5) * M_PI;
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)
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
271 * = -atan(X) - (lam+1/2) Pi, X = tanh(eta Pi)/tan(lam Pi)
273 * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
277 coulomb_FG_series(const double lam, const double eta, const double x,
278 gsl_sf_result * F, gsl_sf_result * G)
280 const int max_iter = 800;
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);
290 double uA_mm2 = 1.0; /* uA sum is for F_{lam} */
291 double uA_mm1 = x*eta/(lam+1.0);
293 double uB_mm2 = 1.0; /* uB sum is for F_{-lam-1} */
294 double uB_mm1 = -x*eta/lam;
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;
303 int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);
305 if(stat_conn == GSL_EUNDRFLW) {
306 F->val = 0.0; /* FIXME: should this be set to Inf too like G? */
311 while(m < max_iter) {
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));
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.
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
335 A_abs_del_prev = abs_dA;
336 B_abs_del_prev = abs_dB;
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);
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);
356 GSL_ERROR ("error", GSL_EMAXITER);
358 return GSL_ERROR_SELECT_2(stat_A, stat_B);
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.
368 coulomb_FG0_series(const double eta, const double x,
369 gsl_sf_result * F, gsl_sf_result * G)
371 const int max_iter = 800;
372 const double x2 = x*x;
373 const double tex = 2.0*eta*x;
375 int stat_CL = CLeta(0.0, eta, &C0);
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 */
381 double v_mm2 = 1.0; /* nu_0 */
382 double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val); /* nu_1 */
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);
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);
393 while(m < max_iter) {
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;
403 u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
404 v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
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.
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
420 u_abs_del_prev = abs_du;
421 v_abs_del_prev = abs_dv;
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);
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);
440 GSL_ERROR ("error", GSL_EMAXITER);
442 return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
446 /* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
451 coulomb_FGmhalf_series(const double eta, const double x,
452 gsl_sf_result * F, gsl_sf_result * G)
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 */
463 double v_mm2, v_mm1, v_m;
466 gsl_sf_result rpsi_1pe;
467 gsl_sf_result rpsi_1p2e;
470 gsl_sf_psi_1piy_e(eta, &rpsi_1pe);
471 gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);
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);
476 f_sum = u_mm2 + u_mm1;
477 g_sum = v_mm2 + v_mm1;
479 while(m < max_iter) {
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;
487 && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;
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);
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);
503 GSL_ERROR ("error", GSL_EMAXITER);
509 /* Evolve the backwards recurrence for F,F'.
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)
514 * R_lam = sqrt(1 + (eta/lam)^2)
515 * S_lam = lam/x + eta/lam
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
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;
533 for(k=kmax-1; k>=0; k--) {
535 double rl = hypot(1.0, el);
536 double sl = el + lam*x_inv;
538 fc_lm1 = (fcl*sl + fpl)/rl;
539 fpl = fc_lm1*sl - fcl*rl;
550 /* Evolve the forward recurrence for G,G'.
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}
555 * where S_lam and R_lam are as above in the F recursion.
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
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;
571 for(k=1; k<=kmax; k++) {
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;
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.
595 coulomb_CF1(double lambda,
596 double eta, double x,
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;
608 double pk = lambda + 1.0;
609 double F = eta/pk + pk*x_inv;
616 if(fabs(F) < CF1_small) F = CF1_small;
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);
627 if(fabs(C) < CF1_small) C = CF1_small;
628 if(fabs(D) < CF1_small) D = CF1_small;
633 /* sign of result depends on sign of denominator */
634 *fcl_sign = - *fcl_sign;
639 GSL_ERROR ("error", GSL_ERUNAWAY);
643 while(fabs(df-1.0) > CF1_acc);
653 old_coulomb_CF1(const double lambda,
654 double eta, double x,
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;
664 double pk = lambda + 1.0;
680 F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;
682 if(fabs(eta*x + pk*pk1) > CF1_acc) break;
683 fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));
687 D = 1.0/((pk + pk1)*(x_inv + ek/pk1));
688 df = -fcl*(1.0 + ek*ek)*D;
690 if(fcl != 1.0) fcl = -1.0;
691 if(D < 0.0) fcl = -fcl;
700 tk = (pk + pk1)*(x_inv + ek/pk1);
701 D = tk - D*(1.0+ek*ek);
702 if(fabs(D) < sqrt(CF1_acc)) {
705 printf("HELP............\n");
710 /* sign of result depends on sign of denominator */
713 df = df*(D*tk - 1.0);
716 GSL_ERROR ("error", GSL_ERUNAWAY);
719 while(fabs(df) > fabs(F)*CF1_acc);
728 /* Evaluate the second continued fraction to
730 * (G' + i F')/(G + i F) := P + i Q
731 * at the specified lambda value.
735 coulomb_CF2(const double lambda, const double eta, const double x,
736 double * result_P, double * result_Q, int * count
739 int status = GSL_SUCCESS;
741 const double CF2_acc = 4.0*GSL_DBL_EPSILON;
742 const double CF2_abort = 2.0e+05;
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);
751 double br = 2.0*(x - eta);
754 double dr = br/(br*br + bi*bi);
755 double di = -bi/(br*br + bi*bi);
757 double dp = -x_inv*(ar*di + ai*dr);
758 double dq = x_inv*(ar*dr - ai*di);
764 double Q = 1.0 - eta*x_inv;
775 D = ar*dr - ai*di + br;
776 di = ai*dr + ar*di + bi;
777 C = 1.0/(D*D + di*di);
780 A = br*dr - bi*di - 1.;
786 status = GSL_ERUNAWAY;
791 while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);
793 if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {
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:
807 * result_F = fjwkb * exp(-exponent)
808 * result_G = gjwkb * exp( exponent)
810 * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]
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.
822 coulomb_jwkb(const double lam, const double eta, const double x,
823 gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,
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));
832 const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);
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));
837 double F = prefactor * 3.0/zeta_half;
838 double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */
842 const double airy_scale_exp = phi;
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);
849 F_exp = log(F) - airy_scale_exp;
850 G_exp = log(G) + airy_scale_exp;
852 if(G_exp >= GSL_LOG_DBL_MAX) {
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);
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);
871 /* Asymptotic evaluation of F and G below the minimal turning point.
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.
883 coulomb_AS_xlt2eta(const double lam, const double eta, const double x,
884 gsl_sf_result * f_AS, gsl_sf_result * g_AS,
887 /* no time to do this now... */
893 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
896 gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
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)
903 const double lam_G = lam_F - k_lam_G;
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);
912 GSL_ERROR ("domain error", GSL_EDOM);
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 */
924 GSL_SF_RESULT_SET(Fp, C0.val, C0.err);
927 GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val));
929 GSL_ERROR ("domain error", GSL_EDOM);
930 /* After all, since we are asking for G, this is a domain error... */
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.
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;
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;
955 double F_scale_frac_err;
956 double F_unnorm_frac_err;
958 /* Determine F'/F at lam_F. */
960 int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
966 /* Recurse down with unnormalized F,F' values. */
968 Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
970 stat_Fr = coulomb_F_recur(lam_min, span, eta, x,
972 &F_lam_min_unnorm, &Fp_lam_min_unnorm
976 F_lam_min_unnorm = F_lam_F;
977 Fp_lam_min_unnorm = Fp_lam_F;
978 stat_Fr = GSL_SUCCESS;
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);
985 else if(lam_min == 0.0) {
986 stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min);
988 else if(lam_min == 0.5) {
989 /* This cannot happen. */
991 F->err = 2.0 * GSL_DBL_EPSILON * fabs(F->val);
993 Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val);
995 G->err = 2.0 * GSL_DBL_EPSILON * fabs(G->val);
997 Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val);
1000 GSL_ERROR ("error", GSL_ESANITY);
1003 stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min);
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;
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);
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);
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,
1028 F->err = F_lam_F_err;
1029 F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F);
1032 Fp->err = Fp_lam_F_err;
1033 Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F);
1035 Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val);
1038 G->err = Gerr_frac * fabs(G_lam_G);
1039 G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val);
1042 Gp->err = Gerr_frac * fabs(Gp->val);
1043 Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val);
1048 return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr);
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'.
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;
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;
1068 stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F);
1070 stat_lam_G = stat_lam_F;
1073 exp_lam_G = exp_lam_F;
1076 stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G);
1079 stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
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;
1086 stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count);
1089 F->val = F_lam_F.val;
1090 F->err = F_lam_F.err;
1092 G->val = G_lam_G.val;
1093 G->err = G_lam_G.err;
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);
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);
1106 if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) {
1107 GSL_ERROR ("overflow", GSL_EOVRFLW);
1110 return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G);
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.
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
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;
1143 int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
1153 F_lam_F = F_sign_lam_F * SMALL; /* unnormalized */
1154 Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
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,
1160 &F_lam_min_unnorm, &Fp_lam_min_unnorm
1162 Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
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;
1169 F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ;
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;
1176 /* Apply scale to values of F,Fp at lam_F (the top). */
1177 F_scale = F_lam_min / F_lam_min_unnorm;
1179 Fp_lam_F *= F_scale;
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,
1188 err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1;
1191 F->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val);
1194 Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val);
1197 G->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val);
1200 Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val);
1205 return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr);
1211 gsl_sf_coulomb_wave_F_array(double lam_min, int kmax,
1212 double eta, double x,
1219 for(k=0; k<=kmax; k++) {
1224 CLeta(0.0, eta, &f_0);
1225 fc_array[0] = f_0.val;
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;
1236 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0,
1237 &F, &Fp, &G, &Gp, F_exp, &G_exp);
1240 double fpl = Fp.val;
1241 double lam = lam_max;
1244 fc_array[kmax] = F.val;
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;
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)
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;
1273 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
1274 &F, &Fp, &G, &Gp, F_exp, G_exp);
1277 double fpl = Fp.val;
1278 double lam = lam_max;
1283 fc_array[kmax] = F.val;
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;
1290 fc_lm1 = (fcl*sl + fpl)/rl;
1291 fc_array[k] = fc_lm1;
1292 fpl = fc_lm1*sl - fcl*rl;
1299 lam = lam_min + 1.0;
1301 gc_array[0] = G.val;
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;
1309 gpl = rl*gcl - sl*gcl1;
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)
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;
1331 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
1332 &F, &Fp, &G, &Gp, F_exp, G_exp);
1335 double fpl = Fp.val;
1336 double lam = lam_max;
1341 fc_array[kmax] = F.val;
1342 fcp_array[kmax] = Fp.val;
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;
1349 fc_lm1 = (fcl*sl + fpl)/rl;
1350 fc_array[k] = fc_lm1;
1351 fpl = fc_lm1*sl - fcl*rl;
1359 lam = lam_min + 1.0;
1361 gc_array[0] = G.val;
1362 gcp_array[0] = Gp.val;
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;
1370 gpl = rl*gcl - sl*gcl1;
1381 gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
1382 double eta, double x,
1386 if(x < 0.0 || lam_min < -0.5) {
1387 GSL_ERROR ("domain error", GSL_EDOM);
1389 else if(x < 10.0/GSL_DBL_MAX) {
1391 for(k=0; k<=kmax; k++) {
1394 if(lam_min == 0.0) {
1395 fc_array[0] = sqrt(C0sq(eta));
1401 GSL_ERROR ("underflow", GSL_EUNDRFLW);
1405 int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax,
1410 for(k=0; k<=kmax; k++) {
1411 fc_array[k] = fc_array[k] / x;