coulomb.c

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

    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;
    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  = SMALL;
    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_lam_min  = F_sign_lam_F / 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 = sqrt(1.0 + el*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 = sqrt(1.0 + el*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 = sqrt(1.0 + el*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 = sqrt(1.0 + el*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 = sqrt(1.0 + el*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)
{
  int k;

  if(x < 0.0 || lam_min < -0.5) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(x < 10.0/GSL_DBL_MAX) {
    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;
  }
}