359 lines
12 KiB
C++
359 lines
12 KiB
C++
#include "fourdst/composition/io/standard_compositions.h"
|
|
#include "fourdst/composition/io/StandardAbundancesBinary.h"
|
|
|
|
#include "fourdst/composition/composition.h"
|
|
#include "fourdst/atomic/atomicSpecies.h"
|
|
#include "fourdst/atomic/species.h"
|
|
#include "fourdst/composition/utils.h"
|
|
|
|
#include <string>
|
|
#include <vector>
|
|
#include <algorithm>
|
|
#include <fstream>
|
|
#include <numeric>
|
|
#include <stdexcept>
|
|
#include <unordered_map>
|
|
#include <print>
|
|
#include <ranges>
|
|
#include <cctype>
|
|
|
|
namespace fourdst:: composition::io {
|
|
namespace {
|
|
inline void ltrim(std::string &s) {
|
|
s.erase(
|
|
s.begin(),
|
|
std::ranges::find_if(s,
|
|
[](const unsigned char ch) {
|
|
return !std::isspace(ch);
|
|
})
|
|
);
|
|
}
|
|
|
|
inline void rtrim(std::string &s) {
|
|
s.erase(
|
|
std::find_if(
|
|
s.rbegin(),
|
|
s.rend(),
|
|
[](const unsigned char ch) {
|
|
return !std::isspace(ch);
|
|
}).base(),
|
|
s.end()
|
|
);
|
|
}
|
|
|
|
inline void trim(std::string &s) {
|
|
ltrim(s);
|
|
rtrim(s);
|
|
}
|
|
|
|
|
|
}
|
|
|
|
bool to_bool(std::string s) {
|
|
std::transform(s.begin(), s.end(), s.begin(),
|
|
[](unsigned char c){ return std::tolower(c); });
|
|
|
|
return s == "true";
|
|
}
|
|
|
|
CompositionData ChemicalFileParser::parse_compositon_data(const std::vector<char>& data,const std::string& scheme) const {
|
|
|
|
// get file and metal_fraction_scheme
|
|
// Load the file
|
|
// find the metal_fraction_scheme
|
|
// return abundances
|
|
|
|
// LOG_TRACE_L1(m_logger, "Parsing chemical abundance for: {}", scheme);
|
|
|
|
bool debug = false;
|
|
if (debug){
|
|
std::println("Parsing chemical abundance for: {}", scheme);
|
|
}
|
|
|
|
std::istringstream stream(std::string(data.begin(), data.end()));
|
|
|
|
// add error message if something goes wrong
|
|
|
|
std::string line;
|
|
int start_line = 0;
|
|
int i = 0;
|
|
|
|
CompositionData comp;
|
|
|
|
while (std::getline(stream, line)) {
|
|
|
|
// find where the end of the scheme block is
|
|
auto end_pos = std::ranges::search(line,std::format("END {}", scheme));
|
|
|
|
// exit if have reached the end of block
|
|
if (!end_pos.empty()) {
|
|
break;
|
|
}
|
|
|
|
if (start_line>0){
|
|
const size_t colon_pos = line.find(':');
|
|
line = line.substr(colon_pos+1);
|
|
|
|
line.erase(std::remove_if(line.begin(), line.end(),
|
|
[](char c){ return c == '[' || c == ']'; }),
|
|
line.end());
|
|
|
|
trim(line);
|
|
std::string item;
|
|
std::stringstream ss(line);
|
|
double val;
|
|
switch(i-start_line){
|
|
case 1:
|
|
comp.comment_str = line;
|
|
break;
|
|
case 2:
|
|
comp.he_abundance = std::pow(10.0,std::stod(line));
|
|
break;
|
|
case 3:
|
|
comp.requires_atomic_weight = to_bool(line);
|
|
break;
|
|
case 4:
|
|
while(std::getline(ss, item, ',')) {
|
|
comp.elements.push_back(item);
|
|
}
|
|
break;
|
|
case 5:
|
|
while(std::getline(ss, item, ',')) {
|
|
val = std::pow(10.0,std::stod(item));
|
|
comp.abundances.push_back(val);
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
// find where the start of the scheme block is
|
|
auto start_pos = std::ranges::search(line, std::format("BEGIN {}", scheme));
|
|
|
|
if (!start_pos.empty()) {
|
|
start_line = i;
|
|
}
|
|
i+=1;
|
|
}
|
|
|
|
// if (start_pos==0):
|
|
// raise error ("Scheme {} not found", scheme)
|
|
|
|
|
|
if (debug){
|
|
std::println("he_abundance: {}", comp.he_abundance);
|
|
std::println("requires_atomic_weight: {}", comp.requires_atomic_weight);
|
|
std::println("elements: {}",comp.elements);
|
|
std::println("abundances: {}", comp.abundances);
|
|
}
|
|
|
|
return comp;
|
|
|
|
}
|
|
|
|
IsotopicPercentage ChemicalFileParser::parse_isotopic_percentage(const std::vector<char>& data,const std::string& scheme) const {
|
|
|
|
// get file and iso_scheme
|
|
// Load the file
|
|
// find the iso_scheme
|
|
// get iso_comp data
|
|
// IsotopicPercentage object
|
|
|
|
bool debug = false;
|
|
if (debug){
|
|
std::println("Parsing Isotopic Percentage for: {}", scheme);
|
|
}
|
|
|
|
std::istringstream stream(std::string(data.begin(), data.end()));
|
|
|
|
// add error message if something goes wrong
|
|
ParsedChemicalData parsed;
|
|
|
|
std::string line;
|
|
int start_line = 0;
|
|
int i = 0;
|
|
|
|
IsotopicPercentage iso;
|
|
|
|
while (std::getline(stream, line)) {
|
|
|
|
// find where the end of the scheme block is
|
|
auto end_pos = std::ranges::search(line,std::format("END {}", scheme));
|
|
|
|
// exit if have reached the end of block
|
|
if (!end_pos.empty()) {
|
|
break;
|
|
}
|
|
|
|
if (start_line>0){
|
|
const size_t colon_pos = line.find(':');
|
|
line = line.substr(colon_pos+1);
|
|
|
|
line.erase(std::remove_if(line.begin(), line.end(),
|
|
[](char c){ return c == '[' || c == ']'; }),
|
|
line.end());
|
|
|
|
trim(line);
|
|
std::string item;
|
|
std::stringstream ss(line);
|
|
parsed.push_back(line);
|
|
switch(i-start_line){
|
|
case 1:
|
|
iso.comment_str = line;
|
|
break;
|
|
case 2:
|
|
while(std::getline(ss, item, ',')) {
|
|
iso.atomic_numbers.push_back(std::stoi(item));
|
|
}
|
|
break;
|
|
case 3:
|
|
while(std::getline(ss, item, ',')) {
|
|
iso.elements.push_back(item);
|
|
}
|
|
break;
|
|
case 4:
|
|
while(std::getline(ss, item, ',')) {
|
|
iso.mass_numbers.push_back(std::stoi(item));
|
|
}
|
|
break;
|
|
case 5:
|
|
while(std::getline(ss, item, ',')) {
|
|
iso.percentages.push_back(std::stod(item));
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
// find where the start of the scheme block is
|
|
auto start_pos = std::ranges::search(line, std::format("BEGIN {}", scheme));
|
|
|
|
if (!start_pos.empty()) {
|
|
start_line = i;
|
|
}
|
|
i+=1;
|
|
}
|
|
|
|
if (debug){
|
|
std::println("atomic_numbers: {}", iso.atomic_numbers);
|
|
std::println("elements: {}",iso.elements);
|
|
std::println("mass_numbers: {}", iso.mass_numbers);
|
|
std::println("percentages: {}", iso.percentages);
|
|
}
|
|
|
|
return iso;
|
|
}
|
|
}
|
|
|
|
namespace fourdst::composition {
|
|
Composition get_composition_record(const std::string& metal_fraction_scheme,
|
|
const std::string& isotopic_percentage_scheme,
|
|
double initial_z, double initial_y) {
|
|
|
|
std::vector<char> data;
|
|
|
|
io::ChemicalFileParser parser;
|
|
io::CompositionData compositions;
|
|
|
|
io::IsotopicPercentage isotopes;
|
|
|
|
data = std::ranges::to<std::vector<char>>(StandardAbundances);
|
|
|
|
compositions = parser.parse_compositon_data(data,metal_fraction_scheme);
|
|
isotopes = parser.parse_isotopic_percentage(data,isotopic_percentage_scheme);
|
|
|
|
std::string name;
|
|
std::vector<fourdst::atomic::Species> species;
|
|
|
|
// construct name of the isotopes
|
|
for (const auto [E,A] : std::ranges::views::zip(isotopes.elements, isotopes.mass_numbers)){
|
|
name = std::format("{}-{}",E,A);
|
|
// std::println("{}", name);
|
|
auto SpeciesObject = fourdst::atomic::species.at(name);
|
|
species.push_back(SpeciesObject);
|
|
|
|
// std::println("Species: {} has mass: {}", SpeciesObject.name(), SpeciesObject.mass());
|
|
|
|
}
|
|
|
|
std::vector<double> massFracs, metal_fractions;
|
|
|
|
// hydrogen and helium are treated separately since they are not metals
|
|
// H1
|
|
massFracs.push_back(std::max(0.0, std::min(1.0, 1.0 - (initial_z + initial_y))));
|
|
// H2
|
|
massFracs.push_back(0.0);
|
|
// He3
|
|
// anders & grevesse 1989 solar mass fractions
|
|
double xsol_he3,xsol_he4;
|
|
xsol_he3=2.9291e-05;
|
|
xsol_he4=2.7521e-01;
|
|
massFracs.push_back(initial_y*xsol_he3/(xsol_he3 + xsol_he4));
|
|
// He4
|
|
massFracs.push_back(initial_y*xsol_he4/(xsol_he3 + xsol_he4));
|
|
// Metals
|
|
// multiply by atomic weight if needed
|
|
|
|
if (compositions.requires_atomic_weight){
|
|
// get isotope with max abundance for each metal
|
|
std::vector<double> element_atomic_weight;
|
|
|
|
for (size_t i = 0; i < isotopes.atomic_numbers.size();) {
|
|
size_t Z = isotopes.atomic_numbers[i];
|
|
if (Z<2) {
|
|
continue;
|
|
}
|
|
double max_iso = isotopes.percentages[i];
|
|
size_t loc = i;
|
|
size_t j = i ;
|
|
while (j<isotopes.atomic_numbers.size() && isotopes.atomic_numbers[j] == Z) {
|
|
if (isotopes.percentages[j]>max_iso) {
|
|
loc = j;
|
|
}
|
|
++j;
|
|
}
|
|
element_atomic_weight.push_back(species[loc].mass());
|
|
i=j;
|
|
}
|
|
|
|
for(size_t i=0;i < compositions.abundances.size();++i){
|
|
metal_fractions.push_back(element_atomic_weight[i]*compositions.abundances[i]);
|
|
}
|
|
} else {
|
|
metal_fractions = compositions.abundances;
|
|
}
|
|
|
|
double sum = std::accumulate(metal_fractions.begin(),metal_fractions.end(),0.0);
|
|
|
|
for(size_t i=0;i < metal_fractions.size();++i){
|
|
std::println("testing: {} , {}", i,compositions.elements[i]);
|
|
metal_fractions[i] = metal_fractions[i]/sum;
|
|
}
|
|
|
|
// Z = max(0d0, min(1d0, 1d0 - (h1 + h2 + he3 + he4)))
|
|
|
|
for (size_t i = 0; i < isotopes.atomic_numbers.size();++i) {
|
|
if (isotopes.atomic_numbers[i]<2) {
|
|
continue;
|
|
}
|
|
double frac = 1e-2*isotopes.percentages[i]*species[i].mass();
|
|
|
|
}
|
|
|
|
//
|
|
Composition comp = buildCompositionFromMassFractions(species, massFracs);
|
|
return comp;
|
|
|
|
}
|
|
|
|
// Composition get_composition_record(const io::SolarCompositions metal_fraction_scheme,
|
|
// const io::IsotopicPercentages isotopic_percentage_scheme,
|
|
// double initial_z,
|
|
// double initial_y) {
|
|
// return get_composition_record(
|
|
// io::SolarComposition_to_string_map.at(metal_fraction_scheme),
|
|
// io::IsotopicPercentages_to_string.at(isotopic_percentage_scheme),
|
|
// initial_z,
|
|
// initial_y
|
|
// );
|
|
// }
|
|
}
|