Otho/industry-website
1
0
Fork 0
mirror of https://asciireactor.com/otho/industry-website.git synced 2024-11-22 19:45:05 +00:00
Code Issues Packages Projects Releases Wiki Activity
78a0c12fb9
industry-website/htdocs/examples/sed.hpp

518 lines
15 KiB
C++
Raw Blame History

#ifndef sed_hpp
#define sed_hpp
#include "agn.hpp"
#include "spline.h"
namespace agn {
// 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;
// Continuum domain, step size constant in log space
const double CONT_MIN_ENERGY=CLOUDY_MIN_EV; // eV
const double CONT_MAX_ENERGY=CLOUDY_MAX_EV; // 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;
// SEDs are represented by 2d histogram tables.
struct sed_table {
std::string header;
table1d table;
};
// To account for the four main powerlaws in a typical
// AGN SED.
// Hardcoded infrared and gamma ray power laws and cutoffs.
const double IR_POWER = 3;
const double GAMMA_POWER = -5;
const double RADIO_CUTOFF = 1e-4; // IN KEV
const double GAMMA_CUTOFF = 1e4; // IN KEV
struct powerlaw_bounds {
double ir_min;
double ir_max;
double uv_min;
double uv_max;
double xray_min;
double xray_max;
double gamma_min;
double gamma_max;
};
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(exp(log(x0.second)-(_power*log(x0.first))))
{}
powerlaw(coord2d x0,double power):
_power(power),
_normal(exp(log(x0.second)-(_power*log(x0.first))))
{}
double eval(double hnu) {
return
_normal
* pow(hnu,_power)
* exp(-(hnu)/GAMMA_CUTOFF)
* exp(-(RADIO_CUTOFF/hnu))
;
}
};
class sed {
public:
// 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 eval(double hnu) {};
sed() {};
};
class sed_powerlaw_spline : public sed {
private:
Spline<double,double> _output_model;
powerlaw _ir_powerlaw;
powerlaw _uv_powerlaw;
powerlaw _xray_powerlaw;
powerlaw _gamma_powerlaw;
powerlaw_bounds _bounds;
// These parameters might still be useful for rolling off various quantities, but aren't used in the strict-spline case.
// 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;
public:
double eval(double hnu);
sed_powerlaw_spline(agn::sed_table& samples,
agn::sed_table& powerlaw_coords);
powerlaw * getpowerlaw(double hnu);
};
class sed_pow_law : public sed {
public:
double eval(double hnu);
// 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
);
};
// Returns coord in eV for given relative coord.
double logspace_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
agn::sed_powerlaw_spline::sed_powerlaw_spline(
agn::sed_table& samples,
agn::sed_table& powerlaw_coords
)
{
// coordinate vectors will be used to construct spline sed model
std::vector<double> x0;
std::vector<double> x1;
// powerlaws are evaluated across four regions of the sed, first
// we construct the powerlaws, here, and locate them
iterator1d table_it = powerlaw_coords.table.begin();
double ir_power = agn::IR_POWER;
double gamma_power = agn::GAMMA_POWER;
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;
_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);
_bounds.ir_min = agn::CLOUDY_MIN_EV;
_bounds.ir_max = ir_high_point.first;
_bounds.uv_min = uv_low_point.first;
_bounds.uv_max = uv_high_point.first;
_bounds.xray_min = xray_low_point.first;
_bounds.xray_max = xray_high_point.first;
_bounds.gamma_min = gamma_low_point.first;
_bounds.gamma_max = agn::CLOUDY_MAX_EV;
if(agn::debug) {
std::cout << "[Constructor] Powerlaw Boundaries: \n";
std::cout << _bounds.ir_min << std::endl;
std::cout << _bounds.ir_max << std::endl;
std::cout << _bounds.uv_min << std::endl;
std::cout << _bounds.uv_max << std::endl;
std::cout << _bounds.xray_min << std::endl;
std::cout << _bounds.xray_max << std::endl;
std::cout << _bounds.gamma_min << std::endl;
std::cout << _bounds.gamma_max << std::endl;
}
// here we inject the powerlaws into the samples
int segments=100;
double hnu = 0;
double value = 0;
for (int i=0; i<=segments; i++) {
hnu =
_bounds.ir_min +
(i/(double)segments)*(_bounds.ir_max - _bounds.ir_min);
value = _ir_powerlaw.eval(hnu);
coord2d point = coord2d(hnu,value);
samples.table.insert(samples.table.end(),point);
}
for (int i=0; i<=segments; i++) {
hnu =
_bounds.uv_min +
(i/(double)segments)*(_bounds.uv_max - _bounds.uv_min);
value = _uv_powerlaw.eval(hnu);
coord2d point = coord2d(hnu,value);
samples.table.insert(samples.table.end(),point);
}
for (int i=0; i<=segments; i++) {
hnu =
_bounds.xray_min +
(i/(double)segments)*(_bounds.xray_max - _bounds.xray_min);
value = _xray_powerlaw.eval(hnu);
coord2d point = coord2d(hnu,value);
samples.table.insert(samples.table.end(),point);
}
for (int i=0; i<=segments; i++) {
hnu =
_bounds.gamma_min +
(i/(double)segments)*(_bounds.gamma_max - _bounds.gamma_min);
value = _gamma_powerlaw.eval(hnu);
coord2d point = coord2d(hnu,value);
samples.table.insert(samples.table.end(),point);
}
if(agn::debug) {
std::cout
<< "[Constructor] Samples after evaluating powerlaws: \n";
std::cout << agn::format_sed_table(samples);
}
// 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<double,double> newspline(x0,x1);
_output_model = 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 table;
double max=0,min=1,hnu;
for(int i=0; i<n; i++) {
// evenly space coordinates in log space and save values
hnu = logspace_hnu_at(i,n);
table.table[hnu] = this->eval(hnu);
// Just collects min and max
if (table.table[hnu] > max) max = table.table[hnu];
if (table.table[hnu] < min) min = table.table[hnu];
}
return table;
}
// sed_powerlaw_spline evaluation
double agn::sed_powerlaw_spline::eval(double hnu) {
double magnitude=0.0;
agn::powerlaw * here = this->getpowerlaw(hnu);
if (here == NULL)
magnitude += this->_output_model[hnu];
else
magnitude += here->eval(hnu);
if (magnitude < agn::CONT_MIN_VAL) return agn::CONT_MIN_VAL;
return magnitude;
}
agn::powerlaw * agn::sed_powerlaw_spline::getpowerlaw(double hnu) {
if (hnu <= _bounds.gamma_max && hnu >= _bounds.gamma_min )
return &_gamma_powerlaw;
if (hnu <= _bounds.uv_max && hnu >= _bounds.uv_min )
return &_uv_powerlaw;
if (hnu <= _bounds.xray_max && hnu >= _bounds.xray_min )
return &_xray_powerlaw;
if (hnu <= _bounds.ir_max && hnu >= _bounds.ir_min )
return &_ir_powerlaw;
return NULL;
}
// sed_pow_law evaluations
double agn::sed_pow_law::eval(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
agn::sed_table agn::read_sed_table(std::ifstream& table_file) {
sed_table resultant;
std::string line;
double hnu;
if(!isdigit(table_file.peek()) && table_file.peek() != '#') {
std::getline(table_file,resultant.header);
}
while(!table_file.eof()) {
std::getline(table_file,line);
if (line[0] == '#') continue;
std::stringstream scratch(line);
scratch >> hnu;
scratch >> resultant.table[hnu];
}
return resultant;
}
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,table;
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.table[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::table1d::iterator table_iterator;
table_iterator=table.table.begin();
while(table_iterator != table.table.end()) {
output
<< std::fixed
<< std::scientific
<< 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::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;
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.table.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::logspace_hnu_at(int i,int n) {
double relative_coord_logspace=(double)(i)/n;
double abs_coord_logspace = relative_coord_logspace*CONT_WIDTH_LOGX + CONT_MIN_LOGX;
return pow(10,abs_coord_logspace);
}
#endif
Reference in New Issue View Git Blame Copy Permalink