coulomb.c

多个常用的特殊函数的例子

  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 = sqrt(1.0 + el*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 = sqrt(1.0 + el*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 + sqrt(1.0 + sinh_arg*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;