Cloudy & Associates

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

Identify most opaque line in 'save fine continuum'

Anna Ogorzalek posted a query on the Cloudy group about identifying
absorption features in the output of the 'save xspec mtable' command.

A partial answer would be to identify features due to the lines alone.

The present adds two columns to the output of 'save fine continuum'.
The first identifies the species whose line is centered at that bin,
and the second reports the optical depth due that line alone (i.e., not
the total including any neighboring lines).
parent 20030830
......@@ -65,6 +65,7 @@
#include <algorithm>
#include <fstream>
#include <bitset>
#include <unordered_map>
// Workaround for Windows...
#if defined(_MSC_VER) && !defined(SYS_CONFIG)
......
......@@ -69,6 +69,45 @@ void LineStackCreate()
if( trace.lgTrace )
fprintf( ioQQQ, "%7ld lines printed in main line array\n",
LineSave.nsum );
if( false )
{
fprintf( ioQQQ, "Overlap in fine continuum:\n" );
fprintf( ioQQQ, "==========================\n" );
int nlines_per_bin = 0;
while( nlines_per_bin < 20 )
{
int ntot = 0;
for( auto it = rfield.fine_lstack.begin();
it != rfield.fine_lstack.end(); ++it )
{
if( int( it->second.size() ) == nlines_per_bin )
{
ntot++;
if( nlines_per_bin > 10 )
{
TransitionProxy tr;
for( auto itl = it->second.begin();
itl != it->second.end(); ++itl )
{
tr = LineSave.lines[*itl].getTransition();
fprintf( ioQQQ, "%s\n", tr.chLabel().c_str() );
}
fprintf( ioQQQ, "----------------------------\n" );
}
}
}
if( ntot > 0 )
{
fprintf( ioQQQ, "%d lines per bin: %d bins\n",
nlines_per_bin, ntot );
fprintf( ioQQQ, "----------------------------\n" );
}
nlines_per_bin++;
}
}
}
/*eina convert a gf into an Einstein A */
......@@ -364,7 +403,7 @@ STATIC LinSv* lincom(
}
}
}
else if( LineSave.ipass == 0 )
{
LineSave.init(LineSave.nsum,(char) chInfo,chComment,chLab,lgAdd,wavelength,tr);
......@@ -381,8 +420,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\tSpecies\tOpt Dep Spec\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,17 +1176,72 @@ void SaveDo(
do
{
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;
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 );
}
}
}
realnum sum1 = rfield.fine_opt_depth[j];
realnum xnu = rfield.fine_anu[j];
for( long jj=1; jj<nskip; ++jj )
{
xnu += rfield.fine_anu[j+jj];
sum1 += rfield.fine_opt_depth[j+jj];
fine_index_no_wind = j + jj - rfield.ipFineConVelShift;
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 );
}
}
}
}
// identify most opaque line in bin
//
string label = "";
realnum max_odep = 0.;
for( auto &lst_ind : all_stack_lines )
{
realnum this_odep =
LineSave.lines[lst_ind].getTransition()
.Emis().TauInSpecific();
if( this_odep > max_odep )
{
max_odep = this_odep;
label = string( LineSave.lines[lst_ind].chALab() );
}
}
double transm = sexp(sum1/nskip);
fprintf( save.params[ipPun].ipPnunit,
"%.6e\t%.3e\n",
AnuUnit(xnu/nskip),
sexp(sum1/nskip) );
"%.6e\t%.3e",
AnuUnit(xnu/nskip),
transm );
if( label != "" and max_odep > 0.01 )
fprintf( save.params[ipPun].ipPnunit,
"\t%s\t%.3e", label.c_str(), max_odep );
fprintf( save.params[ipPun].ipPnunit, "\n" );
j += nskip;
} 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