sed program starting to work

sed/mehdipour2013/mehdipour2013_powlaw_coords.tab

@@ -1,5 +1,6 @@
+12.00       -9
 0.0004  3.2e10
-0.006   4,9e13
+0.006   4.9e13
 .4      2.2e13
 2       1.4e13
 100     7.7e13

sed/mehdipour2013/mehdipour2013_samples.tab

@@ -1,4 +1,6 @@
+0.0013   9e12
 0.002   2e13
 0.005   5e13
 1.3     1.3e13
-240     8.2e13
+240     8.2e13
+300     8e13les

src/sed.hpp

@@ -35,6 +35,19 @@ struct sed_table {
 	table1d table;
 };
 
+// To account for the four main powerlaws in a typical
+// AGN SED.
+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
@@ -44,13 +57,13 @@ 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))))
+        _normal(exp(log(x0.second)-(_power*log(x0.first))))
         {}
-    powerlaw(coord2d x0,double slope):
-        _power(slope),
-        _normal((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 0; }
+    double eval(double hnu) { return _normal*pow(hnu,_power); }
 };
 
 class sed {
@@ -60,7 +73,7 @@ public:
     sed_table histogram_table(int n);
 
     // Argument is photon energy in eV
-    virtual double table(double hnu) {};
+    virtual double eval(double hnu) {};
 
     sed() {};
 };
@@ -72,6 +85,7 @@ private:
     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.
 
@@ -84,14 +98,15 @@ private:
     double _xray_coefficient;
 
 public:
-    double table(double hnu);
+    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 table(double hnu);
+    double eval(double hnu);
     // Argument is photon energy in eV
     double eval_uv(double hnu);
     double eval_xray(double hnu);
@@ -133,7 +148,7 @@ public:
 };
 
 //  Returns coord in eV for given relative coord.
-double hnu_at(int i,int n);
+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
@@ -168,37 +183,79 @@ agn::sed_powerlaw_spline::sed_powerlaw_spline(
     // 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 = 3;
+    double ir_power = (*table_it).first;
+    double gamma_power = (*table_it).second;
+    std::cout << ir_power << ", " << gamma_power;
+    table_it++;
     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);
 
-    double ir_lowerbound = agn::CLOUDY_MIN_EV;
-    double ir_upperbound = ir_high_point.first;
-    double uv_lowerbound = uv_low_point.first;
-    double uv_upperbound = uv_high_point.first;
-    double xray_lowerbound = xray_low_point.first;
-    double xray_upperbound = xray_high_point.first;
-    double gamma_lowerbound = gamma_low_point.first;
-    double gamma_upperbound = agn::CLOUDY_MAX_EV;
+    _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=10;
+    int segments=100;
+    double hnu = 0;
+    double value = 0;
     for (int i=0; i<=segments; i++) {
-        double hnu =
-            ir_lowerbound + 
-            (i/segments)*(ir_upperbound-ir_lowerbound);
-        double value = _ir_powerlaw.eval(hnu);
+        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(point);
+        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) {
@@ -249,30 +306,44 @@ agn::sed_pow_law::sed_pow_law (
 
 // writes log-space histogram with n data
 agn::sed_table agn::sed::histogram_table(int n){
-    agn::sed_table output;
+    agn::sed_table table;
     double max=0,min=1,hnu;
     for(int i=0; i<n; i++) {
-        hnu = hnu_at(i,n);
-        output.table[hnu] = this->table(hnu);
-        if (output.table[hnu] > max) max = output.table[hnu];
-        if (output.table[hnu] < min) min = output.table[hnu];
+        // 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];
     }
-    // Add a final point at 100 KeV
-    hnu = 1e5;
-    output.table[hnu] = this->table(hnu);
-    return output;
+    return table;
 }
 
 // sed_powerlaw_spline evaluation
-double agn::sed_powerlaw_spline::table(double hnu) {
+double agn::sed_powerlaw_spline::eval(double hnu) {
     double magnitude=0.0;
-    magnitude += this->_output_model[hnu];
+    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::table(double hnu) {
+double agn::sed_pow_law::eval(double hnu) {
     double magnitude=0.0;
     magnitude += this->eval_uv(hnu);
     magnitude += this->eval_xray(hnu);
@@ -315,14 +386,31 @@ double agn::sed_pow_law::SED_at_2KeV() {
 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++;
     }
+    else
+        table_file.seekg(0);
     while(!table_file.eof()) {
+        char next;
+        next = table_file.peek();
+        if (next == '\n'){
+                    table_file.get();
+                    next = table_file.peek();
+        }
+        if (next == '#') {
+            char commented;
+            while (commented != '\n') {
+                if(table_file.eof() || table_file.fail())
+                    break;
+                table_file.get(commented);
+            }
+            if(table_file.eof() || table_file.fail())
+                break;
+            next = table_file.peek();
+        }
         table_file >> hnu;
         table_file >> resultant.table[hnu];
     }
@@ -425,10 +513,10 @@ 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_LOGX + CONT_MIN_LOGX;
-    return pow(10,x_coord);
+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);
 }
 
 

src/table_powerlaw_spline_sed.cpp

@@ -28,7 +28,7 @@ int main(int argc, char const *argv[])
                             std::ofstream::out
                             );
     std::ifstream powerlaw_table(
-                            sample_filename,
+                            powerlaw_filename,
                             std::ofstream::out
                             );
     std::ofstream output_table(
@@ -45,11 +45,15 @@ int main(int argc, char const *argv[])
 
     // Read in sampling table and construct a spline model.
     samples = agn::read_sed_table(sample_table);
-    powerlaw_coords = agn::read_sed_table(powerlaw_table);
 
     if(agn::debug) debug_file 
             << "Read samples:\n"
             << format_sed_table(samples);
+    powerlaw_coords = agn::read_sed_table(powerlaw_table);
+
+    if(agn::debug) debug_file 
+            << "Read power coords:\n"
+            << format_sed_table(powerlaw_coords);
 
     agn::sed_powerlaw_spline agnsource(samples,powerlaw_coords);
 
@@ -57,7 +61,7 @@ int main(int argc, char const *argv[])
 
 
     if(agn::verbose) std::cout
-        << "Evaluating relative spectral intensity for "
+        << "Evaluating spectral intensity for "
         << n
         << " photon energy bins.\n";