#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 _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: // '\t\n' for each energy-flux pair std::string format_sed_table(sed_table table); // Read continuum from file with '\t\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 x0; std::vector 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 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; ieval(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