Cloudy & Associates

Commit 59441169 authored by Chatzikos, Marios's avatar Chatzikos, Marios
Browse files

Merge branch 'hfs' into 'master'

Make small improvements re HFS lines

See merge request !41
parents 633f05c0 9770ef0d
......@@ -50,10 +50,10 @@ SILICON 14 30 3.0872
#
PHOSPHORUS 15 31 100.0
#
SULFUR 16 32 94.93
SULFUR 16 33 0.76
SULFUR 16 34 4.29
SULFUR 16 36 0.02
SULPHUR 16 32 94.93
SULPHUR 16 33 0.76
SULPHUR 16 34 4.29
SULPHUR 16 36 0.02
#
CHLORINE 17 35 75.78
CHLORINE 17 37 24.22
......
......@@ -50,10 +50,10 @@ SILICON 14 30 3.08716
#
PHOSPHORUS 15 31 100.0
#
SULFUR 16 32 95.018
SULFUR 16 33 0.75
SULFUR 16 34 4.215
SULFUR 16 36 0.017
SULPHUR 16 32 95.018
SULPHUR 16 33 0.75
SULPHUR 16 34 4.215
SULPHUR 16 36 0.017
#
CHLORINE 17 35 75.771
CHLORINE 17 37 24.229
......
......@@ -51,10 +51,10 @@ SILICON 14 30 3.087
#
PHOSPHORUS 15 31 100.0
#
SULFUR 16 32 95.018
SULFUR 16 33 0.75
SULFUR 16 34 4.215
SULFUR 16 36 0.017
SULPHUR 16 32 95.018
SULPHUR 16 33 0.75
SULPHUR 16 34 4.215
SULPHUR 16 36 0.017
#
CHLORINE 17 35 75.771
CHLORINE 17 37 24.229
......
......@@ -49,10 +49,10 @@ SILICON 14 30 3.0872
#
PHOSPHORUS 15 31 100.0
#
SULFUR 16 32 94.93
SULFUR 16 33 0.76
SULFUR 16 34 4.29
SULFUR 16 36 0.02
SULPHUR 16 32 94.93
SULPHUR 16 33 0.76
SULPHUR 16 34 4.29
SULPHUR 16 36 0.02
#
CHLORINE 17 35 75.78
CHLORINE 17 37 24.22
......
......@@ -55,10 +55,10 @@ SILICON 14 30 3.0872
#
PHOSPHORUS 15 31 100.0
#
SULFUR 16 32 94.93
SULFUR 16 33 0.76
SULFUR 16 34 4.29
SULFUR 16 36 0.02
SULPHUR 16 32 94.93
SULPHUR 16 33 0.76
SULPHUR 16 34 4.29
SULPHUR 16 36 0.02
#
CHLORINE 17 35 75.78
CHLORINE 17 37 24.22
......
# Isotope data
#
# List of astronomically relevant isotopes (see abundances/default-iso.abn).
# Spins and magnetic moments are listed mainly for isotopes discussed by
# Chatzikos et al. 2014, ApJ, 787, 96 (originally from Sunyaev & Churazov 1984,
# SvAL, 10, 201).
#
################################################################################
#
# Created: Aug 16, 2013
# @author: M. Chatzikos
#
# Updated: Jan 16, 2022
# @author: M. Chatzikos
# Comment: Added spins and magnetic moments for 3 stable isotopes: 61-Ni, 63-Cu,
# 65-Cu, and 67-Zn.
#
################################################################################
# Columns:
# (1) Element Symbol
# (2) Atomic Number
......@@ -80,13 +96,13 @@ Fe 26 58 -62148.800 57.9333
Co 27 59 -62223.600 58.9332 3.5 4.627000
Ni 28 58 -60223.000 57.9353
Ni 28 60 -64468.100 59.9308
Ni 28 61 -64216.800 60.9311
Ni 28 61 -64216.800 60.9311 1.5 -0.75002
Ni 28 62 -66742.700 61.9283
Ni 28 64 -67095.900 63.9280
Cu 29 63 -65576.200 62.9296
Cu 29 65 -67259.700 64.9278
Cu 29 63 -65576.200 62.9296 1.5 2.227206
Cu 29 65 -67259.700 64.9278 1.5 2.3816
Zn 30 64 -65999.500 63.9291
Zn 30 66 -68896.300 65.9260
Zn 30 67 -67877.200 66.9271
Zn 30 67 -67877.200 66.9271 2.5 0.875479
Zn 30 68 -70004.000 67.9248
Zn 30 70 -69559.000 69.9253
......@@ -2255,6 +2255,7 @@ This gives the depth [cm], the densities of \hO\ and \hplus\ [cm$^{-3}$],
followed by the level population densities [cm$^{-3}$]
for the levels in the order $n = 1, 2s, 2p, 3, 4, 5,\dots$ for each zone.
%------------------------------------------------------------------------------
\section{Save ionization means}
\label{sec:CommandSaveIonizationMeans}
......@@ -2741,32 +2742,66 @@ end of lines
This produces a list of the most important contributors to the line
radiation pressure for each zone.
\subsection{Save line populations, limit$=-$2}
This will output some information concerning the atomic parameters and
level populations for all lines that are transferred.
%------------------------------------------------------------------------------
\subsection{Save line populations
[row [lower$\vert$upper$\vert$ratio$\vert$tspin]]}
\label{sec:CommandSaveHyperfinePops}
This saves the upper and lower level populations of the requested spectral
lines.\footnote{
Until 2022, the command reported the level populations of (almost) all lines
known to the code.
This caused the output to be several MB even for a single zone model.
The current output was chosen as more user-friendly.
}
The command expects two filenames to be given, each enclosed in double quotes.
The first specifies the output file, while the second contains the spectral
lines to be tracked.
For example,
%
\begin{verbatim}
save line populations "output.txt" "linelist.dat"
\end{verbatim}
%
will read lines from \cdFilename{linelist.dat} and save the output in
\cdFilename{output.txt}.
By default, there is one line of output per spectral line per zone.
Each output line lists the depth, the transition line label,
the lower and upper level densities, the ratio $(n_u/g_u)/(n_l/g_l)$,
and the corresponding excitation temperature.
An example drawn from work on hyperfine lines in galaxy clusters is the
following:
%
\begin{verbatim}
#depth emline nl nu (nu/gu)/(nl/gl) Tspin
5.00000e-01 Fe26 348.000m 9.3587e-11 9.0239e-15 3.2141e-05 3.9964e+00
5.00000e-01 Fe24 3068.00m 1.2614e-08 1.0044e-08 2.6540e-01 3.5353e+00
\end{verbatim}
By default all transitions with upper level densities greater than zero
will be included. The lower limit to the density threshold can be reset
with the optional number than can appear on the line---this is the log of
the smallest population density [cm$^{-3}$] to be printed. This can be used
to make the print out somewhat shorter.
If the keyword \cdCommand{off} appears then
the limit to the smallest population to print will be turned off.
All
populations will be reported.
If the \cdCommand{row} keyword appears, the command will generate one line of
output for each zone, listing the depth, and for each spectral line, one of the
columns above.
By default, the population ratio is reported, but the lower level, upper level
populations, or the spin temperature may be reported, instead; the relevant
keywords are \cdCommand{ratio}, \cdCommand{lower}, \cdCommand{upper}, and
\cdCommand{tspin}, respectively.
For example, using the \cdCommand{row} option with the previous command,
modifies the output to
%
\begin{verbatim}
#depth Fe26 348.000m Fe24 3068.00m
5.00000e-01 3.2141e-05 2.6540e-01
\end{verbatim}
The first block of information gives an index to identify each line.
This is followed by the line's label. This will usually be the spectrum,
as in ``H 1'' followed by the wavelength, as in ``1215A''. The lower and
upper statistical weights are next, followed by the energy of the line in
wavenumbers and the $gf$ value for the transition.
These results may be used in post-processing.
For instance, the population ratio, $(n_u/g_u)/(n_l/g_l)$, may be used by
external programs to correct absorption coefficients for stimulated emission.
The population densities for each zone follow this block of information.
Each line begins with the index used for that line in the atomic parameter
list. This is followed by the populations of the lower and upper level
of the transition [cm$^{-3}$].
Note that results for lines that are not associated with an atomic model (e.g.,
blends, recombinations lines) are always given as zeros.
%------------------------------------------------------------------------------
\subsection{Save line RT}
This produces a file containing information concerning line radiative
......
......@@ -253,17 +253,8 @@ void H21_cm_pops( void )
* NB: continuum subtraction is performed within PutLine() */
set_xIntensity(HFLines[0]);
/* finally save the spin temperature */
hyperfine.Tspin21cm = phycon.te;
if( (*HFLines[0].Hi()).Pop() > SMALLFLOAT )
{
hyperfine.Tspin21cm = TexcLine( HFLines[0] );
/* this line must be non-zero - it does strongly mase in limit_compton_hi_t sim -
* in that sim pop ratio goes to unity for a float and TexcLine ret zero */
if( hyperfine.Tspin21cm == 0. )
hyperfine.Tspin21cm = phycon.te;
}
hyperfine.Tspin21cm = HyperfineTspin( HFLines[0] );
return;
}
......@@ -555,6 +546,7 @@ void HyperfineCreate(void)
colstr[j].cs.resize(Ntemp);
colstr[j].cs2d.resize(Ntemp);
}
hyperfine.HFLabundance.resize(HFLines.size());
/* now rewind the file so we can read it a second time*/
......@@ -632,7 +624,6 @@ void HyperfineCreate(void)
ASSERT(hyperfine.HFLabundance[j] >= 0.0 && hyperfine.HFLabundance[j] <= 1.0);
HFLines[j].Emis().Aul() = (realnum) atof(data[4].c_str());
HFLines[j].Emis().damp() = 1e-20f;
(*HFLines[j].Hi()).g() = (realnum) (2*(abund.IsoAbn[nelem-1].getSpin( Aiso ) + .5) + 1);
(*HFLines[j].Lo()).g() = (realnum) (2*(abund.IsoAbn[nelem-1].getSpin( Aiso ) - .5) + 1);
......@@ -649,6 +640,10 @@ void HyperfineCreate(void)
HFLines[j].WLAng() = (realnum)(wavelength * 1e8f);
HFLines[j].EnergyWN() = (realnum) fenergyWN;
HFLines[j].Emis().dampXvel() = (realnum)( HFLines[j].Emis().Aul()
/ HFLines[j].EnergyWN() / PI4 );
HFLines[j].Emis().damp() = 1e-20f;
HFLines[j].Emis().gf() = (realnum)(GetGF(HFLines[j].Emis().Aul(), fenergyWN, (*HFLines[j].Hi()).g()));
ASSERT(HFLines[j].Emis().gf() > 0.0);
......@@ -757,3 +752,27 @@ double HyperfineCS( size_t i )
return upsilon;
}
double HyperfineTspin( const TransitionProxy &t )
{
DEBUG_ENTRY( "HyperfineTspin()" );
double Tspin = phycon.te;
if( (*t.Hi()).Pop() > SMALLFLOAT )
{
Tspin = TexcLine( t );
/* this line MUST be non-zero
*
* H I 21cm does strongly mase in limit_compton_hi_t;
* in that sim pop ratio goes to unity for a float and
* TexcLine() returns zero */
if( Tspin == 0. )
Tspin = phycon.te;
}
ASSERT( Tspin > 0. or Tspin < 0. );
return Tspin;
}
......@@ -1338,7 +1338,7 @@ STATIC void conorm()
{
/* specify flux density
* option to use arbitrary frequency or range */
f = ffun1(rfield.range[i][0]);
f = ffun1(rfield.range[i][0]);
/* make sure this is positive, could be zero if out of range of table,
* or for various forms of insanity */
......
......@@ -1433,7 +1433,7 @@ STATIC void CoolHyperfine (void)
double cs = (electron_rate_21cm * dense.eden +
atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] +
proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) *
3./ dense.cdsqte;
3. / dense.cdsqte;
PutCS( cs , HFLines[0] );
/* do level pops for H 21 cm which is a special case since Lya pumping in included */
......
......@@ -16,6 +16,21 @@ void HyperfineCreate(void);
/*double HyperfineCS( long nelem , long ion );*/
double HyperfineCS( size_t i );
class TransitionProxy;
/** HyperfineTspin
*
* Compute spin (excitation) temperature for HFS line.
* The spin temperature is always finite (> 0, or < 0 if masing).
* The kinetic temperature is used if the spin temperature evaluates to 0,
* or if the upper limit has a population that is too small (~< 2e-36 cm^-3)
* to reliably compute the spin temperature.
*
* \param t [in] Transition
* \returns Tspin
*/
double HyperfineTspin( const TransitionProxy &t );
/** H21_cm_pops - fine level populations for 21 cm with Lya pumping included*/
void H21_cm_pops( void );
......
......@@ -8,7 +8,7 @@
/*LoadIsotopes read in the nuclear isotope data and allocate space */
void LoadIsotopes ( )
{
DEBUG_ENTRY( "SetIsotopeFractions()" );
DEBUG_ENTRY( "LoadIsotopes()" );
string chFile = "isotopes.dat";
FILE* ioDATA = open_data( chFile, "r" ); // will abort if not found
......
......@@ -399,6 +399,7 @@ public:
{
return m_type == QH;
}
#ifndef NDEBUG
void checkEmergent( const long ipEmType ) const
{
......
......@@ -981,21 +981,6 @@ void set_xIntensity( const TransitionProxy& t )
t.Emis().xObsIntensity() =
t.Emis().xIntensity() = nphot * t.EnergyErg();
if( 0 && t.chLabel() == "H 1 1215.67A" )
{
fprintf(ioQQQ,
"\"%s\"\t %.4e\t %.4e\t %.4e\t %.4e\t %.4e\t %.4e\t %.4e\t %.4e\n",
t.chLabel().c_str(),
t.Emis().Aul(),
(*t.Hi()).Pop(),
(*t.Lo()).Pop(),
rfield.flux_isotropic[j] * rfield.convoc[j],
Pesc,
nphot,
t.Emis().xIntensity(),
t.Emis().xObsIntensity() );
}
if( ! save.lgSubtrCont )
return;
......@@ -1011,5 +996,21 @@ void set_xIntensity( const TransitionProxy& t )
t.Emis().xObsIntensity() = nphot * t.EnergyErg();
if( 0 && t.chLabel() == "Fe24 3068.00m" )
{
fprintf(ioQQQ,
"\"%s\"\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
t.chLabel().c_str(),
t.Emis().Aul(),
(*t.Hi()).Pop(),
(*t.Lo()).Pop(),
rfield.flux_isotropic[j] * rfield.convoc[j],
Pesc,
nphot - dnphot, // undo the correction
dnphot,
t.Emis().xIntensity(),
t.Emis().xObsIntensity() );
}
return;
}
......@@ -1719,24 +1719,86 @@ void ParseSave(Parser& p)
else if( p.nMatch("POPU") )
{
/* save line populations command - first give index and inforamtion
* for all lines, then populations for lines as a function of
* depth, using this index */
/* save line populations "output" "LineList" */
strcpy( save.chSave[save.nsave], "LINP" );
sncatf( chHeader,
"#population information\n" );
/* this is optional limit to smallest population to save - always
* interpreted as a log */
save.punarg[save.nsave][0] = (realnum)exp10(p.FFmtRead());
/* this is default - all positive populations */
if( p.lgEOL() )
save.punarg[save.nsave][0] = 0.f;
/* we parsed off the second file name at start of this routine
* check if file was found, use it if it was, else abort */
if( !lgSecondFilename )
{
fprintf(ioQQQ ,
"There must be a second file name between double"
" quotes on the SAVE LINE POPULATIONS command\n."
"This second file contains the input list of"
" spectral lines.\n"
"I did not find it.\nSorry.\n");
cdEXIT(EXIT_FAILURE);
}
/* actually get the lines, and allocate the space in the arrays
* cdGetLineList will look on path, only do one time in grid */
if( save.params[save.nsave].ipPnunit == NULL )
{
/* make sure we free any allocated space from a previous call */
save.SaveLineListFree(save.nsave);
save.nLineList[save.nsave] = cdGetLineList(chSecondFilename, save.LineList[save.nsave]);
if( save.nLineList[save.nsave] < 0 )
{
fprintf(ioQQQ,
"DISASTER could not open"
" SAVE LINE POPULATIONS file %s \n",
chSecondFilename.c_str() );
cdEXIT(EXIT_FAILURE);
}
}
if( p.nMatch(" OFF") )
// by default give numbers in columns
// ROW keyword says to write the numbers across as one long row
// subsequent keyword tells which column to write
if( p.nMatch(" ROW") )
{
/* no lower limit - print all lines */
save.punarg[save.nsave][0] = -1.f;
save.punarg[save.nsave][0] = 1;
if( p.nMatch( " LOW" ) )
/* lower level population */
save.punarg[save.nsave][1] = 0;
else if( p.nMatch( " UPP" ) )
/* upper level population */
save.punarg[save.nsave][1] = 1;
else if( p.nMatch( " TSP" ) )
/* spin temperature */
save.punarg[save.nsave][1] = 3;
else
/* DEFAULT: population ratio: (nu/gu)/(nl/gl) */
save.punarg[save.nsave][1] = 2;
}
else
// the default, one line per row, multiple columns
save.punarg[save.nsave][0] = 0;
sncatf( chHeader, "#depth\t" );
if( save.punarg[save.nsave][0] )
{
for( long int j=0; j<save.nLineList[save.nsave]; ++j )
{
sncatf( chHeader, "%s ",
save.LineList[save.nsave][j].chLabel.c_str() );
string chTemp;
sprt_wl( chTemp, save.LineList[save.nsave][j].wave );
sncatf( chHeader, "%s", chTemp.c_str() );
if( j != save.nLineList[save.nsave] )
{
sncatf( chHeader, "\t" );
}
}
sncatf( chHeader, "\n" );
}
else
{
sncatf( chHeader,
"emline\tnl\tnu\t(nu/gu)/(nl/gl)\tTspin\n" );
}
}
......
......@@ -2514,15 +2514,86 @@ void SaveDo(
{
if( ! lgLastOnly )
{
static bool lgFirst=true;
/* save line populations, need to do this twice if very first
* call since first call to SaveLineStuff generates atomic parameters
* rather than level pops, routine is below, file static */
SaveLineStuff(save.params[ipPun].ipPnunit,"populat" , save.punarg[ipPun][0]);
if( lgFirst )
if( save.punarg[ipPun][0] )
{
lgFirst = false;
SaveLineStuff(save.params[ipPun].ipPnunit,"populat" , save.punarg[ipPun][0]);
fprintf( save.params[ipPun].ipPnunit,
"%.5e", radius.depth_mid_zone );
for( auto &specline : save.LineList[ipPun] )
{
long ipobs = LineSave.findline(specline);
TransitionProxy tr = LineSave.lines[ipobs].getTransition();
double quant = 0.;
/* NB NB
*
* Recombination lines are not associated with an atomic
* model. Weed them out to avoid segfault.
*/
if( LineSave.lines[ipobs].chSumTyp() == 't'
&& tr.associated() )
{
if( save.punarg[ipPun][1] == 0 )
quant = tr.Lo()->Pop();
else if( save.punarg[ipPun][1] == 1 )
quant = tr.Hi()->Pop();
else if( save.punarg[ipPun][1] == 2 )
quant = safe_div( tr.Hi()->Pop() * tr.Lo()->g(),
tr.Lo()->Pop() * tr.Hi()->g(), 0. );
else if( save.punarg[ipPun][1] == 3 )
quant = TexcLine( tr );
else
{
TotalInsanity();
}
}
fprintf( save.params[ipPun].ipPnunit, "\t%.4e", quant );
}
fprintf( save.params[ipPun].ipPnunit, "\n" );
}
else
{
for( auto &specline : save.LineList[ipPun] )
{
fprintf( save.params[ipPun].ipPnunit,
"%.5e", radius.depth_mid_zone );
fprintf( save.params[ipPun].ipPnunit,
"\t%s ",
specline.chLabel.c_str() );
string chTemp;
sprt_wl( chTemp, specline.wave );
fprintf( save.params[ipPun].ipPnunit,
"%s", chTemp.c_str() );
double lower = 0.,
upper = 0.,
ratio = 0.,
tspin = 0.;
long ipobs = LineSave.findline(specline);
TransitionProxy tr = LineSave.lines[ipobs].getTransition();
/* NB NB
*
* Recombination lines are not associated with an atomic
* model. Weed them out to avoid segfault.
*/
if( LineSave.lines[ipobs].chSumTyp() == 't'
&& tr.associated() )
{
lower = tr.Lo()->Pop();
upper = tr.Hi()->Pop();
ratio = safe_div( tr.Hi()->Pop() * tr.Lo()->g(),
tr.Lo()->Pop() * tr.Hi()->g(), 0. );
tspin = TexcLine( tr );
}
fprintf( save.params[ipPun].ipPnunit,
"\t%10.4e\t%10.4e\t%.4e\t%11.4e\n",
lower, upper, ratio, tspin );
}
}
}
}
......
......@@ -40,6 +40,14 @@ save continuum "limit_caseb_he_den2_temp4.con" last no title units microns
save monitors "limit_caseb_he_den2_temp4.asr"
save line list ratio column "limit_caseb_he_den2_temp4.rat" "linelist_he1.dat" last no hash
#
save line populations "limit_caseb_he_den2_temp4.lpop" "linelist_he1.dat" last no hash
# default row option: ratio
save line populations "limit_caseb_he_den2_temp4.row.lpop" "linelist_he1.dat" row last no hash
save line populations "limit_caseb_he_den2_temp4.row_lower.lpop" "linelist_he1.dat" row lower last no hash
save line populations "limit_caseb_he_den2_temp4.row_upper.lpop" "linelist_he1.dat" row upper last no hash
save line populations "limit_caseb_he_den2_temp4.row_ratio.lpop" "linelist_he1.dat" row ratio last no hash
save line populations "limit_caseb_he_den2_temp4.row_tspin.lpop" "linelist_he1.dat" row tspin last no hash
#
# commands giving the monitors =========
##
## 15 may 22 add this monitor, which botches 10830
......
Supports Markdown
0% or .
You are about to add