Skip to content

Commit

Permalink
Update input.cc
Browse files Browse the repository at this point in the history
  • Loading branch information
lukeshingles committed Aug 30, 2024
1 parent 720d3c5 commit 1574f9d
Showing 1 changed file with 8 additions and 8 deletions.
16 changes: 8 additions & 8 deletions input.cc
Original file line number Diff line number Diff line change
Expand Up @@ -186,8 +186,8 @@ void read_phixs_data_table(std::fstream &phixsfile, const int nphixspoints_input

const double nu_edge = (epsilon(element, upperion, 0) - epsilon(element, lowerion, lowerlevel)) / H;

auto nutable = std::vector<double>(nphixspoints_inputtable);
auto phixstable = std::vector<double>(nphixspoints_inputtable);
auto nugrid_in = std::vector<double>(nphixspoints_inputtable);
auto phixs_in = std::vector<double>(nphixspoints_inputtable);

for (int i = 0; i < nphixspoints_inputtable; i++) {
double energy = -1.;
Expand All @@ -196,23 +196,23 @@ void read_phixs_data_table(std::fstream &phixsfile, const int nphixspoints_input
assert_always(std::stringstream(phixsline) >> energy >> phixs);
assert_always(energy >= 0);
assert_always(phixs >= 0);
nutable[i] = nu_edge + (energy * 13.6 * EV) / H;
nugrid_in[i] = nu_edge + (energy * 13.6 * EV) / H;
// the photoionisation cross-sections in the database are given in Mbarn=1e6 * 1e-28m^2
// to convert to cgs units multiply by 1e-18
phixstable[i] = phixs * 1e-18;
phixs_in[i] = phixs * 1e-18;
}
const double nu_max = nutable[nphixspoints_inputtable - 1];
const double nu_max = nugrid_in.back();

// Now interpolate these cross-sections
levelphixstable[0] = phixstable[0];
levelphixstable[0] = phixs_in[0];

gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_spline *spline = gsl_spline_alloc(gsl_interp_linear, nphixspoints_inputtable);
gsl_spline_init(spline, nutable.data(), phixstable.data(), nphixspoints_inputtable);
gsl_spline_init(spline, nugrid_in.data(), phixs_in.data(), nphixspoints_inputtable);
for (int i = 1; i < globals::NPHIXSPOINTS; i++) {
const double nu = nu_edge * (1. + i * globals::NPHIXSNUINCREMENT);
if (nu > nu_max) {
levelphixstable[i] = phixstable[nphixspoints_inputtable - 1] * pow(nu_max / nu, 3);
levelphixstable[i] = phixs_in[nphixspoints_inputtable - 1] * pow(nu_max / nu, 3);
} else {
levelphixstable[i] = gsl_spline_eval(spline, nu, acc);
}
Expand Down

0 comments on commit 1574f9d

Please sign in to comment.