Cloudy & Associates

Commit 916e8e7a authored by Chatzikos, Marios's avatar Chatzikos, Marios
Browse files

Merge branch 'abs_feat' into 'master'

Report absorption features on output of 'save fine continuum' command

See merge request !37
parents ebd1ad4e 1de9b393
......@@ -448,6 +448,7 @@ appear, then only the intrinsic or the emergent, respectively,
block of lines will be inhibited.
\subsection{Print line optical depths [\_off, faint]}
\label{sec:PrintLineOptDep}
Mean line optical depths are not printed by default\footnote{Line center optical
depths were printed through version C10. Mean line optical depths
......@@ -1438,18 +1439,27 @@ The fine continuum is designed to account for line transfer.
It shows a normalized attenuated incident continuum and
does not include continuous
emission or absorption from the cloud.
This multi-grid approach is needed to combine
precision and speed. This command reports the line transmission coefficient,
$I_{transmitted}/I_{incident}$, for the fine continuum.
This multi-grid approach is needed to combine precision and speed.
This command reports the energy (see below) of the fine continuum bin in the
first column, and the continuum transmission coefficient,
$I_{transmitted}/I_{incident}$, for the fine continuum in the second column.
What follows is a list of absorption lines centered on that energy bin.
Each line is identified by its spectral label (species and wavelength), and is
followed by its mean optical depth.
The latter is the same as the last entry reported by
\cdCommand{print line optical depths}, see Section \ref{sec:PrintLineOptDep}.
All lines centered on that bin are reported, sorted from most-to-least opaque.
The resulting output will be huge if the entire fine continuum is saved.
The command accepts a \cdCommand{range} option to limit the size of the output file.
{\it Users are advised to make use of the {\rm\cdCommand{range}} option
to limit the size of the output file.}
If the keyword \cdCommand{range} appears then
the lower and upper limits to the range
of the fine continuum must be entered.
The command also accepts
a \cdCommand{units} option,
described on \pageref{output_units}, to change the units used
described on page \pageref{output_units}, to change the units used
in the specified energy range, and in resulting output.
The optional third numerical parameter gives the number
......@@ -1458,6 +1468,8 @@ reducing the size of the output file. The default is to average over 10
cells. If the number of cells to be combined is specified then it must
be the third number on the command line, following the limits
of the range.
Note that in a given energy range, the same number of features is reported,
irrespective of the number of cells merged.
The resolving power of the fine continuum is adjusted when the calculation begins so
that lines of the heavy elements can be resolved.
......
......@@ -65,6 +65,8 @@
#include <algorithm>
#include <fstream>
#include <bitset>
#include <unordered_map>
#include <numeric>
// Workaround for Windows...
#if defined(_MSC_VER) && !defined(SYS_CONFIG)
......
......@@ -69,6 +69,70 @@ void LineStackCreate()
if( trace.lgTrace )
fprintf( ioQQQ, "%7ld lines printed in main line array\n",
LineSave.nsum );
if( false )
{
//
// Go over the map of fine continuum indices just generated
// and quantify the amount of line overlap present.
// Currently (Jun 03, 2022) there exist bins with up to 14
// lines centered on them -- these are mainly due to hyperfine
// molecular lines.
//
// Create an ordered map of the number of lines per bin (key)
// and store a list (value) of all fine continuum indices that
// have that many lines overlap.
//
map<size_t, vector<long> > overlap;
for( auto it = rfield.fine_lstack.begin();
it != rfield.fine_lstack.end(); ++it )
{
if( overlap.find( it->second.size() ) == overlap.end() )
{
overlap.insert( make_pair( it->second.size(), vector<long>() ) );
}
overlap[ it->second.size() ].emplace_back( it->first );
}
fprintf( ioQQQ, "Overlap in fine continuum:\n" );
fprintf( ioQQQ, "==========================\n" );
const size_t nlines_per_bin = 10;
// loop over line overlap (number of lines in a fine cont bin)
//
for( auto it = overlap.begin(); it != overlap.end(); ++it )
{
if( it->first > nlines_per_bin )
{
sort( it->second.begin(), it->second.end() );
// loop over fine continuum indices
//
for( auto itf = it->second.begin();
itf != it->second.end(); ++itf )
{
TransitionProxy tr;
// loop over line stack indices
//
for( auto itl = rfield.fine_lstack[ *itf ].begin();
itl != rfield.fine_lstack[ *itf ].end(); ++itl )
{
tr = LineSave.lines[*itl].getTransition();
fprintf( ioQQQ, "%d\t%s\n",
*itl, tr.chLabel().c_str() );
}
fprintf( ioQQQ, "----------------------------\n" );
}
}
fprintf( ioQQQ, "%d lines per bin: %d bins\n",
int( it->first ), int( it->second.size() ) );
fprintf( ioQQQ, "~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n" );
}
}
}
/*eina convert a gf into an Einstein A */
......@@ -364,7 +428,7 @@ STATIC LinSv* lincom(
}
}
}
else if( LineSave.ipass == 0 )
{
LineSave.init(LineSave.nsum,(char) chInfo,chComment,chLab,lgAdd,wavelength,tr);
......@@ -381,8 +445,18 @@ STATIC LinSv* lincom(
fabs( rfield.anu(ipnt-1) - RYDLAM / wavelength) < error );
# endif
}
if( tr.associated() )
{
auto got = rfield.fine_lstack.find( tr.Emis().ipFine() );
if( got == rfield.fine_lstack.end() )
{
rfield.fine_lstack.insert( make_pair( tr.Emis().ipFine(), vector<int>() ) );
}
rfield.fine_lstack[ tr.Emis().ipFine() ].emplace_back( LineSave.nsum );
}
}
/* increment the line counter */
++LineSave.nsum;
......
......@@ -779,7 +779,7 @@ void ParseSave(Parser& p)
strcpy( save.chSave[save.nsave], "CONf" );
sncatf( chHeader,
"#Energy/%s\tTransmitted\n",
"#Energy/%s\tTransmitted\tSpecLine\tSingle-Line Opt Depth\n",
save.chConSavEnr[save.nsave] );
/* range option - important since so much data */
......
......@@ -508,6 +508,10 @@ public:
/** energies at center of each bin for fine continuum */
vector<realnum> fine_anu;
/** map of list of line stack indices (value)
* centered on fine opacity vector index (key) */
unordered_map<long, vector<int> > fine_lstack;
/** shift in fine continuum rest frame scale due to velocity gradient,
* rest frame is velocity of first zone, positive means that first zone
* is blue shifted relative to current zone. This is a decelerating
......
......@@ -305,7 +305,7 @@ STATIC void RT_line_fine_opacity(
/* line center opacity - type realnum since will add to fine opacity array,
* which is realnum */
realnum opac_line = (realnum)t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth;
realnum opac_line = (realnum)t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth;
// this is effective optical depth to this point. Do not do line if
// this product is less than SMALLFLOAT
......
......@@ -1166,7 +1166,7 @@ void SaveDo(
/* upper limit set with range option */
if( save.punarg[ipPun][1]> 0. )
nu_hi = ipFineCont( save.punarg[ipPun][1]);
nu_hi = ipFineCont( save.punarg[ipPun][1] );
else
nu_hi = rfield.nfine;
......@@ -1176,18 +1176,69 @@ void SaveDo(
do
{
realnum sum1 = rfield.fine_opt_depth[j];
realnum xnu = rfield.fine_anu[j];
for( long jj=1; jj<nskip; ++jj )
vector<int> all_stack_lines;
// when winds are present, line center in opacity
// vector is offset from its original position
// -- rfield.fine_lstack works with original position
//
long fine_index_no_wind = j - rfield.ipFineConVelShift;
realnum sum1 = 0.,
xnu = 0.;
for( long jj=0; jj<nskip; ++jj, ++j, ++fine_index_no_wind )
{
xnu += rfield.fine_anu[j+jj];
sum1 += rfield.fine_opt_depth[j+jj];
xnu += rfield.fine_anu[j];
sum1 += rfield.fine_opt_depth[j];
if( fine_index_no_wind > 0 )
{
auto got = rfield.fine_lstack.find( fine_index_no_wind );
if( got != rfield.fine_lstack.end() )
{
for( auto &lst_ind : got->second )
{
all_stack_lines.emplace_back( lst_ind );
}
}
}
}
// sort lines in decreasing opacity
//
sort( all_stack_lines.begin(),
all_stack_lines.end(),
[](int i1, int i2)
{
return LineSave.lines[i1].getTransition().Emis().TauInSpecific()
> LineSave.lines[i2].getTransition().Emis().TauInSpecific();
} );
double transm = sexp(sum1/nskip);
fprintf( save.params[ipPun].ipPnunit,
"%.6e\t%.3e\n",
AnuUnit(xnu/nskip),
sexp(sum1/nskip) );
j += nskip;
"%.6e\t%.3e",
AnuUnit(xnu/nskip),
transm );
static const realnum odep_limit = 0.01;
for( auto &ind: all_stack_lines )
{
TransitionProxy tr = LineSave.lines[ind].getTransition();
// Report mean optical depths, rather than line center
// Keeps optical depths consistent with main output
//
realnum odep = tr.Emis().TauInSpecific() * SQRTPI;
if( odep < odep_limit )
break;
fprintf( save.params[ipPun].ipPnunit,
"\t\"%s\"\t%.3e",
tr.chLabel().c_str(), odep );
}
fprintf( save.params[ipPun].ipPnunit, "\n" );
} while( j < nu_hi );
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment