sampling SED

This commit is contained in:
caes 2017-08-08 06:36:00 -04:00
parent 53aebdea81
commit 6e91480647
2 changed files with 49 additions and 48 deletions

View File

@ -14,7 +14,6 @@ int main(int argc, char const *argv[])
int n = 1000; int n = 1000;
agn::sed_table SED; agn::sed_table SED;
agn::sed_table samples; agn::sed_table samples;
agn::sed_spline agnsource;
const char* sample_filename = argv[1]; const char* sample_filename = argv[1];
const char* output_filename = argv[2]; const char* output_filename = argv[2];
@ -38,14 +37,12 @@ int main(int argc, char const *argv[])
// Read in sampling table and construct a spline model. // Read in sampling table and construct a spline model.
samples = agn::read_sed_table(sample_table); samples = agn::read_sed_table(sample_table);
agnsource = agn::sed_spline(samples); agn::sed_spline agnsource(samples);
if(agn::debug) debug_file { if(agn::debug) debug_file
std::cout
<< "Read samples:\n" << "Read samples:\n"
<< format_sed_table(samples); << format_sed_table(samples);
}
if(agn::verbose) std::cout if(agn::verbose) std::cout
<< "Evaluating relative spectral intensity for " << "Evaluating relative spectral intensity for "
@ -61,11 +58,6 @@ int main(int argc, char const *argv[])
output_table << agn::format_sed_table(SED); output_table << agn::format_sed_table(SED);
if(agn::verbose) std::cout
<< "Printing CLOUDY interpolate command syntax to file "
<< cloudyscript_filename
<< "\n";
if(agn::verbose) std::cout if(agn::verbose) std::cout
<< "Closing files. Goodbye.\n"; << "Closing files. Goodbye.\n";

View File

@ -36,23 +36,27 @@ const double IN_EV_2500A=12398.41929/2500;
// SEDs are represented by 2d histogram tables. // SEDs are represented by 2d histogram tables.
struct sed_table { struct sed_table {
std::string header; std::string header;
table_1d value; table1d value;
}; };
class sed { class sed {
public:
// Continuum output functions // Continuum output functions
// Returns histogram with n bins evenly space in log space // Returns histogram with n bins evenly space in log space
sed_table histogram_table(int n); sed_table histogram_table(int n);
// Argument is photon energy in eV // Argument is photon energy in eV
virtual double sed(double hnu); virtual double value(double hnu) {};
}
class sed_spline : sed { sed() {};
};
class sed_spline : public sed {
private: private:
Spline _spline; Spline<double,double> _spline;
public: public:
double value(double hnu);
sed_spline(agn::sed_table& samples); sed_spline(agn::sed_table& samples);
// These parameters might still be useful for rolling off various quantities, but aren't used in the strict-spline case. // These parameters might still be useful for rolling off various quantities, but aren't used in the strict-spline case.
@ -75,8 +79,9 @@ public:
double _xray_coefficient; double _xray_coefficient;
}; };
class sed_pow_law : sed { class sed_pow_law : public sed {
public: public:
double value(double hnu);
// Argument is photon energy in eV // Argument is photon energy in eV
double eval_uv(double hnu); double eval_uv(double hnu);
double eval_xray(double hnu); double eval_xray(double hnu);
@ -117,18 +122,39 @@ public:
); );
}; };
// Returns coord in eV for given relative coord.
double hnu_at(int i,int n);
// Takes an SED table as input and returns a string with format:
// '<h*nu>\t<flux>\n' for each energy-flux pair
std::string format_sed_table(sed_table table);
// Read continuum from file with '<h*nu>\t<flux>\n' formatting.
// Will ignore up to 1 header.
sed_table read_sed_table(std::ifstream& table_file);
// Does the same but converts hnu from rydberg to eV.
sed_table read_and_convert_sed_table(std::ifstream& table_file);
// Cloudy takes the SED density as input. This function outputs
// the corresponding SED table's SED density function in the form
// of a cloudy input script "interpolate" command.
std::string cloudy_interpolate_str(sed_table SED);
} // end namespace agn
// Constructors // Constructors
agn::sed_spline::sed_spline(agn::sed_table& samples) { agn::sed_spline::sed_spline(agn::sed_table& samples) {
std::vector<double> x; std::vector<double> x;
std::vector<double> y; std::vector<double> y;
iterator2d table_it = samples.begin(); iterator1d table_it = samples.value.begin();
while(table_it != samples.end()) { while(table_it != samples.value.end()) {
x.push(table_it->first); x.push_back(table_it->first);
y.push(table_it->second); y.push_back(table_it->second);
} }
Spline newspline(x,y); Spline<double,double> newspline(x,y);
_spline = newspline; _spline = newspline;
} }
@ -167,18 +193,18 @@ agn::sed_table agn::sed::histogram_table(int n){
double max=0,min=1,hnu; double max=0,min=1,hnu;
for(int i=0; i<n; i++) { for(int i=0; i<n; i++) {
hnu = hnu_at(i,n); hnu = hnu_at(i,n);
output.value[hnu] = this->sed(hnu); output.value[hnu] = this->value(hnu);
if (output.value[hnu] > max) max = output.value[hnu]; if (output.value[hnu] > max) max = output.value[hnu];
if (output.value[hnu] < min) min = output.value[hnu]; if (output.value[hnu] < min) min = output.value[hnu];
} }
// Add a final point at 100 KeV // Add a final point at 100 KeV
hnu = 1e5; hnu = 1e5;
output.value[hnu] = this->sed(hnu); output.value[hnu] = this->value(hnu);
return output; return output;
} }
// sed_spline evaluation // sed_spline evaluation
double agn::sed_spline::sed(double hnu) { double agn::sed_spline::value(double hnu) {
double magnitude=0.0; double magnitude=0.0;
magnitude += this->_spline[hnu]; magnitude += this->_spline[hnu];
if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL;
@ -186,7 +212,7 @@ double agn::sed_spline::sed(double hnu) {
} }
// sed_pow_law evaluations // sed_pow_law evaluations
double agn::sed_pow_law::sed(double hnu) { double agn::sed_pow_law::value(double hnu) {
double magnitude=0.0; double magnitude=0.0;
magnitude += this->eval_uv(hnu); magnitude += this->eval_uv(hnu);
magnitude += this->eval_xray(hnu); magnitude += this->eval_xray(hnu);
@ -220,29 +246,12 @@ double agn::sed_pow_law::SED_at_2KeV() {
// Utilities // Utilities
// Returns coord in eV for given relative coord.
double hnu_at(int i,int n);
// Takes an SED table as input and returns a string with format:
// '<h*nu>\t<flux>\n' for each energy-flux pair
std::string format_sed_table(sed_table table);
// Read continuum from file with '<h*nu>\t<flux>\n' formatting.
// Will ignore up to 1 header.
sed_table read_sed_table(std::ifstream& table_file);
// Does the same but converts hnu from rydberg to eV.
sed_table read_and_convert_sed_table(std::ifstream& table_file);
// Cloudy takes the SED density as input. This function outputs
// the corresponding SED table's SED density function in the form
// of a cloudy input script "interpolate" command.
std::string cloudy_interpolate_str(sed_table SED);
} // end namespace agn
agn::sed_table agn::read_sed_table(std::ifstream& table_file) { agn::sed_table agn::read_sed_table(std::ifstream& table_file) {
sed_table resultant; sed_table resultant;
std::string scratch; std::string scratch;
@ -284,7 +293,7 @@ std::string agn::format_sed_table(agn::sed_table table) {
std::stringstream output; std::stringstream output;
if (!table.header.empty()) output << table.header; if (!table.header.empty()) output << table.header;
output << std::setprecision(5); output << std::setprecision(5);
agn::table2d::iterator table_iterator; agn::table1d::iterator table_iterator;
table_iterator=table.value.begin(); table_iterator=table.value.begin();
while(table_iterator != table.value.end()) { while(table_iterator != table.value.end()) {
output output
@ -301,7 +310,7 @@ std::string agn::format_sed_table(agn::sed_table table) {
std::string agn::cloudy_interpolate_str(agn::sed_table table) { std::string agn::cloudy_interpolate_str(agn::sed_table table) {
std::stringstream output; std::stringstream output;
agn::table2d::iterator table_iterator = table.value.begin(); agn::table1d::iterator table_iterator = table.value.begin();
// Lead in to uv bump at slope=2 in log(energy [rydberg]) space // Lead in to uv bump at slope=2 in log(energy [rydberg]) space
double energy_in_rydbergs = table_iterator->first double energy_in_rydbergs = table_iterator->first
/ agn::RYDBERG_UNIT_EV; / agn::RYDBERG_UNIT_EV;