industry-website/htdocs/examples/spectral_lines.hpp

#ifndef spectral_lines_hpp
#define spectral_lines_hpp

#include "agn.hpp"
//#include "png++-0.2.5/gray_pixel.hpp"
//#include "png++-0.2.5/image.hpp"


namespace agn {

const double EQWIDTH_MIN_VAL = 1E-7;
const double EQWIDTH_MIN_VAL_LOG = -7;
const double SIGMA_THRESHOLD_MULTIPLIER = 1.0;
const double RATIO_THRESHOLD_MULTIPLIER = .06;

typedef std::list<std::string> line_list;
line_list read_line_list(std::ifstream&);

// Trims whitespace from line labels
std::string label_trim(std::string);

// Emission line flux contours are represented by a singly-nested sorted map.
struct eqwidth_table {
	std::string header[2];
	table2d value;
};
bool is_zero(const eqwidth_table&);
eqwidth_table read_eqwidth_table(std::ifstream& table_file);
std::string format_eqwidth_table(eqwidth_table table);
std::string format_eqwidth_table_slice(eqwidth_table table,iterator2d x);
gridcoordlist find_outliers(eqwidth_table table);
gridcoordlist known_outliers();
eqwidth_table smooth(eqwidth_table&, gridcoordlist coordlist);
std::list<eqwidth_table> compile_eqwidth_tables(cloudy_grid,line_list,double = 1.0);
// Operator<< prints table formatted for a "fort file".
std::ostream& operator<< (std::ostream&,eqwidth_table);

} // end namespace agn


































// Definitions

bool agn::is_zero(const eqwidth_table& table_to_check) {
    agn::table2d::const_iterator x_it = table_to_check.value.begin();
    while (x_it != table_to_check.value.end()) {
        agn::table1d::const_iterator y_it= x_it->second.begin();
        while (y_it != x_it->second.end()) {
            if (y_it->second != EQWIDTH_MIN_VAL)
                return false;
                y_it++;
        }
        x_it++;
    }
    return true;
}


agn::eqwidth_table agn::read_eqwidth_table(std::ifstream& table_file) {
    agn::eqwidth_table resultant;
    double hden, phi, flux;
    int number_of_lines = 0;
    std::string null;
    while (std::getline(table_file, null))
        ++number_of_lines;
    int table_entries = number_of_lines - 2;
    table_file.clear();
    table_file.seekg(0);
    getline(table_file,resultant.header[0]);
    getline(table_file,resultant.header[1]);
    for (int i=0; i<table_entries; i++) {
        table_file >> hden >> phi >> flux;
        resultant.value[hden][phi]=flux;
    }
    return resultant;
}

std::string agn::format_eqwidth_table(agn::eqwidth_table table) {
    std::stringstream output;
    output
        << table.header[0]
        << "\n"
        << table.header[1]
        << "\n";
    agn::iterator2d hden_it = table.value.begin();
    while(hden_it != table.value.end()) {
        agn::iterator1d phi_it = table.value[hden_it->first].begin();
        while(phi_it != table.value[hden_it->first].end()) {
            output
                << std::fixed
                << std::setprecision(3)
                << hden_it->first
                << "  "
                << phi_it->first
                << "  "
                << std::scientific
                << std::setprecision(4)
                << phi_it->second
                << "\n";
            phi_it++;
        }
        hden_it++;
    }
    return output.str();
}

std::string agn::format_eqwidth_table_slice(agn::eqwidth_table table,agn::iterator2d fixed) {
    std::stringstream output("");
//    output
//        << table.header[0]
//        << " - constant="
//        << fixed->first
//        << "\n";
    iterator1d x=fixed->second.begin();
    while(x!=fixed->second.end()) {
        output
            << std::fixed
            << std::setprecision(3)
            << x->first
            << "\t"
            << std::scientific
            << std::setprecision(4)
            << x->second
            << "\n";
        x++;
    }
    return output.str();
}


std::ostream& agn::operator<< ( std::ostream& outstream,
                                agn::eqwidth_table table) {
    std::string header0_trimmed = agn::label_trim(table.header[0]);
	outstream
		<< header0_trimmed
		<< std::endl
		<< table.header[1]
		<< std::endl;
	agn::table2d::iterator x_it = table.value.begin();
	while (x_it != table.value.end()) {
		agn::table1d::iterator y_it= x_it->second.begin();
		while (y_it != x_it->second.end()) {
			outstream.width(6);
			outstream
				<< std::fixed
				<< std::setprecision(3)
				<< x_it->first
				<< "\t"
				<< y_it->first
				<< "\t"
				<< std::scientific
				<< std::setprecision(4)
				<< y_it->second
				<< "\n";
				y_it++;
		}
		x_it++;
	}
	return outstream;
}


// Returns a table as the difference between two tables
//agn::eqwidth_table main(int argc, char const *argv[])
//{  // Copied from fortfile subtractor
//    std::ifstream table1;
//    table1.open(argv[1]);
//    std::ifstream table2;
//    table2.open(argv[2]);
//
//    std::string line;
//
//    int table1_numlines = 0;
//    while (std::getline(table1, line))
//        ++table1_numlines;
//    table1.clear();
//    table1.seekg(0);
//
//    int table_entries = table1_numlines - 2;
//
////    int table2_numlines = 0;
////    while (std::getline(table1, line))
////        ++table2_numlines;
////
////    int table2_entries = table2_numlines - 2;
////    table2.clear();
////    table2.seekg(0);
//
//    double hden1[table_entries];
//    double phi1[table_entries];
//    double magnitude1[table_entries];
//    double resultant_magnitude[table_entries];
//
//    std::string first_line, second_line, scratch;
//    getline(table1,first_line);
//    getline(table1,second_line);
//    getline(table2,scratch);
//    getline(table2,scratch);
//
//    first_line[0] = argv[3][0];
//    first_line[1] = argv[3][1];
//    first_line[2] = argv[3][2];
//    first_line[3] = argv[3][3];
//
//    for (int i=0; i<table_entries; i++) {
//        table1 >> hden1[i] >> phi1[i] >> magnitude1 [i];
//    }
//
//
//    int i=0;
//    double hden2, phi2, magnitude2;
//    for (int i=0; i<table_entries; i++) {
//        if ( table2.eof() ) {
//            resultant_magnitude[i]=magnitude1[i];
//        }
//        else {
//        table2 >> hden2 >> phi2 >> magnitude2;
//        resultant_magnitude[i] = magnitude1[i] - magnitude2;
//        }
//        if ( resultant_magnitude[i] < 1e-10 )
//            resultant_magnitude[i]=1e-10;
//    }
//
//
//
//    std::cout
//            << first_line
//            << std::endl
//            //<< "HDen   Phi(H)  ScaInten"
//            << second_line
//            << std::endl;
//
//    for (int i=0; i<table_entries; i++) {
//        std::cout
//            << std::fixed
//            << std::setprecision(3)
//            << hden1[i]
//            << "  "
//            << phi1[i]
//            << "  "
//            << std::scientific
//            << std::setprecision(4)
//            << resultant_magnitude[i]
//            << std::endl;
//    }
//
//    table1.close();
//    table2.close();
//    return 0;
//}


agn::line_list agn::read_line_list(std::ifstream& inputfile) {
	line_list list;
    line_list nfnulist;
    line_list inwtlist;
    line_list inwclist;
	inputfile.clear();
	inputfile.seekg(0);
	std::string line,label;
	while(!inputfile.eof()) {
		getline(inputfile,line);
		if (line[0] == '#' || line.size() < 2) continue;
		if (line[0] == ' ' && line[1] == ' ') continue;
		if (line[0] == ' ')
			label = line.substr(1,agn::LABELSTR_LEN);
		else
			label = line.substr(0,agn::LABELSTR_LEN);
        // above is only good for c17?
		if(agn::debug) std::cout
			<< "Adding label "
			<< label
			<< " to capture list.\n";
        std::string applabel = label.substr(4);
        if (label.substr(0,4).compare(std::string("nFnu")) == 0) {
            nfnulist.push_back(std::string("nFnu").append(applabel));
            //list.push_back(std::string("nInu").append(applabel));
            inwtlist.push_back(std::string("InwT").append(applabel));
            inwclist.push_back(std::string("InwC").append(applabel));
            continue;
        }
		list.push_back(label);
        //if (label.substr(0,4) != "Inci" && applabel != "            0 ")
        //    list.push_back(std::string("Inwd").append(applabel));
	}
    list.insert(list.end(),nfnulist.begin(),nfnulist.end());
    list.insert(list.end(),inwtlist.begin(),inwtlist.end());
    list.insert(list.end(),inwclist.begin(),inwclist.end());
	return list;
}


std::list<agn::eqwidth_table> agn::compile_eqwidth_tables(agn::cloudy_grid grid,agn::line_list line_labels,double scale_factor) {
	std::list<eqwidth_table> table_list_eq_width,table_list_radiated_energy;
	agn::line_list::iterator line_label_it = line_labels.begin();
    // Iterate over labels, creating table for each:
	while (line_label_it != line_labels.end()) {
		double x,y;
		agn::cloudy_line_data data;
		agn::eqwidth_table new_table;
		std::string label = *line_label_it;
        if(agn::debug) std::cout
            << "Processing label "
            << label
            << std::endl;
        std::stringstream header0;
        header0
            << label;
        new_table.header[0] = header0.str();
        new_table.header[1] = "Hden    Phi(H)  F(line)/Inci(1215) ";
        agn::cloudy_grid::iterator result_it = grid.begin();
        while (result_it != grid.end()) {
            x = result_it->first.first;
            y = result_it->first.second;
            // if no data exist under this label, initialize some:
			if ( result_it->second.intrinsic_line_intensity.count(label) == 0 ) {
				data.radiated_energy=-35.0;
				data.eq_width=0.0;
			}
			else {
				data = result_it->second.intrinsic_line_intensity[label];
				data.eq_width /= scale_factor;
				if(agn::line_debug) std::cout
					<< "Added  "
					<< std::setprecision(2)
					<< std::fixed
					<< x
					<< ", "
					<< y
					<< " -- "
					<< std::scientific
					<< data.radiated_energy
					<< ", "
					<< data.eq_width
					<< "\n";
				if ( data.has_duplicates ) {
					//	This needs to be expanded to catch all duplicate values j=1,2,etc.
					// For now, just blindly use the value the program finds...
                    // the functionality is written in the agn.hpp sections, but is not
                    // implemented here, and as of 2017, the agn.hpp function writes the
                    // last-encountered instance of the emission line value to the label
                    // probably needs to be fixed for the inwd part too
				}
			}
			if ( data.radiated_energy < EQWIDTH_MIN_VAL_LOG || data.eq_width < EQWIDTH_MIN_VAL) {
				data.radiated_energy = EQWIDTH_MIN_VAL_LOG;
				data.eq_width = EQWIDTH_MIN_VAL;
			}
			new_table.value[x][y] = data.eq_width;
			result_it++;
		}
		table_list_eq_width.push_back(new_table);
        // Add inward table for emission lines, and only works for c17
        if (label.substr(agn::LABELSTR_POS,4) != "Inci" &&
            label.substr(agn::LABELSTR_POS,4) != "FeKa" &&
            //label.substr(agn::LABELSTR_POS,4) != "Blnd" &&
            label.substr(agn::LABELSTR_POS,3) != "nFn" &&
            label.substr(agn::LABELSTR_POS,3) != "Inw" &&
            label.substr(LABELSTR_LEN - 14) != "            0 ") {
            // These lines are known to not have inward fractions, so
            // should be skipped.
            if (label == "Blnd      1657.00A" ||
                label == "Blnd      2326.00A" ||
                label == "Blnd      2798.00A" ||
                label == "Blnd      1335.00A" ||
                label == "Blnd      977.000A" ||
                label == "Blnd      1909.00A" ||
                label == "Blnd      1176.00A" ||
                label == "Blnd      2141.00A" ||
                label == "Blnd      1085.00A" ||
                label == "Blnd      1750.00A" ||
                label == "Blnd      990.000A" ||
                label == "Blnd      1486.00A" ||
                label == "Blnd      765.000A" ||
                label == "Blnd      1304.00A" ||
                label == "Blnd      8446.00A" ||
                label == "Blnd      1666.00A" ||
                label == "Blnd      835.000A" ||
                label == "Blnd      1402.00A" ||
                label == "Blnd      630.000A" ||
                label == "Blnd      1218.00A" ||
                label == "Blnd      774.000A" ||
                label == "Blnd      5892.00A" ||
                label == "Blnd      615.000A" ||
                label == "Blnd      2665.00A" ||
                label == "Blnd      1860.00A" ||
                label == "Blnd      2335.00A" ||
                label == "Blnd      1888.00A" ||
                label == "Blnd      1397.00A" ||
                label == "Blnd      1256.00A" ||
                label == "Blnd      1720.00A" ||
                label == "Blnd      1198.00A" ||
                label == "Blnd      1406.00A" ||
                label == "Blnd      1199.00A" ||
                label == "Blnd      8579.00A" ||
                label == "Blnd      3933.00A"
                ) {
                line_label_it++;
                continue;
            }

            if(agn::debug) std::cout
                << "Processing Inwd for label "
                << label
                << std::endl;
            std::stringstream inwdheader0;
            agn::eqwidth_table inwd_table;
            inwdheader0                
                << "Inwd "
                << label.substr(5);
            inwd_table.header[0] = inwdheader0.str();
            inwd_table.header[1] = "Hden    Phi(H)  F(line)/Inci(1215)";
            result_it = grid.begin();
            agn::cloudy_line_data inwddata;
            while (result_it != grid.end()) {
                x = result_it->first.first;
                y = result_it->first.second;
                // if no data exist under this label, initialize some:
                if ( result_it->second.intrinsic_line_inward_intensity.count(label) == 0 ) {
                    inwddata.radiated_energy=-35.0;
                    inwddata.eq_width=0.0;
                }
                else {
                    inwddata = result_it->second.intrinsic_line_inward_intensity[label];
                    inwddata.eq_width /= scale_factor;
                }
                if(agn::line_debug) std::cout
                    << "Added  "
                    << std::setprecision(2)
                    << std::fixed
                    << x
                    << ", "
                    << y
                    << " -- "
                    << std::scientific
                    << inwddata.radiated_energy
                    << ", "
                    << inwddata.eq_width
                    << "\n";
                if ( inwddata.has_duplicates ) {
                    //  This needs to be expanded to catch all duplicate values j=1,2,etc.
                    // For now, just blindly use the value the program finds...
                    // the functionality is written in the agn.hpp sections, but is not
                    // implemented here, and as of 2017, the agn.hpp function writes the
                    // last-encountered instance of the emission line value to the label
                    // probably needs to be fixed for the inwd part too
                }
                if ( inwddata.radiated_energy < EQWIDTH_MIN_VAL_LOG || 
                    inwddata.eq_width < EQWIDTH_MIN_VAL) {
                    inwddata.radiated_energy = EQWIDTH_MIN_VAL_LOG;
                    inwddata.eq_width = EQWIDTH_MIN_VAL;
                }
                inwd_table.value[x][y] = inwddata.eq_width;
                result_it++;
            }
            table_list_eq_width.push_back(inwd_table);
        }
		line_label_it++;
	}
	return table_list_eq_width;
}


#endif