Proposed minimal banding change to prevent fewer than 7 ionization bands

source/bands.c

@@ -68,8 +68,6 @@ struct xbands
 	int nbands;           // Actual number of bands in use
 }
 xband;
-
-
 */
 
 
@@ -600,18 +598,13 @@ bands_init (imode, band)
          band->f1[nband] * PLANCK / (BOLTZMANN * tmax), band->f2[nband] * PLANCK / (BOLTZMANN * tmax), band->min_fraction[nband]);
   }
 
+  check_appropriate_banding (band, mode);
+
   /* we used to call freqs init here, but now the photon generation bands are 
      also tied to the ionization bands, see e.g. gh issue #1084 */
-
-  for (nband = 0; nband < band->nbands; nband++)
-  {
-    geo.xfreq[nband] = band->f1[nband];
-  }
-  geo.nxfreq = band->nbands;
-  geo.xfreq[band->nbands] = band->f2[band->nbands - 1];
+  ion_bands_init (mode, band->f1[0], band->f2[band->nbands - 1], band);
 
   /* Now define the freqquency boundaries for the cell spectra */
-
   geo.cell_log_freq_min = log10 (band->f1[0]);
   geo.cell_log_freq_max = log10 (band->f2[band->nbands - 1]);
   geo.cell_delta_lfreq = (geo.cell_log_freq_max - geo.cell_log_freq_min) / NBINS_IN_CELL_SPEC;
@@ -626,9 +619,11 @@ bands_init (imode, band)
 /** 
  * @brief      This is the routine where the frequency
  * 	binning for coarse spectra in each plasma cell is established
- *
+ * 
+ * @param [in] int  mode   The banding mode
  * @param [in] double  freqmin   The minimum frequency
  * @param [in] double  freqmax   The maximum frequency
+ * @param [in] struct xbands *band bands structure to check
  * @return     Always returns 0.  
  *
  * The frequency intervals are stored
@@ -661,14 +656,17 @@ bands_init (imode, band)
  * used to create photons.  
  *
  **********************************************************/
+#define MIN_N_IONBANDS 7
 
 int
-freqs_init (freqmin, freqmax)
+ion_bands_init (mode, freqmin, freqmax, band)
+     int mode;
      double freqmin, freqmax;
+     struct xbands *band;
 {
   int i, n, ngood, good[NXBANDS];
   double xfreq[NXBANDS];
-  int nxfreq;
+  int nxfreq, nband;
 
 
 /* At present set up a single energy band for 2 - 10 keV */
@@ -683,18 +681,27 @@ freqs_init (freqmin, freqmax)
  xfreq[6] = 10000 / HEV;
  xfreq[7] = 50000 / HEV;*/
 
-/* bands to match the cloudy table spectrum - needed to cover all frequencies to let induced compton work OK */
-
-  if (geo.agn_ion_spectype == SPECTYPE_CL_TAB)
+  /* bands to match the cloudy table spectrum - needed to cover all frequencies to let induced compton work OK */
+  if (mode == CLOUDY_TEST_BAND)
   {
     nxfreq = 3;
     xfreq[0] = 0.0001 / HEV;
     xfreq[1] = geo.agn_cltab_low / HEV;
     xfreq[2] = geo.agn_cltab_hi / HEV;
     xfreq[3] = 100000000 / HEV;
   }
+  /* if we have a sufficient number of user-selected bands, then just use those band boundaries */
+  else if (band->nbands >= MIN_N_IONBANDS)
+  {
+    for (nband = 0; nband < band->nbands; nband++)
+    {
+      xfreq[nband] = band->f1[nband];
+    }
+    nxfreq = band->nbands;
+    xfreq[band->nbands] = band->f2[band->nbands - 1];
+  }
 
-/* bands try to deal with a blackbody spectrum */
+  /* if we don't have enough bands, then use the old hardwired bands */
   else
   {
     nxfreq = 10;
@@ -709,13 +716,13 @@ freqs_init (freqmin, freqmax)
     xfreq[8] = 3.162e18;
     xfreq[9] = 1.2e19;          //This is the highest frequency defined in our ionization data
     xfreq[10] = freqmax;
-  }
-
-
-
 
+    Log ("ion_bands_init: only %d generation bands, so adopting default 10 ionization bands from %8.4e Hz to %8.4e Hz\n", band->nbands,
+         freqmin, freqmax);
+  }
 
-  Log ("freqs_init: Photons will be generated between %8.2f (%8.2e) and %8.2f (%8.2e)\n", freqmin * HEV, freqmin, freqmax * HEV, freqmax);
+  Log ("ion_bands_init: %d ionization bands from %8.2feV (%8.2eHz) to %8.2feV (%8.2eHz)\n", nxfreq, freqmin * HEV, freqmin, freqmax * HEV,
+       freqmax);
 
   ngood = 0;
   for (i = 0; i < nxfreq; i++)
@@ -736,7 +743,7 @@ freqs_init (freqmin, freqmax)
     }
   }
 
-  Log_silent ("freqs_init: Of %d starting intervals, %d will have photons\n", nxfreq, ngood);
+  Log_silent ("ion_bands_init: Of %d starting intervals, %d will have photons\n", nxfreq, ngood);
 
   n = 0;
   for (i = 0; i < nxfreq; i++)
@@ -750,10 +757,6 @@ freqs_init (freqmin, freqmax)
   }
   geo.nxfreq = n;
 
-
-
-
-
   /* OK at this point we know at least some photons will be generated in each interval, but we still don't know
    * that the we are going to have a possibilty of photons throughout the first and last intervals.
    */
@@ -769,10 +772,10 @@ freqs_init (freqmin, freqmax)
   }
 
 
-  Log_silent ("freqs_init: There were %d final intervals\n", geo.nxfreq);
+  Log_silent ("ion_bands_init: There were %d final intervals\n", geo.nxfreq);
   for (n = 0; n < geo.nxfreq; n++)
   {
-    Log_silent ("freqs_init: %8.2f (%8.2e)    %8.2f (%8.2e)  \n", geo.xfreq[n] * HEV, geo.xfreq[n], geo.xfreq[n + 1] * HEV,
+    Log_silent ("ion_bands_init: %8.2f (%8.2e)    %8.2f (%8.2e)  \n", geo.xfreq[n] * HEV, geo.xfreq[n], geo.xfreq[n + 1] * HEV,
                 geo.xfreq[n + 1]);
   }
 
@@ -781,3 +784,51 @@ freqs_init (freqmin, freqmax)
   return (0);
 
 }
+
+/**********************************************************/
+/** 
+ * @brief Helper routine for checking the photon generation
+ * banding is reasonable. At this stage fairly rudimentary
+ *
+ * @param [in] struct xbands *band bands structure to check
+ * @param [in] int  mode   The banding mode
+ *
+ * @details
+ * 
+ **********************************************************/
+
+void
+check_appropriate_banding (band, mode)
+     struct xbands *band;
+     int mode;
+{
+  if (geo.system_type == SYSTEM_TYPE_AGN)
+  {
+    if (mode == CV_BAND)
+      Error ("Using CV banding for AGN system. Not recommended!\n");
+    else if (mode == YSO_BAND)
+      Error ("Using YSO banding for AGN system. Not recommended!\n");
+    else if (mode == T_STAR_BAND)
+      Error ("Using Tstar banding for AGN system. Not recommended!\n");
+    else if (band->nbands < 4)
+      Error ("You only have %d photon generation bands for AGN system. Not recommended!\n");
+  }
+  else if (geo.system_type == SYSTEM_TYPE_CV)
+  {
+    if (mode == YSO_BAND)
+      Error ("Using YSO banding for CV system. Not recommended!\n");
+    else if (band->nbands < 4)
+      Error ("You only have %d photon generation bands for CV system. Not recommended!\n");
+  }
+  else if (geo.system_type == SYSTEM_TYPE_BH)
+  {
+    if (mode == CV_BAND)
+      Error ("Using CV banding for BH system. Not recommended!\n");
+    else if (mode == YSO_BAND)
+      Error ("Using YSO banding for BH system. Not recommended!\n");
+    else if (mode == T_STAR_BAND)
+      Error ("Using Tstar banding for BH system. Not recommended!\n");
+    else if (band->nbands < 4)
+      Error ("You only have %d photon generation bands for BH system. Not recommended!\n");
+  }
+}

source/parse.c

@@ -154,13 +154,13 @@ parse_command_line (argc, argv)
         Log ("Using special relativity for Doppler shifts, etc., and including co-moving frame effects.\n");
         j = i;
       }
-      else if (strcmp (argv[i], "-srclassic") == 0)
+      else if (strcmp (argv[i], "-sr_doppler_only") == 0)
       {
         rel_mode = REL_MODE_SR_FREQ;
         Log ("Using full special relativity only for Doppler shifts, etc., and not considering co-moving frame effects.");
         j = i;
       }
-      else if (strcmp (argv[i], "-classic") == 0)
+      else if (strcmp (argv[i], "-nonrel") == 0)
       {
         rel_mode = REL_MODE_LINEAR;
         Log ("Using only old approach with linear Doppler shifts, etc. and not considering co-moving frame effects.\n");
@@ -352,7 +352,7 @@ help ()
 \n\
 This program simulates radiative transfer in a (biconical) CV, YSO, quasar, TDE or (spherical) stellar wind \n\
 \n\
-Usage:  py [-h] [-r] [-t time_max] [-v n] [--dry-run] [-i] [--version] [--rseed] [-p n_steps] [-classic] [-srclassic] xxx  or simply py \n\
+Usage:  py [-h] [-r] [-t time_max] [-v n] [--dry-run] [-i] [--version] [--rseed] [-p n_steps] [-nonrel] [-sr_doppler_only] xxx  or simply py \n\
 \n\
 where xxx is the rootname or full name of a parameter file, e. g. test.pf \n\
 \n\
@@ -387,9 +387,9 @@ These are largely diagnostic or for special cases. These include\n\
                         Range is in powers of 10, the difference beween the number of photons in the first cycle \n\
                         compared to the last. If range is missing, range is assumed to be 1, in which case the  \n\
                         number of photons will in the first cycle will be one order of magniude less than in the last cycle \n\
- -classic               Use Python in its classic configuration, with linear Doppler shifts, etc., and where co-moving frame\n\
+ -nonrel                Use Python in its old non-relativistic configuration, with linear Doppler shifts, etc., and where co-moving frame\n\
                         effects are not taken into account.\n\
- -srclassic             Use Python with full special relativity for Doppler shifts, etc., but do not include any co-moving frame\n\
+ -sr_doppler_only       Use Python with full special relativity for Doppler shifts, etc., but do not include any co-moving frame\n\
                         effects.\n\
  -ignore_partial_cells  Ignore wind cells that are only partially filled by the wind (This is now the default)  \n\
  -include_partial_cells Include wind cells that are only partially filled by the wind   \n\

source/py_optd.c

@@ -345,7 +345,7 @@ print_help (void)
     "A utility program to analyse the optical depth in a Python model.\n\n"
     "usage: py_optical_depth [-h] [-d ndom] [-p tau_stop] [-cion nion] \n"
     "                        [-freq_min min] [-freq_max max] [-i i1 i2 ...]\n"
-    "                        [--classic] [--smax frac] [--no-es] [--version]\n"
+    "                        [--nonrel] [--smax frac] [--no-es] [--version]\n"
     "                        root\n\n"
     "This program can be used in multiple ways. By default, the integrated continuum\n"
     "optical depth along the defined observer lines of sight for the model are returned.\n"
@@ -360,8 +360,8 @@ print_help (void)
     "-freq_min min    The lower frequency boundary for optical depth spectra\n"
     "-freq_max max    The upper frequency boundary for optical depth spectra\n"
     "-i i1 i2 i3 ...  Calculate the optical depth for the given space seperated list of sight lines\n"
-    "--classic        Use linear frequency transforms, to be used when Python was run\n"
-    "                 in classic mode\n"
+    "--nonrel         Use linear frequency transforms, to be used when Python was run\n"
+    "                 in non-relativistic mode\n"
     "--smax frac      Set the maximum fraction a photon can move in terms of cell distances\n"
     "--no-es          Do not include opacity contributions from electron scattering\n"
     "--version        Print the version information and exit.\n";
@@ -426,7 +426,7 @@ get_arguments (int argc, char *argv[])
       print_help ();
       exit (EXIT_SUCCESS);
     }
-    else if (!strcmp (argv[i], "--classic"))    //NOTE: use linear frequency transforms
+    else if (!strcmp (argv[i], "--nonrel"))    //NOTE: use linear frequency transforms
     {
       printf ("Using linear approximations for Doppler shifts\n");
       rel_mode = REL_MODE_LINEAR;

source/setup_line_transfer.c

@@ -73,58 +73,58 @@ get_line_transfer_mode ()
 
   if (user_line_mode == LINE_MODE_ABSORB)
   {
-    Log ("Line_transfer mode:  Simple, pure absorption, no scattering\n");
+    Log ("Line_transfer mode:  Classic, pure absorption, no scattering\n");
     geo.line_mode = user_line_mode;
   }
   else if (user_line_mode == LINE_MODE_SCAT)
   {
-    Log ("Line_transfer mode:  Simple, pure scattering, no absoprtion\n");
+    Log ("Line_transfer mode:  Classic, pure scattering, no absoprtion\n");
     geo.line_mode = user_line_mode;
   }
   else if (user_line_mode == LINE_MODE_SINGLE_SCAT)
   {
-    Log ("Line_transfer mode:  Simple, single scattering, with absorption, but without escape probality treatment\n");
+    Log ("Line_transfer mode:  Classic, single scattering, with absorption, but without escape probality treatment\n");
     geo.line_mode = user_line_mode;
   }
   else if (user_line_mode == LINE_MODE_ESC_PROB)
   {
-    Log ("Line_transfer mode:  Simple, isotropic scattering, escape probabilities\n");
+    Log ("Line_transfer mode:  Classic, isotropic scattering, escape probabilities\n");
     geo.line_mode = user_line_mode;
   }
   else if (user_line_mode == 5)
   {
-    Log ("Line_transfer mode:  Simple, thermal trapping, Single scattering \n");
+    Log ("Line_transfer mode: Classic, thermal trapping, Single scattering \n");
     geo.scatter_mode = SCATTER_MODE_THERMAL;    // Thermal trapping model
     geo.line_mode = LINE_MODE_ESC_PROB;
     geo.rt_mode = RT_MODE_2LEVEL;       // Not macro atom (SS) 
   }
   else if (user_line_mode == 6)
   {
-    Log ("Line_transfer mode:  macro atoms, isotropic scattering  \n");
+    Log ("Line_transfer mode: Macro, isotropic scattering  \n");
     geo.scatter_mode = SCATTER_MODE_ISOTROPIC;
     geo.line_mode = LINE_MODE_ESC_PROB;
     geo.rt_mode = RT_MODE_MACRO;        // Identify macro atom treatment (SS)
     geo.macro_simple = FALSE;   // We don't want the all simple case (SS)
   }
   else if (user_line_mode == 7)
   {
-    Log ("Line_transfer mode:  macro atoms, anisotropic  scattering  \n");
+    Log ("Line_transfer mode: Macro, anisotropic  scattering  \n");
     geo.scatter_mode = SCATTER_MODE_THERMAL;    // thermal trapping
     geo.line_mode = LINE_MODE_ESC_PROB;
     geo.rt_mode = RT_MODE_MACRO;        // Identify macro atom treatment (SS)
     geo.macro_simple = FALSE;   // We don't want the all simple case (SS)
   }
   else if (user_line_mode == 8)
   {
-    Log ("Line_transfer mode: simple macro atoms, isotropic  scattering  \n");
+    Log ("Line_transfer mode: Macro, isotropic  scattering  \n");
     geo.scatter_mode = SCATTER_MODE_ISOTROPIC;
     geo.line_mode = LINE_MODE_ESC_PROB;
     geo.rt_mode = RT_MODE_MACRO;        // Identify macro atom treatment i.e. indivisible packets
     geo.macro_simple = TRUE;    // This is for test runs with all simple ions (SS)
   }
   else if (user_line_mode == 9)
   {
-    Log ("Line_transfer mode: simple macro atoms, anisotropic  scattering  \n");
+    Log ("Line_transfer mode: Macro, anisotropic  scattering  \n");
     geo.scatter_mode = SCATTER_MODE_THERMAL;    // thermal trapping
     geo.line_mode = LINE_MODE_ESC_PROB;
     geo.rt_mode = RT_MODE_MACRO;        // Identify macro atom treatment i.e. indivisible packets
@@ -225,8 +225,8 @@ Available line transfer modes and descriptions are: \n\
   8(deprecated) Indivisible energy packets, force all simple-atoms, anisotropic scattering\n\
   9(deprecated) Indivisible energy packets, force all simple-atoms, anisotropic scattering\n\
 \n\
-  Standard mode is thermal_trapping for runs involving weight reduction and no macro-atoms\n\
-  Standard macro-atom mode is macro_atoms_thermal_trapping\n\
+  Classic mode is thermal_trapping for runs involving weight reduction and no macro-atoms\n\
+  Hybrid macro-atom mode is macro_atoms_thermal_trapping\n\
 \n\
 See this web address for more information: https://github.com/agnwinds/python/wiki/Line-Transfer-and-Scattering\n\
 \n\

source/templates.h

@@ -25,7 +25,8 @@ double upsilon(int n_coll, double u0);
 void skiplines(FILE *fptr, int nskip);
 /* bands.c */
 int bands_init(int imode, struct xbands *band);
-int freqs_init(double freqmin, double freqmax);
+int ion_bands_init(int mode, double freqmin, double freqmax, struct xbands *band);
+void check_appropriate_banding(struct xbands *band, int mode);
 /* bb.c */
 double planck(double t, double freqmin, double freqmax);
 double get_rand_pow(double x1, double x2, double alpha);