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 "gsl__config.h" 34 #include "gsl_math.h" 35 #include "gsl_errno.h" 36 #include "gsl_sf_exp.h" 37 #include "gsl_sf_psi.h" 38 #include "gsl_sf_airy.h" 39 #include "gsl_sf_pow_int.h" 40 #include "gsl_sf_gamma.h" 41 #include "gsl_sf_coulomb.h" 42 43 #include "gsl_specfunc__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 CF1_count; 1061 double Fp_over_F_lam_F; 1062 double Fp_over_F_lam_G; 1063 double F_sign_lam_F; 1064 1065 stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F); 1066 if(k_lam_G == 0) { 1067 stat_lam_G = stat_lam_F; 1068 F_lam_G = F_lam_F; 1069 G_lam_G = G_lam_F; 1070 exp_lam_G = exp_lam_F; 1071 } 1072 else { 1073 stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G); 1074 } 1075 1076 (void) coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count); 1077 if(k_lam_G == 0) { 1078 Fp_over_F_lam_G = Fp_over_F_lam_F; 1079 } 1080 1081 F->val = F_lam_F.val; 1082 F->err = F_lam_F.err; 1083 1084 G->val = G_lam_G.val; 1085 G->err = G_lam_G.err; 1086 1087 Fp->val = Fp_over_F_lam_F * F_lam_F.val; 1088 Fp->err = fabs(Fp_over_F_lam_F) * F_lam_F.err; 1089 Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val); 1090 1091 Gp->val = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val; 1092 Gp->err = fabs(Fp_over_F_lam_G) * G_lam_G.err; 1093 Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val); 1094 1095 *exp_F = exp_lam_F; 1096 *exp_G = exp_lam_G; 1097 1098 if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) { 1099 GSL_ERROR ("overflow", GSL_EOVRFLW); 1100 } 1101 else { 1102 return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G); 1103 } 1104 } 1105 else { 1106 /* x > 2 eta, so we know that we can find a lambda value such 1107 * that x is above the turning point. We do this, evaluate 1108 * using Steed's method at that oscillatory point, then 1109 * use recursion on F and G to obtain the required values. 1110 * 1111 * lam_0 = a value of lambda such that x is below the turning point 1112 * lam_min = minimum of lam_0 and the requested lam_G, since 1113 * we must go at least as low as lam_G 1114 */ 1115 const double SMALL = GSL_SQRT_DBL_EPSILON; 1116 const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta)); 1117 const int N = ceil(lam_F - C + 0.5); 1118 const double lam_0 = lam_F - GSL_MAX(N, 0); 1119 const double lam_min = GSL_MIN(lam_0, lam_G); 1120 double F_lam_F, Fp_lam_F; 1121 double G_lam_G, Gp_lam_G; 1122 double F_lam_min_unnorm, Fp_lam_min_unnorm; 1123 double F_lam_min; 1124 double G_lam_min, Gp_lam_min; 1125 double Fp_over_F_lam_F; 1126 double Fp_over_F_lam_min; 1127 double F_sign_lam_F, F_sign_lam_min; 1128 double P_lam_min, Q_lam_min; 1129 double alpha; 1130 double gamma; 1131 double F_scale; 1132 1133 int CF1_count; 1134 int CF2_count; 1135 int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count); 1136 int stat_CF2; 1137 int stat_Fr; 1138 int stat_Gr; 1139 1140 int F_recur_count; 1141 int G_recur_count; 1142 1143 double err_amplify; 1144 1145 F_lam_F = F_sign_lam_F * SMALL; /* unnormalized */ 1146 Fp_lam_F = Fp_over_F_lam_F * F_lam_F; 1147 1148 /* Backward recurrence to get F,Fp at lam_min */ 1149 F_recur_count = GSL_MAX(k_lam_G, N); 1150 stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x, 1151 F_lam_F, Fp_lam_F, 1152 &F_lam_min_unnorm, &Fp_lam_min_unnorm 1153 ); 1154 Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm; 1155 1156 /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */ 1157 stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count); 1158 alpha = Fp_over_F_lam_min - P_lam_min; 1159 gamma = alpha/Q_lam_min; 1160 1161 F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ; 1162 1163 F_lam_min = F_sign_lam_min / sqrt(alpha*alpha/Q_lam_min + Q_lam_min); 1164 //Fp_lam_min = Fp_over_F_lam_min * F_lam_min; 1165 G_lam_min = gamma * F_lam_min; 1166 Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min; 1167 1168 /* Apply scale to values of F,Fp at lam_F (the top). */ 1169 F_scale = F_lam_min / F_lam_min_unnorm; 1170 F_lam_F *= F_scale; 1171 Fp_lam_F *= F_scale; 1172 1173 /* Forward recurrence to get G,Gp at lam_G (the top). */ 1174 G_recur_count = GSL_MAX(N-k_lam_G,0); 1175 stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x, 1176 G_lam_min, Gp_lam_min, 1177 &G_lam_G, &Gp_lam_G 1178 ); 1179 1180 err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1; 1181 1182 F->val = F_lam_F; 1183 F->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val); 1184 1185 Fp->val = Fp_lam_F; 1186 Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val); 1187 1188 G->val = G_lam_G; 1189 G->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val); 1190 1191 Gp->val = Gp_lam_G; 1192 Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val); 1193 1194 *exp_F = 0.0; 1195 *exp_G = 0.0; 1196 1197 return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr); 1198 } 1199 } 1200 1201 1202 int 1203 gsl_sf_coulomb_wave_F_array(double lam_min, int kmax, 1204 double eta, double x, 1205 double * fc_array, 1206 double * F_exp) 1207 { 1208 if(x == 0.0) { 1209 int k; 1210 *F_exp = 0.0; 1211 for(k=0; k<=kmax; k++) { 1212 fc_array[k] = 0.0; 1213 } 1214 if(lam_min == 0.0){ 1215 gsl_sf_result f_0; 1216 CLeta(0.0, eta, &f_0); 1217 fc_array[0] = f_0.val; 1218 } 1219 return GSL_SUCCESS; 1220 } 1221 else { 1222 const double x_inv = 1.0/x; 1223 const double lam_max = lam_min + kmax; 1224 gsl_sf_result F, Fp; 1225 gsl_sf_result G, Gp; 1226 double G_exp; 1227 1228 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0, 1229 &F, &Fp, &G, &Gp, F_exp, &G_exp); 1230 1231 double fcl = F.val; 1232 double fpl = Fp.val; 1233 double lam = lam_max; 1234 int k; 1235 1236 fc_array[kmax] = F.val; 1237 1238 for(k=kmax-1; k>=0; k--) { 1239 double el = eta/lam; 1240 double rl = hypot(1.0, el); 1241 double sl = el + lam*x_inv; 1242 double fc_lm1 = (fcl*sl + fpl)/rl; 1243 fc_array[k] = fc_lm1; 1244 fpl = fc_lm1*sl - fcl*rl; 1245 fcl = fc_lm1; 1246 lam -= 1.0; 1247 } 1248 1249 return stat_FG; 1250 } 1251 } 1252 1253 1254 int 1255 gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax, 1256 double eta, double x, 1257 double * fc_array, double * gc_array, 1258 double * F_exp, double * G_exp) 1259 { 1260 const double x_inv = 1.0/x; 1261 const double lam_max = lam_min + kmax; 1262 gsl_sf_result F, Fp; 1263 gsl_sf_result G, Gp; 1264 1265 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax, 1266 &F, &Fp, &G, &Gp, F_exp, G_exp); 1267 1268 double fcl = F.val; 1269 double fpl = Fp.val; 1270 double lam = lam_max; 1271 int k; 1272 1273 double gcl, gpl; 1274 1275 fc_array[kmax] = F.val; 1276 1277 for(k=kmax-1; k>=0; k--) { 1278 double el = eta/lam; 1279 double rl = hypot(1.0, el); 1280 double sl = el + lam*x_inv; 1281 double fc_lm1; 1282 fc_lm1 = (fcl*sl + fpl)/rl; 1283 fc_array[k] = fc_lm1; 1284 fpl = fc_lm1*sl - fcl*rl; 1285 fcl = fc_lm1; 1286 lam -= 1.0; 1287 } 1288 1289 gcl = G.val; 1290 gpl = Gp.val; 1291 lam = lam_min + 1.0; 1292 1293 gc_array[0] = G.val; 1294 1295 for(k=1; k<=kmax; k++) { 1296 double el = eta/lam; 1297 double rl = hypot(1.0, el); 1298 double sl = el + lam*x_inv; 1299 double gcl1 = (sl*gcl - gpl)/rl; 1300 gc_array[k] = gcl1; 1301 gpl = rl*gcl - sl*gcl1; 1302 gcl = gcl1; 1303 lam += 1.0; 1304 } 1305 1306 return stat_FG; 1307 } 1308 1309 1310 int 1311 gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax, 1312 double eta, double x, 1313 double * fc_array, double * fcp_array, 1314 double * gc_array, double * gcp_array, 1315 double * F_exp, double * G_exp) 1316 1317 { 1318 const double x_inv = 1.0/x; 1319 const double lam_max = lam_min + kmax; 1320 gsl_sf_result F, Fp; 1321 gsl_sf_result G, Gp; 1322 1323 int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax, 1324 &F, &Fp, &G, &Gp, F_exp, G_exp); 1325 1326 double fcl = F.val; 1327 double fpl = Fp.val; 1328 double lam = lam_max; 1329 int k; 1330 1331 double gcl, gpl; 1332 1333 fc_array[kmax] = F.val; 1334 fcp_array[kmax] = Fp.val; 1335 1336 for(k=kmax-1; k>=0; k--) { 1337 double el = eta/lam; 1338 double rl = hypot(1.0, el); 1339 double sl = el + lam*x_inv; 1340 double fc_lm1; 1341 fc_lm1 = (fcl*sl + fpl)/rl; 1342 fc_array[k] = fc_lm1; 1343 fpl = fc_lm1*sl - fcl*rl; 1344 fcp_array[k] = fpl; 1345 fcl = fc_lm1; 1346 lam -= 1.0; 1347 } 1348 1349 gcl = G.val; 1350 gpl = Gp.val; 1351 lam = lam_min + 1.0; 1352 1353 gc_array[0] = G.val; 1354 gcp_array[0] = Gp.val; 1355 1356 for(k=1; k<=kmax; k++) { 1357 double el = eta/lam; 1358 double rl = hypot(1.0, el); 1359 double sl = el + lam*x_inv; 1360 double gcl1 = (sl*gcl - gpl)/rl; 1361 gc_array[k] = gcl1; 1362 gpl = rl*gcl - sl*gcl1; 1363 gcp_array[k] = gpl; 1364 gcl = gcl1; 1365 lam += 1.0; 1366 } 1367 1368 return stat_FG; 1369 } 1370 1371 1372 int 1373 gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax, 1374 double eta, double x, 1375 double * fc_array, 1376 double * F_exp) 1377 { 1378 if(x < 0.0 || lam_min < -0.5) { 1379 GSL_ERROR ("domain error", GSL_EDOM); 1380 } 1381 else if(x < 10.0/GSL_DBL_MAX) { 1382 int k; 1383 for(k=0; k<=kmax; k++) { 1384 fc_array[k] = 0.0; 1385 } 1386 if(lam_min == 0.0) { 1387 fc_array[0] = sqrt(C0sq(eta)); 1388 } 1389 *F_exp = 0.0; 1390 if(x == 0.0) 1391 return GSL_SUCCESS; 1392 else 1393 GSL_ERROR ("underflow", GSL_EUNDRFLW); 1394 } 1395 else { 1396 int k; 1397 int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax, 1398 eta, x, 1399 fc_array, 1400 F_exp); 1401 1402 for(k=0; k<=kmax; k++) { 1403 fc_array[k] = fc_array[k] / x; 1404 } 1405 return stat_F; 1406 } 1407 } 1408 1409 1410