Attempt to speed up the calculation of the gaunt_factor in the Suther…

…land and (#516)

Basically all that has been done here is to speed up the Dopita and Sutherland gaunt factor calculation by breaking out of the loop designed to find the correct range for a particular approximation.   This is associated with issue #470

source/emission.c

@@ -768,7 +768,8 @@ one_ff (one, f1, f2)
  *
  * @details
  * It interpolates simply interpolates between the tabulated
- * factors from Table 3 of Sutherland.
+ * factors from Table 3 of Sutherland (which is a spline fit
+ * to the calculated gaunt factors at various values of gsquared)
  *
  * ### Notes ###
  * The reference for this is Sutherland, R.~S. 1998, MNRAS, 300, 321
@@ -786,22 +787,34 @@ gaunt_ff (gsquared)
   double gaunt;
   double log_g2;
   double delta;                 //The log difference between our G2 and the one in the table
+
   delta = 0.0;                  /* NSH 130605 to remove o3 compile error */
   index = 0;                    /* NSH 130605 to remove o3 compile error */
   log_g2 = log10 (gsquared);    //The data is in log format
   if (log_g2 < gaunt_total[0].log_gsqrd || log_g2 > gaunt_total[gaunt_n_gsqrd - 1].log_gsqrd)
   {
     return (1.0);
   }
-  for (i = 0; i < gaunt_n_gsqrd; i++)   /*first find the pair of parameter arrays that bracket our temperature */
-  {
-    if (gaunt_total[i].log_gsqrd <= log_g2 && gaunt_total[i + 1].log_gsqrd > log_g2)
-    {
-      index = i;                /* the array to use */
-      delta = log_g2 - gaunt_total[index].log_gsqrd;
-    }
+
+//OLD  for (i = 0; i < gaunt_n_gsqrd; i++)   /*first find the pair of parameter arrays that bracket our temperature */
+//OLD  {
+//OLD    if (gaunt_total[i].log_gsqrd <= log_g2 && gaunt_total[i + 1].log_gsqrd > log_g2)
+//OLD    {
+//OLD      index = i;                /* the array to use */
+//OLD      delta = log_g2 - gaunt_total[index].log_gsqrd;
+//OLD    }
+//OLD  }
+
+  i=0;
+  while (gaunt_total[i].log_gsqrd < log_g2) {
+      i++;
   }
 
+  index=i-1;
+  delta = log_g2 - gaunt_total[index].log_gsqrd;
+
+  /* The outherland interpolation data is a spline fit to the gaunt function. */
+
   gaunt = gaunt_total[index].gff + delta * (gaunt_total[index].s1 + delta * (gaunt_total[index].s2 + gaunt_total[index].s3));
   return (gaunt);
 }

source/get_atomicdata.c

@@ -2302,7 +2302,7 @@ would like to have simple lines for macro-ions */
  * from the data on the website, just with the top few lines commented out,
  * and a label prepended to each line
  * The data is a spline fit to the gaunt factor as a function of g - the reduced
- * temperature.
+ * temperature.  The file format is shown below; there are 45 lines
  * @verbatim
  * # log(g^2) <gff(g^2)>       s1           s2            s3
  * # ------------------------------------------------------------

source/matom.c

@@ -797,10 +797,9 @@ kpkt (p, nres, escape, mode)
      generation (see call to one_ff below) */
   if (geo.ioniz_or_extract)
   {
-    /* in spectral cycles, so use the boundaries of the photon generation bands */
+    /* in ionization cycles, so use the boundaries of the photon generation bands */
     freqmin = xband.f1[0];
 
-
     /* JM 1709 -- introduce a maximum frequency based on exp(-h nu / (kT)), 
        see issue #300 */
     freqmax = ALPHA_FF * xplasma->t_e / H_OVER_K;
@@ -1027,6 +1026,9 @@ kpkt (p, nres, escape, mode)
 
   }
 
+/* This is the end of the cooling rate calculations, which is done only once for each cell
+   and once for each cycle
+   */
 
   /* only include adiabatic cooling if we're in the right mode. First set a default 
      where adiabatic cooling is zero. This will be true if the mode isn't KPKT_MODE_ALL,
@@ -1066,6 +1068,14 @@ kpkt (p, nres, escape, mode)
 
   destruction_choice = random_number (0.0, 1.0) * cooling_normalisation;
 
+  /* ksl - This logic of what follows may not be obvious.  For choosing the basic
+   * process, we just look to see if the destruction choice is less than bf, bf+bb, bf+bb+ff
+   * etc, but inside the bhe "basic_choices", we iteratively reduce! the destruction
+   * choice until we get to one that is less than the cooling associated with a specific
+   * tranistion. If we do not find such a transition, within for example the bf if statement
+   * we drop all the way down to the Error at the end.
+   */
+
 
   if (destruction_choice < mplasma->cooling_bftot)
   {                             //destruction by bf