Otho/industry-website
1
0
Fork 0
mirror of https://asciireactor.com/otho/industry-website.git synced 2024-11-22 19:15:07 +00:00
Code Issues Packages Projects Releases Wiki Activity
eec02bd1c2
industry-website/htdocs/examples/agn.hpp

552 lines
14 KiB
C++
Raw Blame History

#ifndef agn_hpp
#define agn_hpp
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>
#include <cmath>
#include <vector>
#include <map>
#include <list>
#include <cctype>
#include <cstdlib>
#include <set>
#include <cmath>
namespace agn {
const bool debug = false;
const bool line_debug = false;
const bool verbose=true;
// General constants
const double PI=3.14159265358979323846;
const double PLANCK_CONST=4.135668e-15; // in eV * s
const double BOLTZMANN_CONST=0.00008617332385; // in eV / K
const double RYDBERG_CONST=1.0973731568539e7; // in 1 / m
const double RYDBERG_UNIT_EV=13.60569252; // in eV
const double RYDBERG_UNIT_ANGSTROM=1e10/RYDBERG_CONST; // in A
// line label header lengths based on cloudy version
// For cloudy 17, the characters reserved for the output quantities:
// Label: 0-9
// Wavelength: 10-20
// radiated energy: 21-26
// eq width: 27-36
// Cloudy 13 had: 1,13; 17,5; 22,11;
// c17: 0,18; 18,9; 27, 10
// c13 settings
//int labelstr_pos = 1;
//int labelstr_len = 13;
//int radiatedenergystr_pos = 14;
//int radiatedenergystr_len = 8;
//int eqwidthstr_pos = 22;
//int eqwidthstr_len = 9;
// c17 settings
int LABELSTR_POS = 0;
int LABELSTR_LEN = 18;
int RADIATEDENERGYSTR_POS = 18;
int RADIATEDENERGYSTR_LEN = 9;
int EQWIDTHSTR_POS = 27;
int EQWIDTHSTR_LEN = 10;
// Trims whitespace from line labels
std::string label_trim(std::string);
// Some useful containers and functions.
typedef std::map<double,double> table1d;
typedef table1d::iterator iterator1d;
std::string format_table1d(table1d);
typedef std::map<double,table1d> table2d;
typedef table2d::iterator iterator2d;
typedef std::pair<double,double> coord2d;
typedef std::list<coord2d> gridcoordlist;
// Seeks an instream to the line after the first occurrence of seek_string.
void seek_to(std::string,std::istream&);
struct cloudy_line_data {
double radiated_energy,eq_width;
int index;
bool has_duplicates;
};
typedef std::map<std::string,cloudy_line_data> cloudy_line_output;
struct cloudy_result {
std::string header, footer;
std::list<std::string> intrinsic_line_raw_text;
std::list<std::string> emergent_line_raw_text; // Not currently collected
std::list<std::string> cautions;
cloudy_line_output emergent_line_intensity;
cloudy_line_output intrinsic_line_intensity;
cloudy_line_output intrinsic_line_inward_intensity;
int iterations;
double phi,hden,colden;
cloudy_result():
header(""),
footer(""),
intrinsic_line_raw_text(),
cautions(),
iterations(0),
phi(0),
hden(0),
colden(0)
{}
};
typedef std::map<coord2d,cloudy_result> cloudy_grid;
// Easiest to read the entire grid from the file at once.
cloudy_grid read_cloudy_grid(std::ifstream&);
// Operator<< prints general info about the run result.
std::ostream& operator<< (std::ostream&, cloudy_result);
} // end namespace agn
// Definitions
std::string agn::label_trim(std::string label) {
std::stringstream newlabel;
newlabel
<< label[0]
<< label[1]
<< label[2]
<< label[3]
<< " "
<< label[agn::LABELSTR_LEN-8]
<< label[agn::LABELSTR_LEN-7]
<< label[agn::LABELSTR_LEN-6]
<< label[agn::LABELSTR_LEN-5]
<< label[agn::LABELSTR_LEN-4]
<< label[agn::LABELSTR_LEN-3]
<< label[agn::LABELSTR_LEN-2]
<< label[agn::LABELSTR_LEN-1];
for (int i=agn::LABELSTR_LEN; i<label.length(); i++)
newlabel << label[i];
return newlabel.str();
}
void agn::seek_to(std::string seek_string,std::istream& stream) {
std::string test_str="";
while (test_str.find(seek_string) == std::string::npos)
getline(stream,test_str);
}
std::string agn::format_table1d(agn::table1d table) {
std::stringstream output;
output
<< std::scientific
<< std::setprecision(5);
agn::table1d::iterator x=table.begin();
while(x!=table.end()) {
output
<< x->first
<< "\t"
<< x->second
<< "\n";
x++;
}
return output.str();
}
std::ostream& agn::operator<< (std::ostream& outstream, cloudy_result output) {
outstream.width(5);
outstream
<< std::fixed
<< std::setprecision(2)
<< output.hden
<< " x "
<< output.phi
<< " (colden "
<< std::fixed
<< std::setprecision(0)
<< output.colden
<< "): "
<< output.intrinsic_line_intensity.size()
<< " intrinsic emission lines, "
<< output.cautions.size()
<< " cautions present.";
return outstream;
}
std::ifstream& operator>> (std::ifstream& inputfile,agn::cloudy_grid& grid) {
grid = agn::read_cloudy_grid(inputfile);
return inputfile;
}
agn::cloudy_grid agn::read_cloudy_grid(std::ifstream& inputfile) {
if(agn::debug) std::cout
<< "Constructing cloudy output grid from file.\n";
inputfile.clear();
inputfile.seekg(0);
std::string inputline;
int line_num=0;
getline(inputfile,inputline);
std::string curstr = inputline.substr(0,23);
std::string seek_string = "Producing grid output.";
int pos = curstr.find(seek_string);
while(line_num < 1000) {
line_num++;
getline(inputfile,inputline);
curstr = inputline.substr(0,30);
pos = inputline.find(seek_string);
if (pos != std::string::npos) break;
}
if(agn::debug) std::cout
<< "Collecting grid meta data starting on line "
<< line_num
<< std::endl;
agn::gridcoordlist coordlist;
std::pair<double,double> xy (0,0);
// collect parameter coordinates
seek_string="**************************************************";
while(!inputfile.eof()) {
while (inputline.find(seek_string) == std::string::npos)
getline(inputfile,inputline);
while(inputfile.peek() == ' ') inputfile.get(); // Skip blank space
// Break for first grid
if (inputline.find("GRID_DELIMIT") != std::string::npos) break;
char nextchar = inputfile.peek();
//bool hascoords = nextchar != '*' && nextchar != ' ' && nextchar != '\n';
bool hascoords = isalnum(nextchar);
if (hascoords) {
//if(agn::debug) std::cout << "Has coords.\n";
while (!isdigit(inputfile.peek())) inputfile.get();
inputfile >> xy.first;
while (!isdigit(inputfile.peek())) inputfile.get();
inputfile >> xy.second;
std::pair<double,double> newentry (xy);
coordlist.push_back(newentry);
getline(inputfile,inputline);
}
else {
getline(inputfile,inputline);
while (inputline.find(seek_string) == std::string::npos)
getline(inputfile,inputline);
}
}
if(agn::debug) std::cout
<< "Reached first grid. "
<< coordlist.size()
<< " coordinate pairs found.\n";
agn::cloudy_grid grid;
agn::gridcoordlist::iterator coords = coordlist.begin();
while(grid.size() < coordlist.size()) {
bool failmode = false;
agn::cloudy_result point;
point.hden = coords->first;
point.phi = coords->second;
seek_string = "Intrinsic line intensities";
std::string fail_string = "PROBLEM DISASTER";
if(agn::debug) std::cout
<< "Grabbing header and cautions.";
std::string header="";
std::list<std::string> cautions;
while (inputline.find(seek_string) == std::string::npos) {
getline(inputfile,inputline);
if (inputline.find(fail_string) != std::string::npos) {
if(agn::debug) std::cout
<< " Found broken model.";
failmode = true;
break;
}
header.append(inputline);
header.append("\n");
if((inputline[1] == 'C' && inputline[2] == '-')
|| (inputline[1] == ' ' && inputline[2] == '!' ))
cautions.push_back(inputline);
}
if (failmode) {
grid[*coords] = point;
coords++;
//agn::seek_to("c ======================",inputfile);
//agn::seek_to("c ======================",inputfile);
continue;
}
point.header = header;
point.cautions = cautions;
std::stringstream headerstr;
headerstr << header;
std::string headerword="";
headerstr.clear();
headerstr.seekg(0);
int iterations;
while(!headerstr.eof()) {
headerstr >> headerword;
//std::cout << "looking for iterations in " << headerword << "\n";
if(headerword == "Iteration") {
headerstr >> iterations;
break;
}
}
if(agn::debug) std::cout
<< " Found " << iterations << " iterations.";
point.iterations = iterations;
headerstr.clear();
headerstr.seekg(0);
double colden;
while(!headerstr.eof()) {
headerstr >> headerword;
//std::cout << "checking "<<headerword<<"\n";
if(headerword == "column") {
headerstr >> headerword;
//std::cout << "found "<<headerword;
if(headerword == "density") {
headerstr >> colden;
break;
}
}
}
if(agn::debug) std::cout << std::setprecision(2)
<< " Found 10^"<<colden<<" colden.";
point.colden = colden;
if(agn::debug) std::cout
<< " Grabbing emission lines.";
agn::seek_to("general properties...................",inputfile);
std::list<std::string> intrinsic_line_raw_text;
getline(inputfile,inputline);
while (inputline != "") {
if(inputline.find("..") != std::string::npos) {
getline(inputfile,inputline);
continue;
}
intrinsic_line_raw_text.push_back(inputline);
getline(inputfile,inputline);
}
point.intrinsic_line_raw_text.swap(intrinsic_line_raw_text);
std::list<std::string>::iterator linetext_it=point.intrinsic_line_raw_text.begin();
std::list<std::string> duplicate_labels;
int index=0;
std::string label;
// For cloudy 17, the characters reserved for the output quantities:
// Label: 0-9
// Wavelength: 10-20
// radiated energy: 21-26
// eq width: 27-36
// Cloudy 13 had: 1,13; 17,5; 22,11;
// c17: 0,18; 18,9; 27, 10
// c13 settings
//int labelstr_pos = 1;
//int labelstr_len = 13;
//int radiatedenergystr_pos = 14;
//int radiatedenergystr_len = 8;
//int eqwidthstr_pos = 22;
//int eqwidthstr_len = 9;
// c17 settings
int labelstr_pos = 0;
int labelstr_len = 18;
int radiatedenergystr_pos = 18;
int radiatedenergystr_len = 9;
int eqwidthstr_pos = 27;
int eqwidthstr_len = 10;
while(linetext_it != point.intrinsic_line_raw_text.end()) {
// This section is also only written for the "final case" as alerted below
if ((*linetext_it).substr(labelstr_pos,4) == "Inwd") {
agn::cloudy_line_data inwddata;
inwddata.index = ++index;
inwddata.radiated_energy = atof((*linetext_it).substr(
radiatedenergystr_pos,
radiatedenergystr_len).c_str()
);
inwddata.eq_width = atof((*linetext_it).substr(
eqwidthstr_pos,
eqwidthstr_len).c_str()
);
if(line_debug) {
std::cout
<< std::setprecision(5)
<< "Inwd "
<< (*linetext_it).substr(5)
<< ": "
<< label
<< "; "
<< inwddata.radiated_energy
<< " from "
<< (*linetext_it).substr(
radiatedenergystr_pos,
radiatedenergystr_len).c_str()
<< "; "
<< inwddata.eq_width
<< " from "
<< (*linetext_it).substr(
eqwidthstr_pos,
eqwidthstr_len).c_str()
<< "\n";
}
point.intrinsic_line_inward_intensity[label] = inwddata;
linetext_it++;
continue;
}
label=(*linetext_it).substr(
labelstr_pos,
labelstr_len
);
agn::cloudy_line_data data;
data.index = ++index;
//std::stringstream values;
data.radiated_energy = atof((*linetext_it).substr(
radiatedenergystr_pos,
radiatedenergystr_len).c_str()
);
data.eq_width = atof((*linetext_it).substr(
eqwidthstr_pos,
eqwidthstr_len).c_str()
);
if(line_debug) {
std::cout
<< std::setprecision(5)
<< *linetext_it
<< ": "
<< label
<< "; "
<< data.radiated_energy
<< " from "
<< (*linetext_it).substr(
radiatedenergystr_pos,
radiatedenergystr_len).c_str()
<< "; "
<< data.eq_width
<< " from "
<< (*linetext_it).substr(
eqwidthstr_pos,
eqwidthstr_len).c_str()
<< "\n";
}
if(point.intrinsic_line_intensity.count(label) == 0) {
data.has_duplicates = false;
point.intrinsic_line_intensity[label] = data;
}
else {
duplicate_labels.push_back(label);
std::stringstream zerolabel(label);
zerolabel << " j=0";
point.intrinsic_line_intensity[zerolabel.str()]
= point.intrinsic_line_intensity[label];
point.intrinsic_line_intensity[zerolabel.str()].has_duplicates = true;
point.intrinsic_line_intensity.erase(label);
int j=1;
data.has_duplicates = true;
while(true) {
std::stringstream newlabel;
newlabel << label;
newlabel << " j=" << j;
if(point.intrinsic_line_intensity.count(newlabel.str()) != 0) {
j++;
// RED ALERT
// To ease work in summer 2017, this line sets the
// final instance of the emission line value
// to the original label, this was to get around
// the 3645.00A and 3647.00A being reported with
// different values in c17 version 1
point.intrinsic_line_intensity[label] = data;
continue;
}
else {
point.intrinsic_line_intensity[newlabel.str()] = data;
break;
}
}
}
linetext_it++;
}
if(agn::debug) std::cout << " Duplicates found: " << duplicate_labels.size();
if(agn::line_debug) {
std::list<std::string>::iterator dup_it = duplicate_labels.begin();
while (dup_it != duplicate_labels.end() ) {
std::cout << std::endl;
std::cout << *dup_it;
dup_it++;
}
std::cout << " Grabbing footer.";
}
while (inputline == "")
getline(inputfile,inputline);
std::string footer="";
seek_string="GRID_DELIMIT";
while (inputline.find(seek_string) == std::string::npos){
footer.append(inputline);
footer.append("\n");
getline(inputfile,inputline);
}
if(agn::debug) std::cout
<< "\nAdding point to grid ("
<< grid.size()
<< " of "
<< coordlist.size()
<< ") at "
<< coords->first
<< ": "
<< point
<< std::endl;
grid[*coords] = point;
coords++;
}
if(agn::debug) std::cout
<< "Grid captured. "
<< grid.size()
<< " points compiled.\n";
return grid;
}
#endif
Reference in New Issue View Git Blame Copy Permalink