Added powerlaw to spline model

This commit is contained in:
caes 2017-08-08 09:07:34 -04:00
parent a530c8054f
commit 9e325673c1
2 changed files with 30 additions and 23 deletions

View File

@ -8,14 +8,6 @@
namespace agn { namespace agn {
// Continuum domain, step size constant in log space
const double CONT_MIN_ENERGY=1e-2; // eV
const double CONT_MAX_ENERGY=1e5; // eV
const double CONT_MIN_X=log10(CONT_MIN_ENERGY);
const double CONT_MAX_X=log10(CONT_MAX_ENERGY);
const double CONT_WIDTH_X=CONT_MAX_X - CONT_MIN_X;
const double CONT_MIN_VAL=1e-35;
// Cloudy's continuum domain, for reference // Cloudy's continuum domain, for reference
// Pulled from cloudy 17.00, first version // Pulled from cloudy 17.00, first version
// rfield.emm = 1.001e-8f; // rfield.emm = 1.001e-8f;
@ -25,14 +17,18 @@ const double CLOUDY_EGAMRY = 7.354e6; // in Rydberg
const double CLOUDY_MIN_EV=CLOUDY_EMM*RYDBERG_UNIT_EV; const double CLOUDY_MIN_EV=CLOUDY_EMM*RYDBERG_UNIT_EV;
const double CLOUDY_MAX_EV=CLOUDY_EGAMRY*RYDBERG_UNIT_EV; const double CLOUDY_MAX_EV=CLOUDY_EGAMRY*RYDBERG_UNIT_EV;
// Continuum domain, step size constant in log space
const double CONT_MIN_ENERGY=1e-2; // eV
const double CONT_MAX_ENERGY=1e5; // eV
const double CONT_MIN_LOGX=log10(CONT_MIN_ENERGY);
const double CONT_MAX_LOGX=log10(CONT_MAX_ENERGY);
const double CONT_WIDTH_LOGX=CONT_MAX_LOGX - CONT_MIN_LOGX;
const double CONT_MIN_VAL=1e-35;
const double IN_EV_2500A=12398.41929/2500; 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;
@ -52,12 +48,16 @@ public:
sed() {}; sed() {};
}; };
class sed_spline : public sed { class sed_powerlaw_spline : public sed {
private: private:
Spline<double,double> _spline; Spline<double,double> _spline;
public: public:
double value(double hnu); double value(double hnu);
sed_spline(agn::sed_table& samples); sed_powerlaw_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.
@ -146,7 +146,7 @@ std::string cloudy_interpolate_str(sed_table SED);
// Constructors // Constructors
agn::sed_spline::sed_spline(agn::sed_table& samples) { agn::sed_powerlaw_spline::sed_powerlaw_spline(agn::sed_table& samples) {
std::vector<double> x; std::vector<double> x;
std::vector<double> y; std::vector<double> y;
iterator1d table_it = samples.value.begin(); iterator1d table_it = samples.value.begin();
@ -204,8 +204,8 @@ agn::sed_table agn::sed::histogram_table(int n){
return output; return output;
} }
// sed_spline evaluation // sed_powerlaw_spline evaluation
double agn::sed_spline::value(double hnu) { double agn::sed_powerlaw_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;
@ -368,7 +368,7 @@ std::string agn::cloudy_interpolate_str(agn::sed_table table) {
double agn::hnu_at(int i,int n) { double agn::hnu_at(int i,int n) {
double relative_coord=(double)(i)/n; double relative_coord=(double)(i)/n;
double x_coord = relative_coord*CONT_WIDTH_X + CONT_MIN_X; double x_coord = relative_coord*CONT_WIDTH_LOGX + CONT_MIN_LOGX;
return pow(10,x_coord); return pow(10,x_coord);
} }

View File

@ -9,20 +9,26 @@ int main(int argc, char const *argv[])
<< "Setting up environment.\n"; << "Setting up environment.\n";
// Create 2d table using n bins, linear values of SED. The // Create 2d table using n bins, linear values of SED. The
// agn sed_spline class has a function for this. A // agn sed_powerlaw_spline class has a function for this. A
// std::map<double,double> represents the table. // std::map<double,double> represents the table.
int n = 1000; int n = 1000;
agn::sed_table SED; agn::sed_table SED;
agn::sed_table samples; agn::sed_table samples;
agn::sed_table powerlaw_coords;
const char* sample_filename = argv[1]; const char* sample_filename = argv[1];
const char* output_filename = argv[2]; const char* powerlaw_filename = argv[2];
const char* output_filename = argv[3];
const char* debug_filename = "spline_sed_debug"; const char* debug_filename = "spline_sed_debug";
std::ifstream sample_table( std::ifstream sample_table(
sample_filename, sample_filename,
std::ofstream::out std::ofstream::out
); );
std::ifstream powerlaw_table(
sample_filename,
std::ofstream::out
);
std::ofstream output_table( std::ofstream output_table(
output_filename, output_filename,
std::ofstream::out std::ofstream::out
@ -37,12 +43,13 @@ 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);
powerlaw_coords = agn::read_sed_table(powerlaw_table);
if(agn::debug) debug_file if(agn::debug) debug_file
<< "Read samples:\n" << "Read samples:\n"
<< format_sed_table(samples); << format_sed_table(samples);
agn::sed_spline agnsource(samples); agn::sed_powerlaw_spline agnsource(samples,powerlaw_coords);