Return to more reliable gsl integration (#1047)

This version has the correct Cloudy/Python comparisons for H and He macro atoms, and avoids the problems of #1045

source/matom.c

@@ -563,7 +563,7 @@ int temp_choice;                //choice of type of calcualation for alpha_sp
 #define ALPHA_SP_CONSTANT 5.79618e-36
 
 double
-xalpha_sp (cont_ptr, xplasma, ichoice)
+alpha_sp (cont_ptr, xplasma, ichoice)
      struct topbase_phot *cont_ptr;
      PlasmaPtr xplasma;
      int ichoice;
@@ -600,89 +600,6 @@ xalpha_sp (cont_ptr, xplasma, ichoice)
   return (alpha_sp_value);
 }
 
-
-double
-alpha_sp (cont_ptr, xplasma, ichoice)
-     struct topbase_phot *cont_ptr;
-     PlasmaPtr xplasma;
-     int ichoice;
-{
-  double alpha_sp_value, integrand;
-  double fthresh, flast;
-  int n, nmax;
-  double freq, dfreq;
-  double temp;
-  double x;
-
-  temp_choice = ichoice;
-  temp = xplasma->t_e;          //external for use in alph_sp_integrand
-  cont_ext_ptr = cont_ptr;      //"
-  fthresh = cont_ptr->freq[0];  //first frequency in list
-  flast = cont_ptr->freq[cont_ptr->np - 1];     //last frequency in list
-  if ((H_OVER_K * (flast - fthresh) / temp_ext) > ALPHA_MATOM_NUMAX_LIMIT)
-  {
-    //flast is currently very far into the exponential tail: so reduce flast to limit value of h nu / k T.
-    flast = fthresh + temp_ext * ALPHA_MATOM_NUMAX_LIMIT / H_OVER_K;
-  }
-  /* This is the line we want to replace */
-  // alpha_sp_value = num_int (alpha_sp_integrand, fthresh, flast, 1e-4);
-  alpha_sp_value = 0;
-  nmax = cont_ptr->np - 2;      // so that there is one past this poit
-
-
-  for (n = 0; n < nmax; n++)
-  {
-
-    freq = 0.5 * (cont_ptr->freq[n] + cont_ptr->freq[n + 1]);
-    x = 0.5 * (cont_ptr->x[n] + cont_ptr->x[n + 1]);
-    dfreq = cont_ptr->freq[n + 1] - cont_ptr->freq[n];
-
-    if (freq > flast)
-    {
-      break;
-    }
-
-    integrand = x * freq * freq * exp (H_OVER_K * (fthresh - freq) / temp);
-
-    if (ichoice == 1)
-    {
-      integrand *= freq / fthresh;
-
-    }
-    else if (ichoice == 2)
-    {
-      integrand *= (freq - fthresh) / fthresh;  // difference case
-    }
-
-    alpha_sp_value += integrand * dfreq;
-
-  }
-
-
-  /* This is the end of the modification */
-
-  /* The lines above evaluate the integral in alpha_sp. Now we just want to multiply
-     through by the appropriate constant. */
-  if (cont_ptr->macro_info == TRUE && geo.macro_simple == FALSE)
-  {
-    alpha_sp_value = alpha_sp_value * xconfig[cont_ptr->nlev].g / xconfig[cont_ptr->uplev].g * pow (xplasma->t_e, -1.5);
-  }
-  else                          //case for simple element
-  {
-    alpha_sp_value = alpha_sp_value * xconfig[cont_ptr->nlev].g / ion[cont_ptr->nion + 1].g * pow (xplasma->t_e, -1.5); //g for next ion up used
-  }
-
-  alpha_sp_value = alpha_sp_value * ALPHA_SP_CONSTANT;
-
-  if (sane_check (alpha_sp_value))
-  {
-    Error ("alpha_sp:sane_check value is %e\n", alpha_sp_value);
-
-  }
-
-  return (alpha_sp_value);
-}
-
 /**********************************************************/
 /**
  *  @brief the matom estimator for the band-limited spontaneous recombination rate.