industry-website/htdocs/examples/sed.hpp

#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