Science/cloudy-agn
1
0
Fork 0
mirror of https://asciireactor.com/otho/cloudy-agn.git synced 2024-12-05 02:25:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

update

Browse Source
This commit is contained in:
caes 2017-08-09 02:57:14 -04:00
parent 1cd013467f
commit 45bc4d7030
3 changed files with 82 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/generate_powlaw_sed.cpp
View File

@ -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[])
{

138
src/sed.hpp
View File

@ -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<double,double> _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<double,double> _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<double> x;
std::vector<double> 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<double> x0;
std::vector<double> 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<double,double> newspline(x,y);
_spline = newspline;
Spline<double,double> newspline(x0,x1);
_output_model = newspline;
}
agn::sed_pow_law::sed_pow_law (
@ -208,7 +214,7 @@ 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);
@ -225,26 +231,26 @@ agn::sed_table agn::sed::histogram_table(int n){
double max=0,min=1,hnu;
for(int i=0; i<n; i++) {
hnu = hnu_at(i,n);
output.value[hnu] = this->value(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

6
src/table_powerlaw_spline_sed.cpp
View File

@ -2,6 +2,9 @@
#include "sed.hpp"
// Syntax: table_powerlaw_spline <samples table> <powerlaw coordinates> <output table>
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<double,double> represents the table.
// agn sed_powerlaw_spline class has a function for this.
int n = 1000;
agn::sed_table SED;
agn::sed_table samples;