k-packet changes for free-free cooling

source/matom.c

@@ -795,29 +795,16 @@ kpkt (p, nres, escape, mode)
 
   /* JM 1511 -- Fix for issue 187. We need band limits for free free packet
      generation (see call to one_ff below) */
-  if (geo.ioniz_or_extract)
-  {
-    /* 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;
+  freqmin = xband.f1[0];
+  freqmax = ALPHA_FF * xplasma->t_e / H_OVER_K;
 
-    /* ksl This is a Bandaid for when the temperatures are very low */
-    if (freqmax < 1.1 * freqmin)
-    {
-      freqmax = 1.1 * freqmin;
-    }
-  }
-  else
+  /* ksl This is a Bandaid for when the temperatures are very low */
+  /* in this case cooling_ff should be low compared to cooling_ff_lofreq anyway */
+  if (freqmax < 1.1 * freqmin)
   {
-    /* in spectral cycles, use the frequency range of the final spectrum */
-    freqmin = geo.sfmin;
-    freqmax = geo.sfmax;
+    freqmax = 1.1 * freqmin;
   }
 
-
   /* ksl 091108 - If the kpkt destruction rates for this cell are not known they are calculated here.  This happens
    * every time the wind is updated */
 
@@ -943,7 +930,8 @@ kpkt (p, nres, escape, mode)
        volume.  Recall however that vol is part of the windPtr */
     if (one->vol > 0)
     {
-      cooling_ff = mplasma->cooling_ff = total_free (one, xplasma->t_e, 0.0, VERY_BIG) / xplasma->vol / xplasma->ne;    // JM 1411 - changed to use filled volume
+      cooling_ff = mplasma->cooling_ff = total_free (one, xplasma->t_e, freqmin, freqmax) / xplasma->vol / xplasma->ne; // JM 1411 - changed to use filled volume
+      cooling_ff += mplasma->cooling_ff_lofreq = total_free (one, xplasma->t_e, 0.0, freqmin) / xplasma->vol / xplasma->ne;
     }
     else
     {
@@ -958,7 +946,7 @@ kpkt (p, nres, escape, mode)
          removed it? We delete this whole "else" if we're sure
          volumes are never zero. */
 
-      cooling_ff = mplasma->cooling_ff = 0.0;
+      cooling_ff = mplasma->cooling_ff = mplasma->cooling_ff_lofreq = 0.0;
       Error ("kpkt: A scattering event in cell %d with vol = 0???\n", one->nwind);
       //Diagnostic      return(-1);  //57g -- Cannot diagnose with an exit
       *escape = 1;
@@ -1014,9 +1002,6 @@ kpkt (p, nres, escape, mode)
 
     cooling_normalisation += cooling_adiabatic;
 
-
-
-
     mplasma->cooling_bbtot = cooling_bbtot;
     mplasma->cooling_bftot = cooling_bftot;
     mplasma->cooling_bf_coltot = cooling_bf_coltot;
@@ -1176,28 +1161,30 @@ kpkt (p, nres, escape, mode)
       }
     }
   }
+
+  /* consult issues #187, #492 regarding free-free */
   else if (destruction_choice < (mplasma->cooling_bftot + mplasma->cooling_bbtot + mplasma->cooling_ff))
   {                             //this is a ff destruction
-
-    /* The limits for ff emission are hard-wired: 40 microns ->
-       twice energy of He II edge. Shouldn't be a problem unless
-       we're in plasma with temperatures that gives significant ff
-       emission outside this range. */
-
     *escape = 1;                //we are making an r-packet not exciting a macro atom
-
     *nres = -2;
-
-    /* used to do one_ff (one, 7.5e12, 2.626e16) here,
-       but now use the band boundaries, see #187. */
     p->freq = one_ff (one, freqmin, freqmax);   //get frequency of resulting energy packet
-
     return (0);
   }
+  else if (destruction_choice < (mplasma->cooling_bftot + mplasma->cooling_bbtot + mplasma->cooling_ff + mplasma->cooling_ff_lofreq))
+  {                             //this is ff at low frequency
+    *escape = 1;
+    *nres = -2;
+    /* we don't bother tracking photons below 1e14 Hz, 
+       so record that this photon was lost to "low frequency free-free" */
+    p->istat = P_LOFREQ_FF;
+    return (0);
+  }
+
 
 
   /* JM 1310 -- added loop to check if destruction occurs via adiabatic cooling */
-  else if (destruction_choice < (mplasma->cooling_bftot + mplasma->cooling_bbtot + mplasma->cooling_ff + cooling_adiabatic))
+  else if (destruction_choice <
+           (mplasma->cooling_bftot + mplasma->cooling_bbtot + mplasma->cooling_ff + mplasma->cooling_ff_lofreq + cooling_adiabatic))
   {
 
     if (geo.adiabatic == 0 || mode != KPKT_MODE_ALL)
@@ -1217,7 +1204,9 @@ kpkt (p, nres, escape, mode)
   else
   {
     /* We want destruction by collisional ionization in a macro atom. */
-    destruction_choice = destruction_choice - mplasma->cooling_bftot - mplasma->cooling_bbtot - mplasma->cooling_ff - cooling_adiabatic;
+    destruction_choice =
+      destruction_choice - mplasma->cooling_bftot - mplasma->cooling_bbtot - mplasma->cooling_ff - mplasma->cooling_ff_lofreq -
+      cooling_adiabatic;
 
     for (i = 0; i < nphot_total; i++)
     {
@@ -1255,8 +1244,9 @@ kpkt (p, nres, escape, mode)
 
   Error ("matom.c: Failed to select a destruction process in kpkt. Abort.\n");
   Error
-    ("matom.c: cooling_bftot %g, cooling_bbtot %g, cooling_ff %g, cooling_bf_coltot %g cooling_adiabatic %g\n",
-     mplasma->cooling_bftot, mplasma->cooling_bbtot, mplasma->cooling_ff, mplasma->cooling_bf_coltot, mplasma->cooling_adiabatic);
+    ("matom.c: choice %8.4e norm %8.4e cooling_bftot %g, cooling_bbtot %g, cooling_ff %g, cooling_ff_lofreq %g, cooling_bf_coltot %g cooling_adiabatic %g cooling_adiabatic %g\n",
+     destruction_choice, cooling_normalisation, mplasma->cooling_bftot, mplasma->cooling_bbtot, mplasma->cooling_ff,
+     mplasma->cooling_ff_lofreq, mplasma->cooling_bf_coltot, mplasma->cooling_adiabatic, cooling_adiabatic);
 
   *escape = 1;
   p->istat = P_ERROR_MATOM;

source/matrix_ion.c

@@ -308,15 +308,19 @@ matrix_ion_populations (xplasma, mode)
         newden[nn] = xplasma->density[nn] / elem_dens[ion[nn].z];
       }
 
-      for (mm = 0; mm < nrows; mm++)    // inner loop over the elements of the population array
+      /* if the ion is "simple" then find it's calculated ionization state in populations array */
+      else
       {
-        if (xion[mm] == nn)     // if this element contains the population of the ion is question
+
+        for (mm = 0; mm < nrows; mm++)  // inner loop over the elements of the population array
         {
-          newden[nn] = populations[mm]; // get the population
+          if (xion[mm] == nn)   // if this element contains the population of the ion is question
+          {
+            newden[nn] = populations[mm];       // get the population
+          }
         }
       }
 
-
       if (newden[nn] < DENSITY_MIN)     // this wil also capture the case where population doesnt have a value for this ion
         newden[nn] = DENSITY_MIN;
     }

source/para_update.c

@@ -376,15 +376,15 @@ communicate_matom_estimators_para ()
   gamma_helper = calloc (sizeof (double), NPLASMA * 4 * size_gamma_est);
   alpha_helper = calloc (sizeof (double), NPLASMA * 2 * size_alpha_est);
   level_helper = calloc (sizeof (double), NPLASMA * nlevels_macro);
-  cell_helper = calloc (sizeof (double), 7 * NPLASMA);
+  cell_helper = calloc (sizeof (double), 8 * NPLASMA);
   cooling_bf_helper = calloc (sizeof (double), NPLASMA * 2 * nphot_total);
   cooling_bb_helper = calloc (sizeof (double), NPLASMA * nlines);
 
   jbar_helper2 = calloc (sizeof (double), NPLASMA * size_Jbar_est);
   gamma_helper2 = calloc (sizeof (double), NPLASMA * 4 * size_gamma_est);
   alpha_helper2 = calloc (sizeof (double), NPLASMA * 2 * size_alpha_est);
   level_helper2 = calloc (sizeof (double), NPLASMA * nlevels_macro);
-  cell_helper2 = calloc (sizeof (double), 7 * NPLASMA);
+  cell_helper2 = calloc (sizeof (double), 8 * NPLASMA);
   cooling_bf_helper2 = calloc (sizeof (double), NPLASMA * 2 * nphot_total);
   cooling_bb_helper2 = calloc (sizeof (double), NPLASMA * nlines);
 
@@ -407,7 +407,8 @@ communicate_matom_estimators_para ()
     cell_helper[mpi_i + 3 * NPLASMA] = macromain[mpi_i].cooling_bf_coltot / np_mpi_global;
     cell_helper[mpi_i + 4 * NPLASMA] = macromain[mpi_i].cooling_bbtot / np_mpi_global;
     cell_helper[mpi_i + 5 * NPLASMA] = macromain[mpi_i].cooling_ff / np_mpi_global;
-    cell_helper[mpi_i + 6 * NPLASMA] = macromain[mpi_i].cooling_adiabatic / np_mpi_global;
+    cell_helper[mpi_i + 6 * NPLASMA] = macromain[mpi_i].cooling_ff_lofreq / np_mpi_global;
+    cell_helper[mpi_i + 7 * NPLASMA] = macromain[mpi_i].cooling_adiabatic / np_mpi_global;
 
 
 
@@ -490,7 +491,8 @@ communicate_matom_estimators_para ()
     macromain[mpi_i].cooling_bf_coltot = cell_helper2[mpi_i + 3 * NPLASMA];
     macromain[mpi_i].cooling_bbtot = cell_helper2[mpi_i + 4 * NPLASMA];
     macromain[mpi_i].cooling_ff = cell_helper2[mpi_i + 5 * NPLASMA];
-    macromain[mpi_i].cooling_adiabatic = cell_helper2[mpi_i + 6 * NPLASMA];
+    macromain[mpi_i].cooling_ff_lofreq = cell_helper2[mpi_i + 6 * NPLASMA];
+    macromain[mpi_i].cooling_adiabatic = cell_helper2[mpi_i + 7 * NPLASMA];
 
 
     for (n = 0; n < nlevels_macro; n++)

source/photo_gen_matom.c

@@ -327,7 +327,8 @@ get_matom_f (mode)
 
                 if (nres > NLINES + nphot_total)
                 {
-                  Error ("Problem in get_matom_f (1). Abort. \n");
+                  Error ("Problem in get_matom_f (1). Abort. nres is %d, NLINES %d, nphot_total %d m %d %8.4e\n",
+                         nres, NLINES, nphot_total, m, macromain[n].matom_abs[m]);
                   Exit (0);
                 }
 

source/python.h

@@ -1002,7 +1002,7 @@ typedef struct macro
      and used to select destruction rates for kpkts */
   double cooling_normalisation;
   double cooling_bbtot, cooling_bftot, cooling_bf_coltot;
-  double cooling_ff;
+  double cooling_ff, cooling_ff_lofreq;
   double cooling_adiabatic;     // this is just cool_adiabatic / vol / ne
 
 
@@ -1071,7 +1071,8 @@ typedef struct photon
     P_HIT_DISK = 7,             //Banged into disk
     P_SEC = 8,                  //Photon hit secondary
     P_ADIABATIC = 9,            //records that a photon created a kpkt which was destroyed by adiabatic cooling
-    P_ERROR_MATOM = 10          //Some kind of error in processing of a photon which excited a macroattom
+    P_ERROR_MATOM = 10,         //Some kind of error in processing of a photon which excited a macroattom
+    P_LOFREQ_FF = 11            //records a photon that had too low a frequency  
   } istat;                      /*status of photon. */
 
   int nscat;                    /*number of scatterings */

source/run.c

@@ -54,10 +54,10 @@ calculate_ionization (restart_stat)
      int restart_stat;
 {
   int n, nn;
-  double zz, zzz, zze, ztot, zz_adiab;
+  double zz, zzz, zze, ztot, zz_adiab, zz_lofreq;
   double zz_abs, zz_scat, zz_star, zz_disk;
   double zz_err, zz_else;
-  int nn_adiab;
+  int nn_adiab, nn_lofreq;
   WindPtr w;
   PhotPtr p;
 
@@ -227,8 +227,8 @@ calculate_ionization (restart_stat)
     trans_phot (w, p, 0);
 
     /*Determine how much energy was absorbed in the wind */
-    zze = zzz = zz_adiab = zz_abs = zz_scat = zz_star = zz_disk = zz_err = zz_else = 0.0;
-    nn_adiab = 0;
+    zze = zzz = zz_adiab = zz_abs = zz_scat = zz_star = zz_disk = zz_err = zz_else = zz_lofreq = 0.0;
+    nn_adiab = nn_lofreq = 0;
     for (nn = 0; nn < NPHOT; nn++)
     {
       zzz += p[nn].w;
@@ -241,6 +241,11 @@ calculate_ionization (restart_stat)
         zz_adiab += p[nn].w;
         nn_adiab++;
       }
+      else if (p[nn].istat == P_LOFREQ_FF)
+      {
+        zz_lofreq += p[nn].w;
+        nn_lofreq++;
+      }
       else if (p[nn].istat == P_ABSORB)
       {
         zz_abs += p[nn].w;
@@ -271,7 +276,10 @@ calculate_ionization (restart_stat)
       ("!!python: Total photon luminosity after transphot  %18.12e (absorbed/lost  %18.12e). Radiated luminosity %18.12e\n",
        zzz, zzz - zz, zze);
     if (geo.rt_mode == RT_MODE_MACRO)
+    {
       Log ("!!python: luminosity lost by adiabatic kpkt destruction %18.12e number of packets %d\n", zz_adiab, nn_adiab);
+      Log ("!!python: luminosity lost to low-frequency free-free    %18.12e number of packets %d\n", zz_lofreq, nn_lofreq);
+    }
     Log ("!!python: luminosity lost by being completely absorbed  %18.12e \n", zz_abs);
     Log ("!!python: luminosity lost by too many scatters          %18.12e \n", zz_scat);
     Log ("!!python: luminosity lost by hitting the star           %18.12e \n", zz_star);
@@ -492,6 +500,15 @@ make_spectra (restart_stat)
 
   kbf_need (freqmin, freqmax);
 
+  /* force recalculation of kpacket rates */
+  if (geo.rt_mode == RT_MODE_MACRO) 
+  {
+    for (n == 0; n < NPLASMA; n++)
+    {
+      macromain[n].kpkt_rates_known = -1;
+    }
+  }
+
   /* BEGIN CYCLES TO CREATE THE DETAILED SPECTRUM */
 
   /* the next section initializes the spectrum array in two cases, for the

source/trans_phot.c

@@ -616,16 +616,9 @@ trans_phot_single (WindPtr w, PhotPtr p, int iextract)
       break;
     }
 
-    if (pp.istat == P_ADIABATIC)
+    if (pp.istat == P_ERROR_MATOM || pp.istat == P_LOFREQ_FF || pp.istat == P_ADIABATIC)
     {
-      istat = pp.istat = p->istat = P_ADIABATIC;
-      stuff_phot (&pp, p);
-      break;
-    }
-
-    if (pp.istat == P_ERROR_MATOM)
-    {
-      istat = pp.istat = p->istat = P_ERROR_MATOM;
+      istat = p->istat = pp.istat;
       stuff_phot (&pp, p);
       break;
     }

source/wind2d.c

@@ -888,9 +888,6 @@ wind_div_v ()
     div += xxx[2] = (v2[2] - v1[2]) / delta;
 
 
-//      printf ("BLAH cell %i div=%e\n",icell,div);
-
-
     /* we have now evaluated the divergence, so can store in the wind pointer */
     wmain[icell].div_v = div;