spline sed generation

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<double,double> 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;
+}
+

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<double> x;
+    std::vector<double> 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; i<n; i++) {
+        hnu = hnu_at(i,n);
+        output.value[hnu] = this->sed(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:
+// '<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
+
+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

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; i<n; i++) {
-        hnu = hnu_at(i,n);
-        output.value[hnu] = this->sed(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();
-}

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;

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; i<n; i++) {
+        hnu = hnu_at(i,n);
+        output.value[hnu] = this->sed(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