From c0bbe38a4b95aef5e4efd2ae54f6d210f905d8f9 Mon Sep 17 00:00:00 2001
From: "Knox S. Long" <long@stsci.edu>
Date: Sun, 12 May 2019 10:29:59 -0500
Subject: [PATCH] Attempt to speed up the calculation of the gaunt_factor in
 the Sutherland 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       | 29 +++++++++++++++++++++--------
 source/get_atomicdata.c |  2 +-
 source/matom.c          | 14 ++++++++++++--
 3 files changed, 34 insertions(+), 11 deletions(-)

diff --git a/source/emission.c b/source/emission.c
index ec408019f..27e205817 100644
--- a/source/emission.c
+++ b/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,6 +787,7 @@ 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
@@ -793,15 +795,26 @@ gaunt_ff (gsquared)
   {
     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);
 }
diff --git a/source/get_atomicdata.c b/source/get_atomicdata.c
index 15d74dfcb..c286f6595 100644
--- a/source/get_atomicdata.c
+++ b/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
  * # ------------------------------------------------------------
diff --git a/source/matom.c b/source/matom.c
index 554a6d04a..92425aa7c 100755
--- a/source/matom.c
+++ b/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