diff --git a/src/generate_powlaw_sed.cpp b/src/generate_powlaw_sed.cpp index b02fbfe..f0c5e9a 100644 --- a/src/generate_powlaw_sed.cpp +++ b/src/generate_powlaw_sed.cpp @@ -1,6 +1,7 @@ #include "agn.hpp" #include "sed.hpp" +// This is old and probably doesn't work. int main(int argc, char const *argv[]) { diff --git a/src/sed.hpp b/src/sed.hpp index 7b6f9f1..3fa43b7 100644 --- a/src/sed.hpp +++ b/src/sed.hpp @@ -32,9 +32,25 @@ const double IN_EV_2500A=12398.41929/2500; // SEDs are represented by 2d histogram tables. struct sed_table { std::string header; - table1d value; + table1d table; }; +class powerlaw { +private: + // f(x) = _normal*x^_power + double _power; + double _normal; +public: + powerlaw(): _power(0), _normal(0) {} + powerlaw(coord2d x0,coord2d x1): + _power((log(x1.second)-log(x0.second))/(log(x1.first)-log(x0.first))), + _normal((log(x0.second)-(_power*log(x0.first)))) + {} + powerlaw(coord2d x0,double slope): + _power(slope), + _normal((log(x0.second)-(_power*log(x0.first)))) + {} +}; class sed { public: @@ -43,30 +59,18 @@ public: sed_table histogram_table(int n); // Argument is photon energy in eV - virtual double value(double hnu) {}; + virtual double table(double hnu) {}; sed() {}; }; class sed_powerlaw_spline : public sed { private: - Spline _spline; - - // powerlaw parameters - double _ir_slope = 3; - double _ir_high_point_x; - double _ir_high_point_y; - double _uv_low_point_x; - double _uv_low_point_y; - double _uv_high_point_x; - double _uv_high_point_y; - double _xray_low_point_x; - double _xray_low_point_y; - double _xray_high_point_x; - double _xray_high_point_y; - double _gamma_low_point_x; - double _gamma_low_point_y; - double _gamma_slope = -2; + Spline _output_model; + powerlaw _ir_powerlaw; + powerlaw _uv_powerlaw; + powerlaw _xray_powerlaw; + powerlaw _gamma_powerlaw; // These parameters might still be useful for rolling off various quantities, but aren't used in the strict-spline case. @@ -79,14 +83,14 @@ private: double _xray_coefficient; public: - double value(double hnu); + double table(double hnu); sed_powerlaw_spline(agn::sed_table& samples, agn::sed_table& powerlaw_coords); }; class sed_pow_law : public sed { public: - double value(double hnu); + double table(double hnu); // Argument is photon energy in eV double eval_uv(double hnu); double eval_xray(double hnu); @@ -154,40 +158,42 @@ std::string cloudy_interpolate_str(sed_table SED); agn::sed_powerlaw_spline::sed_powerlaw_spline( agn::sed_table& samples, agn::sed_table& powerlaw_coords - ) { - std::vector x; - std::vector y; - samples - iterator1d table_it = powerlaw_coords.value.begin(); - _ir_high_point_x = table_it->first; - _ir_high_point_y = table_it->second; - table_it++; - _uv_high_point_x = table_it->first; - _uv_high_point_y = table_it->second; - table_it++; - _uv_high_point_x = table_it->first; - _uv_high_point_y = table_it->second; - table_it++; - _xray_high_point_x = table_it->first; - _xray_high_point_y = table_it->second; - table_it++; - _xray_high_point_x = table_it->first; - _xray_high_point_y = table_it->second; - table_it++; - _gamma_high_point_x = table_it->first; - _gamma_high_point_y = table_it->second; + ) +{ + // coordinate vectors will be used to construct spline sed model + std::vector x0; + std::vector x1; + + // powerlaws are evaluated across four regions of the sed, first + iterator1d table_it = powerlaw_coords.table.begin(); + double ir_power = 3; + coord2d ir_high_point = *table_it; table_it++; + coord2d uv_low_point = *table_it; table_it++; + coord2d uv_high_point = *table_it; table_it++; + coord2d xray_low_point = *table_it; table_it++; + coord2d xray_high_point = *table_it; table_it++; + coord2d gamma_low_point = *table_it; + double gamma_power = -2; + _ir_powerlaw = powerlaw(ir_high_point,ir_power); + _uv_powerlaw = powerlaw(uv_low_point,uv_high_point); + _xray_powerlaw = powerlaw(xray_low_point,xray_high_point); + _gamma_powerlaw = powerlaw(gamma_low_point,gamma_power); + + ir_bounds= + for (int i=0; i<10; i++) { - _ - agn::coord2d uv_point = + _ir_powerlaw.eval() } - table_it = samples.value.begin(); - while(table_it != samples.value.end()) { - x.push_back(table_it->first); - y.push_back(table_it->second); + + // load all samples into coordinate vectors + table_it = samples.table.begin(); + while(table_it != samples.table.end()) { + x0.push_back(table_it->first); + x1.push_back(table_it->second); table_it++; } - Spline newspline(x,y); - _spline = newspline; + Spline newspline(x0,x1); + _output_model = newspline; } agn::sed_pow_law::sed_pow_law ( @@ -208,12 +214,12 @@ agn::sed_pow_law::sed_pow_law ( _cutoff_xray_rydberg(cutoff_xray_rydberg), _log_radius_in_cm(log_radius_in_cm), _scaling_factor(scaling_factor) - { - _cutoff_uv_eV = cutoff_uv_rydberg*RYDBERG_UNIT_EV; - _cutoff_xray_eV = cutoff_xray_rydberg*RYDBERG_UNIT_EV; - _radius_in_cm = pow(10,log_radius_in_cm); - _radius_in_cm_squared = _radius_in_cm*_radius_in_cm; - _xray_coefficient = agn::sed_pow_law::SED_at_2KeV(); +{ + _cutoff_uv_eV = cutoff_uv_rydberg*RYDBERG_UNIT_EV; + _cutoff_xray_eV = cutoff_xray_rydberg*RYDBERG_UNIT_EV; + _radius_in_cm = pow(10,log_radius_in_cm); + _radius_in_cm_squared = _radius_in_cm*_radius_in_cm; + _xray_coefficient = agn::sed_pow_law::SED_at_2KeV(); } @@ -225,26 +231,26 @@ agn::sed_table agn::sed::histogram_table(int n){ double max=0,min=1,hnu; for(int i=0; ivalue(hnu); - if (output.value[hnu] > max) max = output.value[hnu]; - if (output.value[hnu] < min) min = output.value[hnu]; + output.table[hnu] = this->table(hnu); + if (output.table[hnu] > max) max = output.table[hnu]; + if (output.table[hnu] < min) min = output.table[hnu]; } // Add a final point at 100 KeV hnu = 1e5; - output.value[hnu] = this->value(hnu); + output.table[hnu] = this->table(hnu); return output; } // sed_powerlaw_spline evaluation -double agn::sed_powerlaw_spline::value(double hnu) { +double agn::sed_powerlaw_spline::table(double hnu) { double magnitude=0.0; - magnitude += this->_spline[hnu]; + magnitude += this->_output_model[hnu]; if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; return magnitude; } // sed_pow_law evaluations -double agn::sed_pow_law::value(double hnu) { +double agn::sed_pow_law::table(double hnu) { double magnitude=0.0; magnitude += this->eval_uv(hnu); magnitude += this->eval_xray(hnu); @@ -296,7 +302,7 @@ agn::sed_table agn::read_sed_table(std::ifstream& table_file) { } while(!table_file.eof()) { table_file >> hnu; - table_file >> resultant.value[hnu]; + table_file >> resultant.table[hnu]; } return resultant; } @@ -306,7 +312,7 @@ agn::sed_table agn::read_and_convert_sed_table(std::ifstream& table_file) { sed_table resultant; std::string scratch; int current_line=0; - double hnu_in_ryd,hnu_in_ev,value; + double hnu_in_ryd,hnu_in_ev,table; std::getline(table_file,scratch); if(!isdigit(scratch[0])) { resultant.header = scratch; @@ -317,7 +323,7 @@ agn::sed_table agn::read_and_convert_sed_table(std::ifstream& table_file) { //std::cout << c; table_file >> hnu_in_ryd; hnu_in_ev = hnu_in_ryd*agn::RYDBERG_UNIT_EV; - table_file >> resultant.value[hnu_in_ev]; + table_file >> resultant.table[hnu_in_ev]; getline(table_file,scratch); } } @@ -327,8 +333,8 @@ std::string agn::format_sed_table(agn::sed_table table) { if (!table.header.empty()) output << table.header; output << std::setprecision(5); agn::table1d::iterator table_iterator; - table_iterator=table.value.begin(); - while(table_iterator != table.value.end()) { + table_iterator=table.table.begin(); + while(table_iterator != table.table.end()) { output << std::fixed << std::scientific @@ -344,7 +350,7 @@ std::string agn::format_sed_table(agn::sed_table table) { std::string agn::cloudy_interpolate_str(agn::sed_table table) { std::stringstream output; - agn::table1d::iterator table_iterator = table.value.begin(); + agn::table1d::iterator table_iterator = table.table.begin(); // Lead in to uv bump at slope=2 in log(energy [rydberg]) space double energy_in_rydbergs = table_iterator->first / agn::RYDBERG_UNIT_EV; @@ -365,7 +371,7 @@ std::string agn::cloudy_interpolate_str(agn::sed_table table) { << ")"; int count=0; - while(table_iterator != table.value.end()) { + while(table_iterator != table.table.end()) { energy_in_rydbergs = table_iterator->first / agn::RYDBERG_UNIT_EV; double log_SED_density = log10( table_iterator->second diff --git a/src/table_powerlaw_spline_sed.cpp b/src/table_powerlaw_spline_sed.cpp index d816710..7178faf 100644 --- a/src/table_powerlaw_spline_sed.cpp +++ b/src/table_powerlaw_spline_sed.cpp @@ -2,6 +2,9 @@ #include "sed.hpp" +// Syntax: table_powerlaw_spline + + int main(int argc, char const *argv[]) { @@ -9,8 +12,7 @@ int main(int argc, char const *argv[]) << "Setting up environment.\n"; // Create 2d table using n bins, linear values of SED. The - // agn sed_powerlaw_spline class has a function for this. A - // std::map represents the table. + // agn sed_powerlaw_spline class has a function for this. int n = 1000; agn::sed_table SED; agn::sed_table samples;