From 7d2b6eba3e3d2decb01ef07a7cf9475743c815ac Mon Sep 17 00:00:00 2001 From: caes Date: Tue, 8 Aug 2017 05:15:14 -0400 Subject: [PATCH] spline sed generation --- src/generate_spline_sed.cpp | 76 ++++++ src/sed.hpp | 364 +++++++++++++++++++++++++++++ src/sed/mehdipour/generate_sed.cpp | 90 ------- src/sed/mehdipour/powlaw-cutoff.pl | 2 +- src/sed/mehdipour/sed.hpp | 84 +++++++ 5 files changed, 525 insertions(+), 91 deletions(-) create mode 100644 src/generate_spline_sed.cpp create mode 100644 src/sed.hpp diff --git a/src/generate_spline_sed.cpp b/src/generate_spline_sed.cpp new file mode 100644 index 0000000..771d1e0 --- /dev/null +++ b/src/generate_spline_sed.cpp @@ -0,0 +1,76 @@ +#include "agn.hpp" + + + +int main(int argc, char const *argv[]) +{ + + if (agn::verbose) std::cout + << "Setting up environment.\n"; + + // Create 2d table using n bins, linear values of SED. The + // agn sed_spline class has a function for this. A + // std::map represents the table. + int n = 1000; + agn::sed_table SED; + agn::sed_table samples; + agn::sed_spline agnsource; + + const char* sample_filename = argv[1]; + const char* output_filename = argv[2]; + const char* debug_filename = "spline_sed_debug"; + + std::ifstream sample_table( + sample_filename, + std::ofstream::out + ); + std::ofstream output_table( + output_filename, + std::ofstream::out + ); + std::ofstream debug_file( + debug_filename, + std::ofstream::out + ); + + if (agn::verbose) std::cout + << "Creating agn sed object.\n"; + + // Read in sampling table and construct a spline model. + samples = agn::read_sed_table(sample_table); + agnsource = agn::sed_spline(samples); + + + if(agn::debug) debug_file { + std::cout + << "Read samples:\n" + << format_sed_table(samples); + } + + if(agn::verbose) std::cout + << "Evaluating relative spectral intensity for " + << n + << " photon energy bins.\n"; + + SED = agnsource.histogram_table(n); + + if(agn::verbose) std::cout + << "Printing SED table to file " + << output_filename + << "\n"; + + 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 + << "Closing files. Goodbye.\n"; + + debug_file.close(); + output_table.close(); + return 0; +} + diff --git a/src/sed.hpp b/src/sed.hpp new file mode 100644 index 0000000..a6c606b --- /dev/null +++ b/src/sed.hpp @@ -0,0 +1,364 @@ +#ifndef sed_hpp +#define sed_hpp + + + +#include "agn.hpp" +#include "spline.h" + +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 +// Pulled from cloudy 17.00, first version +// rfield.emm = 1.001e-8f; +// rfield.egamry = 7.354e6f; +const double CLOUDY_EMM = 1.001e-8; // in Rydberg +const double CLOUDY_EGAMRY = 7.354e6; // in Rydberg +const double CLOUDY_MIN_EV=CLOUDY_EMM*RYDBERG_UNIT_EV; +const double CLOUDY_MAX_EV=CLOUDY_EGAMRY*RYDBERG_UNIT_EV; + +const double IN_EV_2500A=12398.41929/2500; + + + + + + + +// SEDs are represented by 2d histogram tables. +struct sed_table { + std::string header; + table_1d value; +}; + + +class sed { + // Continuum output functions + // Returns histogram with n bins evenly space in log space + sed_table histogram_table(int n); + + // Argument is photon energy in eV + virtual double sed(double hnu); +} + +class sed_spline : sed { +private: + Spline _spline; +public: + 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. + + // Continuum shape arguments + double _T; //TCut + double _alpha_ox; + double _alpha_x; + double _alpha_uv; + double _cutoff_uv_rydberg; + double _cutoff_xray_rydberg; + double _log_radius_in_cm; + + // Derived values + double _cutoff_uv_eV; // IRCut + double _cutoff_xray_eV; // lowend_cutoff + double _radius_in_cm; + double _radius_in_cm_squared; + double _scaling_factor; + double _xray_coefficient; +}; + +class sed_pow_law : sed { +public: + // Argument is photon energy in eV + double eval_uv(double hnu); + double eval_xray(double hnu); + // Determined differently to be of use as the + // xray coefficient. + double SED_at_2KeV(); + + + // Continuum shape arguments + double _T; //TCut + double _alpha_ox; + double _alpha_x; + double _alpha_uv; + double _cutoff_uv_rydberg; + double _cutoff_xray_rydberg; + double _log_radius_in_cm; + + // Derived values + double _cutoff_uv_eV; // IRCut + double _cutoff_xray_eV; // lowend_cutoff + double _radius_in_cm; + double _radius_in_cm_squared; + double _scaling_factor; + double _xray_coefficient; + + sed_pow_law ( + double T, + double alpha_ox, + double alpha_x, + double alpha_uv, + double cutoff_uv_rydberg, + double cutoff_xray_rydberg, + double log_radius_in_cm, + double scaling_factor = 1.0 + + // EL[e] model scaling factor + // double scaling_factor = 1.39666E44 + ); +}; + + +// Constructors + +agn::sed_spline::sed_spline(agn::sed_table& samples) { + std::vector x; + std::vector y; + iterator2d table_it = samples.begin(); + while(table_it != samples.end()) { + x.push(table_it->first); + y.push(table_it->second); + } + Spline newspline(x,y); + _spline = newspline; +} + +agn::sed_pow_law::sed_pow_law ( + double T, + double alpha_ox, + double alpha_x, + double alpha_uv, + double cutoff_uv_rydberg, + double cutoff_xray_rydberg, + double log_radius_in_cm, + double scaling_factor + ): + _T(T), + _alpha_ox(alpha_ox), + _alpha_x(alpha_x), + _alpha_uv(alpha_uv), + _cutoff_uv_rydberg(cutoff_uv_rydberg), + _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(); +} + + +// Class Functions + +// writes log-space histogram with n data +agn::sed_table agn::sed::histogram_table(int n){ + agn::sed_table output; + double max=0,min=1,hnu; + for(int i=0; ised(hnu); + if (output.value[hnu] > max) max = output.value[hnu]; + if (output.value[hnu] < min) min = output.value[hnu]; + } + // Add a final point at 100 KeV + hnu = 1e5; + output.value[hnu] = this->sed(hnu); + return output; +} + +// sed_spline evaluation +double agn::sed_spline::sed(double hnu) { + double magnitude=0.0; + magnitude += this->_spline[hnu]; + if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; + return magnitude; +} + +// sed_pow_law evaluations +double agn::sed_pow_law::sed(double hnu) { + double magnitude=0.0; + magnitude += this->eval_uv(hnu); + magnitude += this->eval_xray(hnu); + + return magnitude; +} +double agn::sed_pow_law::eval_uv(double hnu) { + double bigbump_kT = _T + * agn::BOLTZMANN_CONST; + double magnitude = pow(hnu,(1+_alpha_uv)) + * exp(-(hnu)/bigbump_kT) + * exp(-(_cutoff_uv_eV/hnu)) + * _scaling_factor; + if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; + return magnitude; +} +double agn::sed_pow_law::eval_xray(double hnu) { + return _xray_coefficient + * pow(hnu/2000,1+_alpha_x) + * exp(-_cutoff_xray_eV/hnu) + * _scaling_factor; +} +double agn::sed_pow_law::SED_at_2KeV() { + double ELe_at_2500A_no_scale = eval_uv(IN_EV_2500A) + / _scaling_factor; + double energy_ratio = 2000/IN_EV_2500A; + // Returns EL[e] at 2 KeV + return ELe_at_2500A_no_scale + * pow(energy_ratio,_alpha_ox + 1); +} + + + +// 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: +// '\t\n' for each energy-flux pair +std::string format_sed_table(sed_table table); + +// Read continuum from file with '\t\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) { + sed_table resultant; + std::string scratch; + int current_line=0; + double hnu; + std::getline(table_file,scratch); + if(!isdigit(scratch[0])) { + resultant.header = scratch; + current_line++; + } + while(!table_file.eof()) { + table_file >> hnu; + table_file >> resultant.value[hnu]; + } +} + + +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; + std::getline(table_file,scratch); + if(!isdigit(scratch[0])) { + resultant.header = scratch; + current_line++; + } + int c=0; + while(!table_file.eof()) { + //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]; + getline(table_file,scratch); + } +} + +std::string agn::format_sed_table(agn::sed_table table) { + std::stringstream output; + if (!table.header.empty()) output << table.header; + output << std::setprecision(5); + agn::table2d::iterator table_iterator; + table_iterator=table.value.begin(); + while(table_iterator != table.value.end()) { + output + << std::fixed + << table_iterator->first + << "\t" + << std::scientific + << table_iterator->second + << "\n"; + table_iterator++; + } + return output.str(); +} + +std::string agn::cloudy_interpolate_str(agn::sed_table table) { + std::stringstream output; + agn::table2d::iterator table_iterator = table.value.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; + double log_uv_bump_start = log10( energy_in_rydbergs ); + double log_lowest_value = log10(table_iterator->second + / table_iterator->first); + double log_min_energy = log10(agn::CLOUDY_EMM) + - 1; + double log_SED_density = log_lowest_value + - 2*(log_uv_bump_start + - log_min_energy); + if ( log_SED_density < 1e-36 ) log_SED_density = 1e-36; + output + << "interpolate (" + << pow(10,log_min_energy) + << " " + << log_SED_density + << ")"; + int count=0; + + while(table_iterator != table.value.end()) { + energy_in_rydbergs = table_iterator->first + / agn::RYDBERG_UNIT_EV; + double log_SED_density = log10( table_iterator->second + / table_iterator->first); + if ((count%5)==0) output << "\n" << "continue "; + else output << " "; + output + << "(" + << energy_in_rydbergs + << " " + << log_SED_density + << ")"; + count++; + table_iterator++; + } + // Trail off at slope=-2 in log(energy [rydberg]) space + while ( energy_in_rydbergs < agn::CLOUDY_EGAMRY ) { + double log_energy = log10(energy_in_rydbergs); + energy_in_rydbergs = pow(10,log_energy+1); + log_SED_density -= 2; + output + << "(" + << energy_in_rydbergs + << " " + << log_SED_density + << ")"; + } + return output.str(); +} + + +double agn::hnu_at(int i,int n) { + double relative_coord=(double)(i)/n; + double x_coord = relative_coord*CONT_WIDTH_X + CONT_MIN_X; + return pow(10,x_coord); +} + + +#endif \ No newline at end of file diff --git a/src/sed/mehdipour/generate_sed.cpp b/src/sed/mehdipour/generate_sed.cpp index f35a9b3..6028a5c 100644 --- a/src/sed/mehdipour/generate_sed.cpp +++ b/src/sed/mehdipour/generate_sed.cpp @@ -98,93 +98,3 @@ int main(int argc, char const *argv[]) return 0; } - - - - - - - - - - - -double agn::hnu_at(int i,int n) { - double relative_coord=(double)(i)/n; - double x_coord = relative_coord*CONT_WIDTH_X + CONT_MIN_X; - return pow(10,x_coord); -} - -agn::sed_table agn::sed_pow_law::histogram_table(int n){ - agn::sed_table output; - double max=0,min=1,hnu; - for(int i=0; ised(hnu); - if (output.value[hnu] > max) max = output.value[hnu]; - if (output.value[hnu] < min) min = output.value[hnu]; - } - // Add a final point at 100 KeV - hnu = 1e5; - output.value[hnu] = this->sed(hnu); - return output; -} - -double agn::sed_pow_law::sed(double hnu) { - double magnitude=0.0; - magnitude += this->eval_uv(hnu); - magnitude += this->eval_xray(hnu); - if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; - return magnitude; -} -double agn::sed_pow_law::eval_uv(double hnu) { - double bigbump_kT = _T - * agn::BOLTZMANN_CONST; - double magnitude = pow(hnu,(1+_alpha_uv)) - * exp(-(hnu)/bigbump_kT) - * exp(-(_cutoff_uv_eV/hnu)) - * _scaling_factor; - if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; - return magnitude; -} -double agn::sed_pow_law::eval_xray(double hnu) { - return _xray_coefficient - * pow(hnu/2000,1+_alpha_x) - * exp(-_cutoff_xray_eV/hnu) - * _scaling_factor; -} - -double agn::sed_pow_law::SED_at_2KeV() { - double ELe_at_2500A_no_scale = eval_uv(IN_EV_2500A) - / _scaling_factor; - double energy_ratio = 2000/IN_EV_2500A; - // Returns EL[e] at 2 KeV - return ELe_at_2500A_no_scale - * pow(energy_ratio,_alpha_ox + 1); -} - -agn::sed_pow_law::sed_pow_law ( - double T, - double alpha_ox, - double alpha_x, - double alpha_uv, - double cutoff_uv_rydberg, - double cutoff_xray_rydberg, - double log_radius_in_cm, - double scaling_factor - ): - _T(T), - _alpha_ox(alpha_ox), - _alpha_x(alpha_x), - _alpha_uv(alpha_uv), - _cutoff_uv_rydberg(cutoff_uv_rydberg), - _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(); -} \ No newline at end of file diff --git a/src/sed/mehdipour/powlaw-cutoff.pl b/src/sed/mehdipour/powlaw-cutoff.pl index 78a496c..6b01e2c 100755 --- a/src/sed/mehdipour/powlaw-cutoff.pl +++ b/src/sed/mehdipour/powlaw-cutoff.pl @@ -1,4 +1,4 @@ -#!/usr/local/bin/perl +#!/usr/bin/env use strict; use warnings; use 5.010; use utf8; use IO::Handle; diff --git a/src/sed/mehdipour/sed.hpp b/src/sed/mehdipour/sed.hpp index df3e05a..2639142 100644 --- a/src/sed/mehdipour/sed.hpp +++ b/src/sed/mehdipour/sed.hpp @@ -23,6 +23,9 @@ const double CLOUDY_MAX_EV=CLOUDY_EGAMRY*RYDBERG_UNIT_EV; const double IN_EV_2500A=12398.41929/2500; +// Pulled from cloudy 17.00, first version +const double emm = 1.001e-8f; +const double egamry = 7.354e6f; // SEDs are represented by 2d histogram tables. @@ -215,5 +218,86 @@ std::string agn::cloudy_interpolate_str(agn::sed_table table) { } +double agn::hnu_at(int i,int n) { + double relative_coord=(double)(i)/n; + double x_coord = relative_coord*CONT_WIDTH_X + CONT_MIN_X; + return pow(10,x_coord); +} + +agn::sed_table agn::sed_pow_law::histogram_table(int n){ + agn::sed_table output; + double max=0,min=1,hnu; + for(int i=0; ised(hnu); + if (output.value[hnu] > max) max = output.value[hnu]; + if (output.value[hnu] < min) min = output.value[hnu]; + } + // Add a final point at 100 KeV + hnu = 1e5; + output.value[hnu] = this->sed(hnu); + return output; +} + +double agn::sed_pow_law::sed(double hnu) { + double magnitude=0.0; + magnitude += this->eval_uv(hnu); + magnitude += this->eval_xray(hnu); + if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; + return magnitude; +} +double agn::sed_pow_law::eval_uv(double hnu) { + double bigbump_kT = _T + * agn::BOLTZMANN_CONST; + double magnitude = pow(hnu,(1+_alpha_uv)) + * exp(-(hnu)/bigbump_kT) + * exp(-(_cutoff_uv_eV/hnu)) + * _scaling_factor; + if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL; + return magnitude; +} +double agn::sed_pow_law::eval_xray(double hnu) { + return _xray_coefficient + * pow(hnu/2000,1+_alpha_x) + * exp(-_cutoff_xray_eV/hnu) + * _scaling_factor; +} + +double agn::sed_pow_law::SED_at_2KeV() { + double ELe_at_2500A_no_scale = eval_uv(IN_EV_2500A) + / _scaling_factor; + double energy_ratio = 2000/IN_EV_2500A; + // Returns EL[e] at 2 KeV + return ELe_at_2500A_no_scale + * pow(energy_ratio,_alpha_ox + 1); +} + +agn::sed_pow_law::sed_pow_law ( + double T, + double alpha_ox, + double alpha_x, + double alpha_uv, + double cutoff_uv_rydberg, + double cutoff_xray_rydberg, + double log_radius_in_cm, + double scaling_factor + ): + _T(T), + _alpha_ox(alpha_ox), + _alpha_x(alpha_x), + _alpha_uv(alpha_uv), + _cutoff_uv_rydberg(cutoff_uv_rydberg), + _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(); +} + + #endif \ No newline at end of file