Merge branch 'master' into coll_bug

docs/latex/QuickStart/VeuszCookbook.tex

@@ -229,6 +229,7 @@ You can also place graphs within a grid widget to get an arrangement of separate
 Figure \ref{fig:orion_hii_pdr_pp} shows an example, 
 taken from the Orion HII region / PDR simulation, of a nested plot.
 
+
 \begin{figure}
 \begin{center}
 \includegraphics[clip=on,width=0.8\columnwidth,height=0.8\textheight,keepaspectratio]{orion_hii_pdr_pp}
@@ -239,4 +240,13 @@ Orion H II region / PDR.}
 \label{fig:orion_hii_pdr_pp}
 \end{figure}
 
+
+\subsection{Documenting your work}
+
+The Properties widget includes a Notes field that can be used to keep a history of the 
+data behind a plot.
+This is useful so that you can come back to a project after some time and recover
+the history of the work leading up to a result or the location of the helper files that were used
+along the way.
+
 %\end{document}

docs/latex/frontis_common.tex

 {\small
 \noindent
-{\em Software:} Copyright \copyright\ 1978-2019 Gary J. Ferland and others. All rights reserved.
+{\em Software:} Copyright \copyright\ 1978-2022 Gary J. Ferland and others. All rights reserved.
 
 \vspace{1em}
 
@@ -44,7 +44,7 @@ to the following FreeBSD format documentation license:
 \vspace{1em}
 
 \noindent
-{\em Documentation:} Copyright \copyright\ 1978--2019 Gary J. Ferland and others. All rights reserved.
+{\em Documentation:} Copyright \copyright\ 1978--2022 Gary J. Ferland and others. All rights reserved.
 
 \vspace{1em}
 

source/cloudy.cpp

@@ -127,6 +127,15 @@ bool cloudy()
 
 	SetPrintLineCol();
 
+        /* Initialize pseudo-continua, if requested */
+	SpeciesPseudoContCreate();
+
+	/* Initialize species bands, if requested 
+	 *
+	 * NB NB
+	 * This function requires that the atomic data be read in. */
+	SpeciesBandsCreate();
+
 	/* print header */
 	PrtHeader();
 

source/continuum.h

@@ -115,6 +115,15 @@ void SpeciesPseudoContCreate();
 /** SpeciesPseudoContAccum - accumulate pseudo-continua */
 void SpeciesPseudoContAccum();
 
+/** addUniqueSpeciesBand - add unique combination of species and bands file to list
+ *
+ *  Note that this command requires that the atomic data be read in.
+ *
+ *  \param [in] filename     name of bands file
+ *  \param [in[ speciesLabel species whose emission to accumulate
+ */
+void addUniqueSpeciesBand( const string &filename, const string &speciesLabel );
+
 /** SpeciesBandsCreate - initialize requested species bands files */
 void SpeciesBandsCreate();
 

source/cool_eval.cpp

@@ -1293,9 +1293,9 @@ STATIC double CoolH2_GA08 ( double Temp )
 {
 	DEBUG_ENTRY( "CoolH2_GA08()" );
 
-	double x_para = h2.para_density / dense.gas_phase[ ipHYDROGEN ];
-	double x_ortho = h2.ortho_density / dense.gas_phase[ ipHYDROGEN ];
-	double f_para = h2.para_density / (h2.para_density + h2.ortho_density);
+	double x_para = safe_div( h2.para_density, double( dense.gas_phase[ ipHYDROGEN ] ) );
+	double x_ortho = safe_div( h2.ortho_density, double( dense.gas_phase[ ipHYDROGEN ] ) );
+	double f_para = safe_div( h2.para_density, h2.para_density + h2.ortho_density );
 	double f_ortho = 1. - f_para;
 
 	/* Collisions with atomic H */
@@ -1354,14 +1354,10 @@ STATIC double CoolH2_GA08 ( double Temp )
 		}
 	}
 
-	/* NB NB NB
-	 * 
-	 * The interpolation below will fail if the high density cooling is 0,
-	 * that is, the total cooling will be 0, even if the low-density cooling
-	 * is not.  This may happen if permitted by the function that computes
-	 * LTE cooling, and/or the tabulated data (see definition of
-	 * interpolate_LTE_Cooling() and associated data file). */
-	return cooling_high / (1. + cooling_high / cooling_low);
+	// originally cooling_high / (1. + cooling_high / cooling_low)
+	// See Eq. (39) therein
+	//
+	return safe_div( cooling_low * cooling_high, cooling_low + cooling_high );
 }
 
 STATIC void CoolHyperfine (void)

source/generic_state.cpp

@@ -35,6 +35,7 @@ double density(const genericState& gs)
 	else
 		return 0.;
 }
+
 double depart(const genericState& gs)
 {
 	if ( gs.sp != NULL )
@@ -46,6 +47,7 @@ double depart(const genericState& gs)
 	else
 		return 1.0;
 }
+
 double energy(const genericState& gs) 
 {
 	if ( gs.sp != NULL )
@@ -57,6 +59,7 @@ double energy(const genericState& gs)
 	else
 		return 0.0;
 }
+
 double levels(const genericState& gs) 
 {
 	if ( gs.sp != NULL && gs.sp->levels != NULL )
@@ -67,9 +70,17 @@ double levels(const genericState& gs)
 		return 0.0;
 }
 
+bool genericState::isSpecies() const
+{
+	if( sp != NULL && sp != null_molezone )
+		return true;
+	else
+		return false;
+}
+
 string genericState::label() const
 {
-	if ( sp != NULL && sp != null_molezone )
+	if( isSpecies() )
 		return sp->global()->label;
 	else if (q.associated())
 		return q.chLabel();
@@ -81,7 +92,7 @@ string genericState::label() const
 
 string genericState::database() const
 {
-	if ( sp != NULL && sp != null_molezone && sp->dbase != NULL )
+	if( isSpecies() && sp->dbase != NULL )
 		return sp->dbase->database;
 	else
 		return "";
@@ -89,7 +100,7 @@ string genericState::database() const
 
 bool genericState::associated() const
 {
-	if (sp != NULL && sp != null_molezone)
+	if( isSpecies() )
 		return true;
 	else if (q.associated())
 		return true;
@@ -342,3 +353,18 @@ vector<genericState> matchGeneric(const string &chLabel, bool lgValidate)
 	}
 	return v;
 }
+
+genericState getSpeciesGeneric(const string &speciesLabel)
+{
+	DEBUG_ENTRY( "getSpeciesGeneric()" );
+
+	vector<genericState> v = matchGeneric( speciesLabel, true );
+	if( v.size() == 0 || ! v[0].isSpecies() )
+	{
+		fprintf( ioQQQ, "PROBLEM: Invalid species '%s' requested.\n",
+				speciesLabel.c_str() );
+		cdEXIT( EXIT_FAILURE );
+	}
+
+	return v[0];
+}

source/generic_state.h

@@ -23,6 +23,7 @@ public:
 	bool associated() const;
 	string label() const;
 	string database() const;
+	bool isSpecies() const;
 };
 
 double column(const genericState&);
@@ -58,4 +59,13 @@ const molezone *getLevelsGeneric(const string &chLabel, bool lgValidate, vector<
  */
 vector<genericState> matchGeneric(const string &chLabel, bool lgValidate);
 
+/** getSpeciesGeneric - acquire the species matching the input string
+ *
+ * \param [in] speciesLabel  input species string
+ *
+ * \returns species     reference to requested species
+ */
+genericState getSpeciesGeneric( const string &speciesLabel );
+
+
 #endif // GENERIC_STATE_H_

source/hextra.h

@@ -109,7 +109,7 @@ struct t_hextra : public module {
 		 hextra.background_density = 1.99e-9f;*/
 		/* >>chng 05 apr 16, to get proper ionization rate in ism_set_cr_rate, where
 		 * H is forced to be fully atomic, no molecules, density from 1.99 to 2.15 */
-		/* >>chng 02 apr 05, update to
+		/* >>chng 12 apr 05, update to
 		 * >>refer	cosmic ray 	ionization	Indriolo, N., Geballe, T., Oka, T., & McCall, B.J. 2007, ApJ, 671, 1736
 		 */
 		background_density = 2.15e-9f*7.9f;

source/init_sim_postparse.cpp

@@ -196,12 +196,6 @@ void InitSimPostparse( void )
 	MonitorResults.SumErrorCaseMonitor = 0.;
 	MonitorResults.nSumErrorCaseMonitor = 0;
 
-	/* Initialize pseudo-continua, if requested */
-	SpeciesPseudoContCreate();
-
-	/* Initialize species bands, if requested */
-	SpeciesBandsCreate();
-
 	mean.setup_molecules();
 
 	return;

source/parse_save.cpp

@@ -22,12 +22,11 @@
 #include "service.h"
 #include "species.h"
 #include "dense.h"
+#include "continuum.h"
 
 /* check for keyword UNITS on line, then scan wavelength or energy units if present */
 STATIC const char* ChkUnits(Parser &p);
 
-STATIC void addUniqueSpeciesBand( const string &filename, const string &speciesLabel );
-
 inline void saveXSPEC(unsigned int option)
 {
 	static const char* labelXSPEC[NUM_OUTPUT_TYPES] = {
@@ -2615,13 +2614,6 @@ void ParseSave(Parser& p)
 		ioMAP = save.params[save.nsave].ipPnunit;
 	}
 
-	/* make sure FeII bands are always processed
-	 * if a 'save species bands' command has not been issued
-	 * the bands will be computed, and printed on main output,
-	 * but no 'save' output file will be created */
-	if( dense.lgElmtOn[ipIRON] )
-		addUniqueSpeciesBand( "FeII_bands.ini", "Fe+" );
-
 	/* if not done already and chTitle has been set to a string then print title
 	 * logic to prevent more than one title in grid calculation */
 	if( save.lgSaveTitle(save.nsave) && chTitle.length() > 0 )
@@ -2823,38 +2815,3 @@ STATIC const char* ChkUnits( Parser &p )
 	}
 	return val;
 }
-
-STATIC bool specBandsExists( const string &filename, const string &speciesLabel )
-{
-	DEBUG_ENTRY( "specBandsExists()" );
-
-	bool exists = false;
-
-	for( vector<save_species_bands>::iterator it = save.specBands.begin();
-		it != save.specBands.end(); ++it )
-	{
-		if( (*it).filename == filename &&
-			(*it).speciesLabel == speciesLabel )
-		{
-			exists = true;
-			break;
-		}
-	}
-
-	return exists;
-}
-
-STATIC void addUniqueSpeciesBand( const string &filename, const string &speciesLabel )
-{
-	DEBUG_ENTRY( "addUniqueSpeciesBand()" );
-
-	if( specBandsExists( filename, speciesLabel ) )
-		return;
-
-	save_species_bands thisSpBand;
-	thisSpBand.filename = filename;
-	thisSpBand.speciesLabel = speciesLabel;
-	save.specBands.push_back( thisSpBand );
-
-	return;
-}

source/save.h

@@ -228,13 +228,6 @@ public:
 	long nBins;
 };
 
-class save_species_bands
-{
-public:
-	string filename;
-	string speciesLabel;
-};
-
 class save_img_matrix : public t_prt_matrix
 {
 public:
@@ -636,7 +629,6 @@ public:
 
 	/** Parameters for species bands */
 	string SpeciesBandFile[LIMPUN];
-	vector<save_species_bands> specBands;
 
 	save_img_matrix img_matrix;
 };

source/save_species.cpp

@@ -451,8 +451,7 @@ void SaveSpeciesOptDep( const long int ipPun, const string &speciesLabel )
 		save.SaveHeaderDone( ipPun );
 	}
 
-	genericState species;
-	getSpecies( speciesLabel, species );
+	genericState species = getSpeciesGeneric( speciesLabel );
 
 	if( species.sp->lines == NULL )
 	{

source/species.cpp

@@ -916,6 +916,24 @@ void chemical_to_spectral( const string &chLabelChem, string &chLabelSpec )
 	}
 }
 
+/* Tell if the given species is enabled */
+bool isSpeciesActive( const string &chLabelChem )
+{
+	DEBUG_ENTRY( "isSpeciesActive()" );
+
+	for( size_t i=0; i<mole_global.list.size(); ++i )
+	{
+		if( mole.species[i].levels == NULL )
+			continue;
+		if( mole.species[i].dbase == NULL )
+			continue;
+		if( mole.species[i].levels->chLabel() == chLabelChem )
+			return mole.species[i].dbase->lgActive;
+	}
+
+	return false;
+}
+
 /*This function fills the nelem and IonStg fields */
 STATIC void states_nelemfill(void)
 {

source/species.h

@@ -135,14 +135,19 @@ bool parse_chemical( const string &chLabelChem,
  */
 void chemical_to_spectral( const string &chLabelChem, string &chLabelSpec );
 
-class genericState;
-
-/** getSpecies - acquire the species matching the input string
+/**
+ * isSpeciesActive - Tell if the species is enabled
+ *
+ * An atomic/ionic species may be inactive if the relevant element is disabled.
+ * Likewise, a molecular species may be inactive if one of its constituents
+ * elements is disabled.  This function reports on the given species.
+ * If not found, the result is false.
+ *
+ * \param chLabelChem [in]	Chemical label, e.g., "C+2", or "HCl"
  *
- * \param speciesLabel	input species string
- * \param species	output reference to requested species
+ * \return bool			True, if species active; false, otherwise.
  */
-void getSpecies( const string &speciesLabel, genericState &species );
+bool isSpeciesActive( const string &chSpecLabel );
 
 /** isAtomicIonValid - Tell if an atomic ion is has a charge consistent
  * 			with its atomic number.

source/species_pseudo_cont.cpp

@@ -45,22 +45,6 @@ STATIC string getIntenTypeStr( const int ipContType )
 	}
 }
 
-void getSpecies( const string &speciesLabel, genericState &species )
-{
-	DEBUG_ENTRY( "getSpecies()" );
-
-	vector<genericState> v = matchGeneric( speciesLabel, false );
-	if( v.size() != 1 )
-	{
-		fprintf( ioQQQ, "Error: Incorrect number of matches"
-				" (%d) for species '%s'\n",
-				int(v.size()), speciesLabel.c_str() );
-		cdEXIT( EXIT_FAILURE );
-	}
-	//	printf( "sp= '%s'\n", v[0].label().c_str() );
-	species = v[0];
-}
-
 
 
 /*==============================================================================*/
@@ -141,7 +125,7 @@ double band_cont::getInten( const long ibin, const int ipContType ) const
 			inten = inten_inward[ ibin ] + inten_outward[ ibin ];
 			break;
 		default:
-			fprintf( ioQQQ, "Error: Illegal continuum type: %d\n",
+			fprintf( ioQQQ, "PROBLEM: Illegal continuum type: %d\n",
 					ipContType );
 			cdEXIT( EXIT_FAILURE );
 			break;
@@ -170,7 +154,7 @@ private:
 	{
 		if( ! check_index( ibin ) )
 		{
-			fprintf( ioQQQ, "Error: Pseudo-continuum bin (%ld) "
+			fprintf( ioQQQ, "PROBLEM: Pseudo-continuum bin (%ld) "
 					"for species '%s' "
 					"out of range (0, %ld)\n",
 					ibin,
@@ -211,7 +195,7 @@ void pseudo_cont::setup( string &label, double wlo, double whi, long nb )
 		//			ibin, wl[ ibin ] );
 	}
 
-	getSpecies( speciesLabel, species );
+	species = getSpeciesGeneric( speciesLabel );
 	//	printf( "sp= '%s'\n", species.label().c_str() );
 }
 
@@ -274,7 +258,7 @@ void pseudo_cont::sumBand( double *sumOutward, double *sumInward ) const
 }
 /*============================================================================*/
 
-static vector< pseudo_cont > PseudoCont;
+static vector<pseudo_cont> PseudoCont;
 
 STATIC void getPseudoIndex( const string &speciesLabel,
 				vector<pseudo_cont>::iterator &this_it )
@@ -418,7 +402,7 @@ void SaveSpeciesPseudoCont( const long ipPun, const string &speciesLabel )
 	if( it == PseudoCont.end() )
 	{
 		fprintf( ioQQQ,
-			"Error: Species continuum data unmatched for species '%s'\n",
+			"PROBLEM: Species continuum data unmatched for species '%s'\n",
 			speciesLabel.c_str() );
 		cdEXIT( EXIT_FAILURE );
 	}
@@ -649,6 +633,8 @@ bool bands_file::load()
 }
 
 /*==============================================================================*/
+
+/** Bands -- list of bands files specified with 'save species bands' commands */
 static vector<bands_file> Bands;
 
 STATIC void findBandsFile( const string &filename,
@@ -689,19 +675,20 @@ STATIC void addBandsFile( const string &filename )
 /*==============================================================================*/
 /*				SPECIES BAND EMISSION				*/
 /*==============================================================================*/
-class species_bands : public band_cont
+class band_emission : public band_cont
 {
 private:
 	string bandLabel,
 		comment;
 	vector<bands_file>::iterator bands_it;
+	bool isInitd = false;
 public:
 	const string inwdLabel = "InwdBnd";
 
 	void setup( const string &splab, vector<bands_file>::iterator it )
 	{
 		speciesLabel = splab;
-		getSpecies( splab, species );
+		species = getSpeciesGeneric( splab );
 		//	printf("species: '%s'\n", species.label().c_str());
 
 		bands_it = it;
@@ -721,13 +708,14 @@ public:
 		bandLabel = spectralLabel + "b";
 		comment = spectralLabel + " emission in bands defined in " +
 				(*bands_it).bandFilename();
+		isInitd = true;
 	}
 private:
 	void check_index_fatal( const long iband ) const
 	{
 		if( ! check_index( iband ) )
 		{
-			fprintf( ioQQQ, "Error: Band (%ld) "
+			fprintf( ioQQQ, "PROBLEM: Band (%ld) "
 					"for species '%s' "
 					"from file '%s' "
 					"out of range (0, %ld)\n",
@@ -760,11 +748,15 @@ public:
 		check_index_fatal( iband );
 		return (*bands_it).getWlHi( iband );
 	}
+	bool initialized() const
+	{
+		return isInitd;
+	}
 };
 
-void species_bands::sumBand( double *sumOutward, double *sumInward ) const
+void band_emission::sumBand( double *sumOutward, double *sumInward ) const
 {
-	DEBUG_ENTRY( "species_bands::sumBand()" );
+	DEBUG_ENTRY( "band_emission::sumBand()" );
 
 	for( long i=0; i<nBins; ++i )
 	{
@@ -801,9 +793,9 @@ void species_bands::sumBand( double *sumOutward, double *sumInward ) const
 	}
 }
 
-void species_bands::insert()
+void band_emission::insert()
 {
-	DEBUG_ENTRY( "species_bands::insert()" );
+	DEBUG_ENTRY( "band_emission::insert()" );
 
 	for( long iband = 0; iband < nBins; iband++ )
 	{
@@ -823,25 +815,57 @@ void species_bands::insert()
 }
 /*============================================================================*/
 
-static vector<species_bands> SpecBands;
-
-STATIC void getSpecBandsIndex( const string &speciesLabel, const string &fileBands,
-				 vector<species_bands>::iterator &this_it )
+class species_band
 {
-	DEBUG_ENTRY( "getSpecBandsIndex()" );
+public:
+	string filename;
+	string speciesLabel;
+	band_emission bandEmission;
+
+	explicit species_band(const string &file, const string &species)
+		: filename(file), speciesLabel(species) {};
+};
+
+static vector<species_band> SpecBands;
+typedef vector<species_band>::iterator sb_itor;
 
-	this_it = SpecBands.end();
+STATIC sb_itor findSpecBand( const string &filename, const string &speciesLabel )
+{
+	DEBUG_ENTRY( "findSpecBand()" );
 
-	for( vector<species_bands>::iterator it = SpecBands.begin();
-		it != SpecBands.end(); ++it )
+	for( auto it = SpecBands.begin(); it != SpecBands.end(); ++it )
 	{
-		if( speciesLabel == (*it).label() &&
-			fileBands == (*it).bandFilename() )
+		if( it->filename == filename && it->speciesLabel == speciesLabel )
 		{
-			this_it = it;
-			break;
+			return it;
 		}
 	}
+
+	return SpecBands.end();
+}
+
+STATIC bool specBandsExists( const string &filename, const string &speciesLabel )
+{
+	DEBUG_ENTRY( "specBandsExists()" );
+
+	auto it = findSpecBand( filename, speciesLabel );
+	if( it == SpecBands.end() )
+		return false;
+	else
+		return true;
+}
+
+void addUniqueSpeciesBand( const string &filename, const string &speciesLabel )
+{
+	DEBUG_ENTRY( "addUniqueSpeciesBand()" );
+
+	if( specBandsExists( filename, speciesLabel ) )
+		return;
+
+	species_band thisSpBand( filename, speciesLabel );
+	SpecBands.push_back( thisSpBand );
+
+	return;
 }
 
 /*============================================================================*/
@@ -850,23 +874,31 @@ void SpeciesBandsCreate()
 {
 	DEBUG_ENTRY( "SpeciesBandsCreate()" );
 
-	// Already initialized
-	if( SpecBands.size() != 0 )
-		return;
+	/* make sure FeII bands are always processed
+	 *
+	 * if a 'save species bands' command has not been issued
+	 * the bands will be computed, and printed on main output,
+	 * but no 'save' output file will be created
+	 */
+	addUniqueSpeciesBand( "FeII_bands.ini", "Fe+" );
 
-	for( auto it = save.specBands.begin(); it != save.specBands.end(); ++it )
+	for( auto it = SpecBands.begin(); it != SpecBands.end(); ++it )
 	{
+		if( ! isSpeciesActive( (*it).speciesLabel ) )
+			continue;
+
 		addBandsFile( (*it).filename );
 	}
 
-	for( auto it = save.specBands.begin(); it != save.specBands.end(); ++it )