Otho/industry-website
1
0
Fork 0
mirror of https://asciireactor.com/otho/industry-website.git synced 2024-11-25 12:15:07 +00:00
Code Issues Packages Projects Releases Wiki Activity
0dd8354144
industry-website/htdocs/examples/spectral_lines.hpp

495 lines
15 KiB
C++
Raw Blame History

#ifndef spectral_lines_hpp
#define spectral_lines_hpp
#include "agn.hpp"
//#include "png++-0.2.5/gray_pixel.hpp"
//#include "png++-0.2.5/image.hpp"
namespace agn {
const double EQWIDTH_MIN_VAL = 1E-7;
const double EQWIDTH_MIN_VAL_LOG = -7;
const double SIGMA_THRESHOLD_MULTIPLIER = 1.0;
const double RATIO_THRESHOLD_MULTIPLIER = .06;
typedef std::list<std::string> line_list;
line_list read_line_list(std::ifstream&);
// Trims whitespace from line labels
std::string label_trim(std::string);
// Emission line flux contours are represented by a singly-nested sorted map.
struct eqwidth_table {
std::string header[2];
table2d value;
};
bool is_zero(const eqwidth_table&);
eqwidth_table read_eqwidth_table(std::ifstream& table_file);
std::string format_eqwidth_table(eqwidth_table table);
std::string format_eqwidth_table_slice(eqwidth_table table,iterator2d x);
gridcoordlist find_outliers(eqwidth_table table);
gridcoordlist known_outliers();
eqwidth_table smooth(eqwidth_table&, gridcoordlist coordlist);
std::list<eqwidth_table> compile_eqwidth_tables(cloudy_grid,line_list,double = 1.0);
// Operator<< prints table formatted for a "fort file".
std::ostream& operator<< (std::ostream&,eqwidth_table);
} // end namespace agn
// Definitions
bool agn::is_zero(const eqwidth_table& table_to_check) {
agn::table2d::const_iterator x_it = table_to_check.value.begin();
while (x_it != table_to_check.value.end()) {
agn::table1d::const_iterator y_it= x_it->second.begin();
while (y_it != x_it->second.end()) {
if (y_it->second != EQWIDTH_MIN_VAL)
return false;
y_it++;
}
x_it++;
}
return true;
}
agn::eqwidth_table agn::read_eqwidth_table(std::ifstream& table_file) {
agn::eqwidth_table resultant;
double hden, phi, flux;
int number_of_lines = 0;
std::string null;
while (std::getline(table_file, null))
++number_of_lines;
int table_entries = number_of_lines - 2;
table_file.clear();
table_file.seekg(0);
getline(table_file,resultant.header[0]);
getline(table_file,resultant.header[1]);
for (int i=0; i<table_entries; i++) {
table_file >> hden >> phi >> flux;
resultant.value[hden][phi]=flux;
}
return resultant;
}
std::string agn::format_eqwidth_table(agn::eqwidth_table table) {
std::stringstream output;
output
<< table.header[0]
<< "\n"
<< table.header[1]
<< "\n";
agn::iterator2d hden_it = table.value.begin();
while(hden_it != table.value.end()) {
agn::iterator1d phi_it = table.value[hden_it->first].begin();
while(phi_it != table.value[hden_it->first].end()) {
output
<< std::fixed
<< std::setprecision(3)
<< hden_it->first
<< " "
<< phi_it->first
<< " "
<< std::scientific
<< std::setprecision(4)
<< phi_it->second
<< "\n";
phi_it++;
}
hden_it++;
}
return output.str();
}
std::string agn::format_eqwidth_table_slice(agn::eqwidth_table table,agn::iterator2d fixed) {
std::stringstream output("");
// output
// << table.header[0]
// << " - constant="
// << fixed->first
// << "\n";
iterator1d x=fixed->second.begin();
while(x!=fixed->second.end()) {
output
<< std::fixed
<< std::setprecision(3)
<< x->first
<< "\t"
<< std::scientific
<< std::setprecision(4)
<< x->second
<< "\n";
x++;
}
return output.str();
}
std::ostream& agn::operator<< ( std::ostream& outstream,
agn::eqwidth_table table) {
std::string header0_trimmed = agn::label_trim(table.header[0]);
outstream
<< header0_trimmed
<< std::endl
<< table.header[1]
<< std::endl;
agn::table2d::iterator x_it = table.value.begin();
while (x_it != table.value.end()) {
agn::table1d::iterator y_it= x_it->second.begin();
while (y_it != x_it->second.end()) {
outstream.width(6);
outstream
<< std::fixed
<< std::setprecision(3)
<< x_it->first
<< "\t"
<< y_it->first
<< "\t"
<< std::scientific
<< std::setprecision(4)
<< y_it->second
<< "\n";
y_it++;
}
x_it++;
}
return outstream;
}
// Returns a table as the difference between two tables
//agn::eqwidth_table main(int argc, char const *argv[])
//{ // Copied from fortfile subtractor
// std::ifstream table1;
// table1.open(argv[1]);
// std::ifstream table2;
// table2.open(argv[2]);
//
// std::string line;
//
// int table1_numlines = 0;
// while (std::getline(table1, line))
// ++table1_numlines;
// table1.clear();
// table1.seekg(0);
//
// int table_entries = table1_numlines - 2;
//
//// int table2_numlines = 0;
//// while (std::getline(table1, line))
//// ++table2_numlines;
////
//// int table2_entries = table2_numlines - 2;
//// table2.clear();
//// table2.seekg(0);
//
// double hden1[table_entries];
// double phi1[table_entries];
// double magnitude1[table_entries];
// double resultant_magnitude[table_entries];
//
// std::string first_line, second_line, scratch;
// getline(table1,first_line);
// getline(table1,second_line);
// getline(table2,scratch);
// getline(table2,scratch);
//
// first_line[0] = argv[3][0];
// first_line[1] = argv[3][1];
// first_line[2] = argv[3][2];
// first_line[3] = argv[3][3];
//
// for (int i=0; i<table_entries; i++) {
// table1 >> hden1[i] >> phi1[i] >> magnitude1 [i];
// }
//
//
// int i=0;
// double hden2, phi2, magnitude2;
// for (int i=0; i<table_entries; i++) {
// if ( table2.eof() ) {
// resultant_magnitude[i]=magnitude1[i];
// }
// else {
// table2 >> hden2 >> phi2 >> magnitude2;
// resultant_magnitude[i] = magnitude1[i] - magnitude2;
// }
// if ( resultant_magnitude[i] < 1e-10 )
// resultant_magnitude[i]=1e-10;
// }
//
//
//
// std::cout
// << first_line
// << std::endl
// //<< "HDen Phi(H) ScaInten"
// << second_line
// << std::endl;
//
// for (int i=0; i<table_entries; i++) {
// std::cout
// << std::fixed
// << std::setprecision(3)
// << hden1[i]
// << " "
// << phi1[i]
// << " "
// << std::scientific
// << std::setprecision(4)
// << resultant_magnitude[i]
// << std::endl;
// }
//
// table1.close();
// table2.close();
// return 0;
//}
agn::line_list agn::read_line_list(std::ifstream& inputfile) {
line_list list;
line_list nfnulist;
line_list inwtlist;
line_list inwclist;
inputfile.clear();
inputfile.seekg(0);
std::string line,label;
while(!inputfile.eof()) {
getline(inputfile,line);
if (line[0] == '#' || line.size() < 2) continue;
if (line[0] == ' ' && line[1] == ' ') continue;
if (line[0] == ' ')
label = line.substr(1,agn::LABELSTR_LEN);
else
label = line.substr(0,agn::LABELSTR_LEN);
// above is only good for c17?
if(agn::debug) std::cout
<< "Adding label "
<< label
<< " to capture list.\n";
std::string applabel = label.substr(4);
if (label.substr(0,4).compare(std::string("nFnu")) == 0) {
nfnulist.push_back(std::string("nFnu").append(applabel));
//list.push_back(std::string("nInu").append(applabel));
inwtlist.push_back(std::string("InwT").append(applabel));
inwclist.push_back(std::string("InwC").append(applabel));
continue;
}
list.push_back(label);
//if (label.substr(0,4) != "Inci" && applabel != " 0 ")
// list.push_back(std::string("Inwd").append(applabel));
}
list.insert(list.end(),nfnulist.begin(),nfnulist.end());
list.insert(list.end(),inwtlist.begin(),inwtlist.end());
list.insert(list.end(),inwclist.begin(),inwclist.end());
return list;
}
std::list<agn::eqwidth_table> agn::compile_eqwidth_tables(agn::cloudy_grid grid,agn::line_list line_labels,double scale_factor) {
std::list<eqwidth_table> table_list_eq_width,table_list_radiated_energy;
agn::line_list::iterator line_label_it = line_labels.begin();
// Iterate over labels, creating table for each:
while (line_label_it != line_labels.end()) {
double x,y;
agn::cloudy_line_data data;
agn::eqwidth_table new_table;
std::string label = *line_label_it;
if(agn::debug) std::cout
<< "Processing label "
<< label
<< std::endl;
std::stringstream header0;
header0
<< label;
new_table.header[0] = header0.str();
new_table.header[1] = "Hden Phi(H) F(line)/Inci(1215) ";
agn::cloudy_grid::iterator result_it = grid.begin();
while (result_it != grid.end()) {
x = result_it->first.first;
y = result_it->first.second;
// if no data exist under this label, initialize some:
if ( result_it->second.intrinsic_line_intensity.count(label) == 0 ) {
data.radiated_energy=-35.0;
data.eq_width=0.0;
}
else {
data = result_it->second.intrinsic_line_intensity[label];
data.eq_width /= scale_factor;
if(agn::line_debug) std::cout
<< "Added "
<< std::setprecision(2)
<< std::fixed
<< x
<< ", "
<< y
<< " -- "
<< std::scientific
<< data.radiated_energy
<< ", "
<< data.eq_width
<< "\n";
if ( data.has_duplicates ) {
// This needs to be expanded to catch all duplicate values j=1,2,etc.
// For now, just blindly use the value the program finds...
// the functionality is written in the agn.hpp sections, but is not
// implemented here, and as of 2017, the agn.hpp function writes the
// last-encountered instance of the emission line value to the label
// probably needs to be fixed for the inwd part too
}
}
if ( data.radiated_energy < EQWIDTH_MIN_VAL_LOG || data.eq_width < EQWIDTH_MIN_VAL) {
data.radiated_energy = EQWIDTH_MIN_VAL_LOG;
data.eq_width = EQWIDTH_MIN_VAL;
}
new_table.value[x][y] = data.eq_width;
result_it++;
}
table_list_eq_width.push_back(new_table);
// Add inward table for emission lines, and only works for c17
if (label.substr(agn::LABELSTR_POS,4) != "Inci" &&
label.substr(agn::LABELSTR_POS,4) != "FeKa" &&
//label.substr(agn::LABELSTR_POS,4) != "Blnd" &&
label.substr(agn::LABELSTR_POS,3) != "nFn" &&
label.substr(agn::LABELSTR_POS,3) != "Inw" &&
label.substr(LABELSTR_LEN - 14) != " 0 ") {
// These lines are known to not have inward fractions, so
// should be skipped.
if (label == "Blnd 1657.00A" ||
label == "Blnd 2326.00A" ||
label == "Blnd 2798.00A" ||
label == "Blnd 1335.00A" ||
label == "Blnd 977.000A" ||
label == "Blnd 1909.00A" ||
label == "Blnd 1176.00A" ||
label == "Blnd 2141.00A" ||
label == "Blnd 1085.00A" ||
label == "Blnd 1750.00A" ||
label == "Blnd 990.000A" ||
label == "Blnd 1486.00A" ||
label == "Blnd 765.000A" ||
label == "Blnd 1304.00A" ||
label == "Blnd 8446.00A" ||
label == "Blnd 1666.00A" ||
label == "Blnd 835.000A" ||
label == "Blnd 1402.00A" ||
label == "Blnd 630.000A" ||
label == "Blnd 1218.00A" ||
label == "Blnd 774.000A" ||
label == "Blnd 5892.00A" ||
label == "Blnd 615.000A" ||
label == "Blnd 2665.00A" ||
label == "Blnd 1860.00A" ||
label == "Blnd 2335.00A" ||
label == "Blnd 1888.00A" ||
label == "Blnd 1397.00A" ||
label == "Blnd 1256.00A" ||
label == "Blnd 1720.00A" ||
label == "Blnd 1198.00A" ||
label == "Blnd 1406.00A" ||
label == "Blnd 1199.00A" ||
label == "Blnd 8579.00A" ||
label == "Blnd 3933.00A"
) {
line_label_it++;
continue;
}
if(agn::debug) std::cout
<< "Processing Inwd for label "
<< label
<< std::endl;
std::stringstream inwdheader0;
agn::eqwidth_table inwd_table;
inwdheader0
<< "Inwd "
<< label.substr(5);
inwd_table.header[0] = inwdheader0.str();
inwd_table.header[1] = "Hden Phi(H) F(line)/Inci(1215)";
result_it = grid.begin();
agn::cloudy_line_data inwddata;
while (result_it != grid.end()) {
x = result_it->first.first;
y = result_it->first.second;
// if no data exist under this label, initialize some:
if ( result_it->second.intrinsic_line_inward_intensity.count(label) == 0 ) {
inwddata.radiated_energy=-35.0;
inwddata.eq_width=0.0;
}
else {
inwddata = result_it->second.intrinsic_line_inward_intensity[label];
inwddata.eq_width /= scale_factor;
}
if(agn::line_debug) std::cout
<< "Added "
<< std::setprecision(2)
<< std::fixed
<< x
<< ", "
<< y
<< " -- "
<< std::scientific
<< inwddata.radiated_energy
<< ", "
<< inwddata.eq_width
<< "\n";
if ( inwddata.has_duplicates ) {
// This needs to be expanded to catch all duplicate values j=1,2,etc.
// For now, just blindly use the value the program finds...
// the functionality is written in the agn.hpp sections, but is not
// implemented here, and as of 2017, the agn.hpp function writes the
// last-encountered instance of the emission line value to the label
// probably needs to be fixed for the inwd part too
}
if ( inwddata.radiated_energy < EQWIDTH_MIN_VAL_LOG ||
inwddata.eq_width < EQWIDTH_MIN_VAL) {
inwddata.radiated_energy = EQWIDTH_MIN_VAL_LOG;
inwddata.eq_width = EQWIDTH_MIN_VAL;
}
inwd_table.value[x][y] = inwddata.eq_width;
result_it++;
}
table_list_eq_width.push_back(inwd_table);
}
line_label_it++;
}
return table_list_eq_width;
}
#endif
Reference in New Issue View Git Blame Copy Permalink