From 44492ffe3342ac0d04b3df68e78315c35552705b Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Wed, 13 Dec 2023 12:56:05 -1000 Subject: [PATCH 01/24] Add functionality to calculate abundance factor --- source/helper.h | 2 ++ source/loop.cpp | 61 +++++++++++++++++++++++++++++++++++++++++++++---- source/loop.h | 14 ++++++++++++ 3 files changed, 73 insertions(+), 4 deletions(-) diff --git a/source/helper.h b/source/helper.h index 52fbddd..c04166a 100644 --- a/source/helper.h +++ b/source/helper.h @@ -52,6 +52,8 @@ struct Parameters { bool save_terms; /* Switch for using the adaptive solver option */ bool use_adaptive_solver; + /* Switch for using time-variable abundances */ + bool use_variable_abundances; /* Path to output file */ std::string output_filename; /* XML node holding DEM calculation parameters */ diff --git a/source/loop.cpp b/source/loop.cpp index a44eddb..60cb088 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -192,7 +192,16 @@ void Loop::CalculateDerivs(const state_type &state, state_type &derivs, double t double f_e = CalculateThermalConduction(state[3],state[2],"electron"); double f_i = CalculateThermalConduction(state[4],state[2],"ion"); - double radiative_loss = CalculateRadiativeLoss(state[3]); + double radiative_loss; + if use_variable_abundances + { + double abundance_factor = CalculateAbundanceFactor(state[2], results.density[0]); + radiative_loss = CalculateRadiativeLoss(state[3], abundance_factor); + } + else + { + radiative_loss = CalculateRadiativeLoss(state[3]); + } double heat = heater->Get_Heating(time); double c1 = CalculateC1(state[3],state[4],state[2]); double c2 = CalculateC2(); @@ -263,7 +272,16 @@ void Loop::SaveTerms(void) double f_e = CalculateThermalConduction(__state[3], __state[2], "electron"); double f_i = CalculateThermalConduction(__state[4], __state[2], "ion"); double c1 = CalculateC1(__state[3], __state[4], __state[2]); - double radiative_loss = CalculateRadiativeLoss(__state[3]); + double radiative_loss; + if use_variable_abundances + { + double abundance_factor = CalculateAbundanceFactor(__state[2], results.density[0]); + radiative_loss = CalculateRadiativeLoss(__state[3], abundance_factor); + } + else + { + radiative_loss = CalculateRadiativeLoss(__state[3]); + } // Save terms terms.f_e.push_back(f_e); @@ -350,6 +368,21 @@ double Loop::CalculateRadiativeLoss(double temperature) return chi * std::pow( 10.0, (alpha*log_temperature) ); } +double Loop::CalculateRadiativeLoss(double temperature, double abundance_factor) +{ + +} + +double Loop::CalculateAbundanceFactor(double density, double initial_density); +{ + double initial_abundance_factor = 4.0; // Assumes initially coronal plasma + + // Calculate using a weighted average of the density + // AF = 1.0 + (AF_0 - 1) * (n_0 / n) + return 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); + +} + double Loop::CalculateCollisionFrequency(double temperature_e,double density) { // TODO: find a reference for this formula @@ -367,7 +400,17 @@ double Loop::CalculateC1(double temperature_e, double temperature_i, double dens double grav_correction = 1.0; double loss_correction = 1.0; double scale_height = CalculateScaleHeight(temperature_e,temperature_i); - double radiative_loss = CalculateRadiativeLoss(temperature_e); + double radiative_loss; + if use_variable_abundances + { + double abundance_factor = CalculateAbundanceFactor(density, results.density[0]); + radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); + } + else + { + radiative_loss = CalculateRadiativeLoss(temperature_e); + } + if(parameters.use_c1_grav_correction) { @@ -425,7 +468,17 @@ double Loop::CalculateVelocity(double temperature_e, double temperature_i, doubl double c4 = CalculateC4(); double density = pressure_e/(BOLTZMANN_CONSTANT*temperature_e); double c1 = CalculateC1(temperature_e,temperature_i,density); - double R_tr = c1*std::pow(density,2)*CalculateRadiativeLoss(temperature_e)*parameters.loop_length; + double radiative_loss; + if use_variable_abundances + { + double abundance_factor = CalculateAbundanceFactor(density, results.density[0]); + radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); + } + else + { + radiative_loss = CalculateRadiativeLoss(temperature_e); + } + double R_tr = c1*std::pow(density,2)*radiative_loss*parameters.loop_length; double fe = CalculateThermalConduction(temperature_e,density,"electron"); double fi = CalculateThermalConduction(temperature_i,density,"ion"); double sc = CalculateScaleHeight(temperature_e,temperature_i); diff --git a/source/loop.h b/source/loop.h index c23d9f4..17944ef 100644 --- a/source/loop.h +++ b/source/loop.h @@ -203,9 +203,13 @@ class Loop { // communication) and twice the coronal abundances of Meyer (1985). This is the same power-law // radiative loss function as is implemented in the HYDRAD code and the EBTEL IDL code. // + // The overloaded function uses the abundance factor to adjust the radiative loss curve for abundances + // that are not strictly coronal, as in the original function. + // // @return radiative loss function (in erg cm$^3$ s$^{-1}$) // static double CalculateRadiativeLoss(double temperature); + static double CalculateRadiativeLoss(double temperature, double abundance_factor); // Calculate derivatives of EBTEL equations // @state current state of the loop @@ -229,6 +233,16 @@ class Loop { // @return time until next change in the loop heating (in s) // static double CalculateTimeNextHeating(double time); + + // Calculate the current abundance factor + // + // Calculates the abundance factor as it varies due to filling and draining of the loop + // from chromospheric evaporation, which is assumed to bring photospheric material (AF = 1.0). It is + // also assumed to be initially coronal (AF = 4.0) + // + // @return the abundance factor (unitless) + static double CalculateAbundanceFactor(double density, double initial_density); + }; // Pointer to the class typedef Loop* LOOP; From 1fb037d92f3a943cc1fd29bec1c71f4db41ae0ff Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 18 Dec 2023 13:13:03 -1000 Subject: [PATCH 02/24] add search function for nearest value in array --- source/loop.cpp | 9 +++++++-- source/loop.h | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/source/loop.cpp b/source/loop.cpp index 60cb088..8c8ab67 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -376,11 +376,16 @@ double Loop::CalculateRadiativeLoss(double temperature, double abundance_factor) double Loop::CalculateAbundanceFactor(double density, double initial_density); { double initial_abundance_factor = 4.0; // Assumes initially coronal plasma - + double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for radiative loss function + int array_length = sizeof(abundance_array) / sizeof(abundance_array[0]); + // Calculate using a weighted average of the density // AF = 1.0 + (AF_0 - 1) * (n_0 / n) - return 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); + double abundance_factor = 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); + // Find the nearest value of AF to the values in our discrete list, then return that value + int index = find_closest(abundance_factor, abundance_array, array_length); + return abundance_array[index]; } double Loop::CalculateCollisionFrequency(double temperature_e,double density) diff --git a/source/loop.h b/source/loop.h index 17944ef..aaced7f 100644 --- a/source/loop.h +++ b/source/loop.h @@ -9,6 +9,7 @@ Loop class definition #include "helper.h" #include "heater.h" #include "util/constants.h" +#include "util/misc.h" // Loop object // From 001871c17fd455840176a381f243944e0c5f43e1 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 18 Dec 2023 13:14:16 -1000 Subject: [PATCH 03/24] add misc.cpp and misc.h --- source/util/misc.cpp | 27 +++++++++++++++++++++++++++ source/util/misc.h | 12 ++++++++++++ 2 files changed, 39 insertions(+) create mode 100644 source/util/misc.cpp create mode 100644 source/util/misc.h diff --git a/source/util/misc.cpp b/source/util/misc.cpp new file mode 100644 index 0000000..97b78f6 --- /dev/null +++ b/source/util/misc.cpp @@ -0,0 +1,27 @@ +/* +misc.cpp +Miscellaneous functions for manipulating data +*/ + +#include "misc.h" + + +int find_closest(double x, double array[], int array_length) +{ + /* Traverses the whole array, so O(N) search. + * A binary search would be more efficient for long arrays O(ln N). + */ + int index = 0; + double closest_value = array[0]; + + for( int i=1; i < array_length; ++i ) + { + if( abs(x - closest_value) > abs(x - array[i]) ) + { + closest_value = array[i]; + index = i; + } + } + + return index; +} diff --git a/source/util/misc.h b/source/util/misc.h new file mode 100644 index 0000000..261b622 --- /dev/null +++ b/source/util/misc.h @@ -0,0 +1,12 @@ +// **** +// * +// * Header file for miscellaneous functions +// * +// **** + + +// Find the closest value in an array +int find_closest(double x, double array[], int array_length); + +// Read csv data + From 4a044aac8fcc2a80bf044dd8aebf15222f69793c Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 18 Dec 2023 13:17:03 -1000 Subject: [PATCH 04/24] add radiative loss data --- data/radiation/abund_10_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_15_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_20_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_25_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_30_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_35_rad_loss.dat | 102 +++++++++++++++++++++++++++ data/radiation/abund_40_rad_loss.dat | 102 +++++++++++++++++++++++++++ 7 files changed, 714 insertions(+) create mode 100644 data/radiation/abund_10_rad_loss.dat create mode 100644 data/radiation/abund_15_rad_loss.dat create mode 100644 data/radiation/abund_20_rad_loss.dat create mode 100644 data/radiation/abund_25_rad_loss.dat create mode 100644 data/radiation/abund_30_rad_loss.dat create mode 100644 data/radiation/abund_35_rad_loss.dat create mode 100644 data/radiation/abund_40_rad_loss.dat diff --git a/data/radiation/abund_10_rad_loss.dat b/data/radiation/abund_10_rad_loss.dat new file mode 100644 index 0000000..bc07666 --- /dev/null +++ b/data/radiation/abund_10_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.309,-23.274,-23.265,-23.259,-23.259,-23.260,-23.260 +04.05,-22.814,-22.770,-22.760,-22.751,-22.751,-22.751,-22.751 +04.10,-22.365,-22.318,-22.307,-22.299,-22.298,-22.298,-22.298 +04.15,-22.007,-21.960,-21.951,-21.943,-21.942,-21.942,-21.942 +04.20,-21.811,-21.768,-21.760,-21.753,-21.753,-21.753,-21.753 +04.25,-21.791,-21.753,-21.746,-21.742,-21.741,-21.742,-21.742 +04.30,-21.866,-21.834,-21.829,-21.827,-21.827,-21.828,-21.829 +04.35,-21.962,-21.937,-21.933,-21.933,-21.935,-21.937,-21.938 +04.40,-22.035,-22.017,-22.015,-22.018,-22.021,-22.025,-22.026 +04.45,-22.067,-22.055,-22.055,-22.061,-22.064,-22.070,-22.072 +04.50,-22.068,-22.061,-22.063,-22.071,-22.075,-22.082,-22.084 +04.55,-22.048,-22.046,-22.048,-22.059,-22.063,-22.071,-22.074 +04.60,-22.007,-22.009,-22.012,-22.025,-22.030,-22.039,-22.043 +04.65,-21.943,-21.948,-21.952,-21.967,-21.973,-21.983,-21.987 +04.70,-21.858,-21.864,-21.869,-21.884,-21.891,-21.901,-21.905 +04.75,-21.757,-21.764,-21.769,-21.784,-21.789,-21.798,-21.801 +04.80,-21.650,-21.656,-21.661,-21.674,-21.678,-21.683,-21.685 +04.85,-21.550,-21.556,-21.560,-21.571,-21.574,-21.576,-21.577 +04.90,-21.477,-21.482,-21.486,-21.495,-21.497,-21.498,-21.499 +04.95,-21.430,-21.435,-21.438,-21.446,-21.448,-21.450,-21.452 +05.00,-21.405,-21.409,-21.412,-21.419,-21.421,-21.425,-21.427 +05.05,-21.413,-21.416,-21.418,-21.423,-21.425,-21.431,-21.435 +05.10,-21.431,-21.433,-21.434,-21.438,-21.440,-21.449,-21.454 +05.15,-21.426,-21.428,-21.429,-21.432,-21.434,-21.446,-21.453 +05.20,-21.405,-21.407,-21.408,-21.410,-21.413,-21.426,-21.434 +05.25,-21.386,-21.388,-21.389,-21.391,-21.394,-21.407,-21.417 +05.30,-21.379,-21.381,-21.381,-21.384,-21.386,-21.399,-21.409 +05.35,-21.376,-21.378,-21.379,-21.381,-21.383,-21.396,-21.406 +05.40,-21.386,-21.387,-21.388,-21.390,-21.391,-21.403,-21.412 +05.45,-21.442,-21.444,-21.444,-21.446,-21.447,-21.456,-21.464 +05.50,-21.565,-21.566,-21.567,-21.568,-21.569,-21.575,-21.580 +05.55,-21.696,-21.697,-21.698,-21.699,-21.700,-21.703,-21.706 +05.60,-21.776,-21.778,-21.778,-21.780,-21.780,-21.782,-21.784 +05.65,-21.807,-21.808,-21.809,-21.810,-21.811,-21.812,-21.814 +05.70,-21.815,-21.816,-21.817,-21.818,-21.819,-21.821,-21.821 +05.75,-21.826,-21.827,-21.828,-21.830,-21.830,-21.832,-21.833 +05.80,-21.849,-21.850,-21.851,-21.853,-21.854,-21.856,-21.857 +05.85,-21.868,-21.869,-21.870,-21.873,-21.875,-21.877,-21.877 +05.90,-21.872,-21.874,-21.875,-21.879,-21.880,-21.883,-21.884 +05.95,-21.869,-21.871,-21.873,-21.877,-21.879,-21.882,-21.883 +06.00,-21.872,-21.874,-21.875,-21.879,-21.881,-21.884,-21.885 +06.05,-21.884,-21.886,-21.887,-21.891,-21.893,-21.896,-21.897 +06.10,-21.902,-21.903,-21.904,-21.908,-21.909,-21.912,-21.913 +06.15,-21.922,-21.923,-21.924,-21.927,-21.928,-21.930,-21.931 +06.20,-21.950,-21.951,-21.952,-21.954,-21.955,-21.956,-21.957 +06.25,-21.994,-21.995,-21.996,-21.998,-21.998,-22.000,-22.000 +06.30,-22.065,-22.066,-22.066,-22.068,-22.069,-22.070,-22.070 +06.35,-22.162,-22.163,-22.163,-22.165,-22.165,-22.167,-22.167 +06.40,-22.265,-22.266,-22.266,-22.268,-22.268,-22.269,-22.269 +06.45,-22.359,-22.359,-22.359,-22.360,-22.361,-22.362,-22.362 +06.50,-22.434,-22.434,-22.434,-22.435,-22.435,-22.436,-22.436 +06.55,-22.487,-22.487,-22.487,-22.487,-22.487,-22.488,-22.488 +06.60,-22.519,-22.519,-22.519,-22.519,-22.519,-22.520,-22.520 +06.65,-22.535,-22.535,-22.535,-22.535,-22.535,-22.536,-22.536 +06.70,-22.539,-22.539,-22.539,-22.539,-22.540,-22.540,-22.540 +06.75,-22.536,-22.536,-22.536,-22.536,-22.536,-22.536,-22.536 +06.80,-22.529,-22.529,-22.529,-22.529,-22.529,-22.529,-22.529 +06.85,-22.521,-22.521,-22.521,-22.521,-22.521,-22.521,-22.521 +06.90,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516 +06.95,-22.518,-22.518,-22.518,-22.518,-22.518,-22.518,-22.518 +07.00,-22.529,-22.530,-22.530,-22.530,-22.530,-22.530,-22.530 +07.05,-22.554,-22.554,-22.554,-22.554,-22.554,-22.554,-22.554 +07.10,-22.592,-22.592,-22.592,-22.592,-22.592,-22.592,-22.593 +07.15,-22.641,-22.641,-22.641,-22.641,-22.641,-22.641,-22.641 +07.20,-22.689,-22.689,-22.689,-22.689,-22.689,-22.689,-22.689 +07.25,-22.727,-22.727,-22.727,-22.727,-22.727,-22.727,-22.727 +07.30,-22.752,-22.752,-22.752,-22.752,-22.752,-22.752,-22.752 +07.35,-22.766,-22.766,-22.766,-22.766,-22.766,-22.766,-22.766 +07.40,-22.771,-22.771,-22.771,-22.771,-22.771,-22.770,-22.770 +07.45,-22.768,-22.768,-22.768,-22.768,-22.768,-22.768,-22.768 +07.50,-22.762,-22.762,-22.762,-22.762,-22.762,-22.762,-22.762 +07.55,-22.751,-22.751,-22.751,-22.751,-22.751,-22.751,-22.751 +07.60,-22.738,-22.738,-22.738,-22.738,-22.738,-22.738,-22.738 +07.65,-22.723,-22.723,-22.723,-22.723,-22.723,-22.723,-22.723 +07.70,-22.706,-22.706,-22.706,-22.706,-22.706,-22.706,-22.706 +07.75,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688 +07.80,-22.670,-22.670,-22.670,-22.670,-22.670,-22.670,-22.670 +07.85,-22.651,-22.651,-22.651,-22.651,-22.651,-22.651,-22.651 +07.90,-22.632,-22.632,-22.632,-22.632,-22.632,-22.632,-22.632 +07.95,-22.612,-22.612,-22.612,-22.612,-22.612,-22.612,-22.612 +08.00,-22.592,-22.592,-22.592,-22.592,-22.592,-22.592,-22.592 +08.05,-22.572,-22.572,-22.572,-22.572,-22.572,-22.572,-22.572 +08.10,-22.552,-22.552,-22.552,-22.552,-22.552,-22.552,-22.552 +08.15,-22.531,-22.531,-22.531,-22.531,-22.531,-22.531,-22.531 +08.20,-22.510,-22.510,-22.510,-22.510,-22.510,-22.510,-22.510 +08.25,-22.488,-22.488,-22.488,-22.488,-22.488,-22.488,-22.488 +08.30,-22.467,-22.467,-22.467,-22.467,-22.467,-22.467,-22.467 +08.35,-22.444,-22.444,-22.444,-22.444,-22.444,-22.444,-22.444 +08.40,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422 +08.45,-22.399,-22.399,-22.399,-22.399,-22.399,-22.399,-22.399 +08.50,-22.376,-22.376,-22.376,-22.376,-22.376,-22.376,-22.376 +08.55,-22.353,-22.353,-22.353,-22.353,-22.353,-22.353,-22.353 +08.60,-22.330,-22.330,-22.330,-22.330,-22.330,-22.330,-22.330 +08.65,-22.306,-22.306,-22.306,-22.306,-22.306,-22.306,-22.306 +08.70,-22.283,-22.283,-22.283,-22.283,-22.283,-22.283,-22.283 +08.75,-22.259,-22.259,-22.259,-22.259,-22.259,-22.259,-22.259 +08.80,-22.235,-22.235,-22.235,-22.235,-22.235,-22.235,-22.235 +08.85,-22.211,-22.211,-22.211,-22.211,-22.211,-22.211,-22.211 +08.90,-22.187,-22.187,-22.187,-22.187,-22.187,-22.187,-22.187 +08.95,-22.163,-22.163,-22.163,-22.163,-22.163,-22.163,-22.163 +09.00,-22.138,-22.138,-22.138,-22.138,-22.138,-22.138,-22.138 diff --git a/data/radiation/abund_15_rad_loss.dat b/data/radiation/abund_15_rad_loss.dat new file mode 100644 index 0000000..a3e9125 --- /dev/null +++ b/data/radiation/abund_15_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.249,-23.220,-23.213,-23.208,-23.208,-23.209,-23.210 +04.05,-22.775,-22.735,-22.726,-22.718,-22.718,-22.718,-22.719 +04.10,-22.340,-22.296,-22.286,-22.278,-22.277,-22.277,-22.278 +04.15,-21.990,-21.945,-21.936,-21.929,-21.928,-21.928,-21.928 +04.20,-21.796,-21.755,-21.747,-21.741,-21.740,-21.740,-21.740 +04.25,-21.775,-21.739,-21.732,-21.728,-21.728,-21.728,-21.728 +04.30,-21.848,-21.818,-21.813,-21.811,-21.811,-21.813,-21.813 +04.35,-21.941,-21.917,-21.914,-21.915,-21.916,-21.919,-21.920 +04.40,-22.010,-21.993,-21.992,-21.995,-21.998,-22.003,-22.005 +04.45,-22.036,-22.025,-22.025,-22.031,-22.035,-22.041,-22.044 +04.50,-22.030,-22.025,-22.026,-22.034,-22.037,-22.045,-22.048 +04.55,-22.005,-22.004,-22.006,-22.015,-22.019,-22.028,-22.032 +04.60,-21.962,-21.965,-21.968,-21.979,-21.983,-21.993,-21.997 +04.65,-21.900,-21.905,-21.909,-21.922,-21.927,-21.938,-21.942 +04.70,-21.817,-21.824,-21.828,-21.842,-21.848,-21.858,-21.863 +04.75,-21.720,-21.727,-21.731,-21.745,-21.751,-21.759,-21.763 +04.80,-21.618,-21.624,-21.628,-21.640,-21.645,-21.650,-21.652 +04.85,-21.525,-21.531,-21.535,-21.545,-21.548,-21.550,-21.551 +04.90,-21.462,-21.467,-21.470,-21.479,-21.481,-21.482,-21.483 +04.95,-21.422,-21.427,-21.430,-21.439,-21.441,-21.443,-21.444 +05.00,-21.401,-21.405,-21.408,-21.415,-21.417,-21.422,-21.424 +05.05,-21.410,-21.414,-21.415,-21.420,-21.422,-21.429,-21.433 +05.10,-21.428,-21.431,-21.432,-21.435,-21.438,-21.446,-21.452 +05.15,-21.423,-21.425,-21.426,-21.429,-21.432,-21.443,-21.450 +05.20,-21.401,-21.403,-21.404,-21.407,-21.410,-21.423,-21.431 +05.25,-21.382,-21.384,-21.385,-21.387,-21.390,-21.403,-21.412 +05.30,-21.373,-21.375,-21.376,-21.378,-21.381,-21.394,-21.403 +05.35,-21.369,-21.371,-21.372,-21.374,-21.376,-21.388,-21.398 +05.40,-21.376,-21.377,-21.378,-21.380,-21.382,-21.393,-21.402 +05.45,-21.427,-21.429,-21.430,-21.431,-21.433,-21.442,-21.449 +05.50,-21.539,-21.541,-21.541,-21.543,-21.544,-21.550,-21.555 +05.55,-21.653,-21.654,-21.655,-21.656,-21.657,-21.661,-21.663 +05.60,-21.713,-21.714,-21.715,-21.717,-21.718,-21.720,-21.721 +05.65,-21.726,-21.728,-21.728,-21.730,-21.731,-21.733,-21.734 +05.70,-21.719,-21.721,-21.721,-21.723,-21.724,-21.726,-21.727 +05.75,-21.715,-21.716,-21.717,-21.719,-21.720,-21.722,-21.723 +05.80,-21.720,-21.722,-21.723,-21.725,-21.727,-21.729,-21.729 +05.85,-21.723,-21.725,-21.726,-21.730,-21.731,-21.734,-21.734 +05.90,-21.718,-21.720,-21.721,-21.725,-21.727,-21.730,-21.731 +05.95,-21.710,-21.712,-21.714,-21.718,-21.720,-21.724,-21.725 +06.00,-21.711,-21.713,-21.715,-21.719,-21.721,-21.725,-21.726 +06.05,-21.723,-21.726,-21.727,-21.731,-21.733,-21.736,-21.737 +06.10,-21.742,-21.744,-21.745,-21.749,-21.750,-21.753,-21.754 +06.15,-21.765,-21.766,-21.767,-21.770,-21.771,-21.773,-21.774 +06.20,-21.795,-21.796,-21.797,-21.800,-21.801,-21.802,-21.803 +06.25,-21.844,-21.845,-21.846,-21.848,-21.849,-21.850,-21.851 +06.30,-21.922,-21.923,-21.924,-21.926,-21.927,-21.928,-21.928 +06.35,-22.030,-22.031,-22.031,-22.033,-22.034,-22.035,-22.035 +06.40,-22.145,-22.145,-22.146,-22.147,-22.148,-22.149,-22.149 +06.45,-22.246,-22.247,-22.247,-22.248,-22.248,-22.249,-22.250 +06.50,-22.325,-22.326,-22.326,-22.326,-22.327,-22.328,-22.328 +06.55,-22.379,-22.379,-22.379,-22.380,-22.380,-22.380,-22.380 +06.60,-22.409,-22.409,-22.410,-22.410,-22.410,-22.410,-22.410 +06.65,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422 +06.70,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422 +06.75,-22.415,-22.415,-22.415,-22.415,-22.415,-22.415,-22.415 +06.80,-22.404,-22.404,-22.404,-22.404,-22.404,-22.404,-22.404 +06.85,-22.394,-22.394,-22.394,-22.394,-22.394,-22.394,-22.394 +06.90,-22.388,-22.388,-22.388,-22.388,-22.388,-22.388,-22.388 +06.95,-22.390,-22.390,-22.390,-22.390,-22.390,-22.390,-22.390 +07.00,-22.404,-22.404,-22.404,-22.404,-22.404,-22.404,-22.404 +07.05,-22.432,-22.432,-22.433,-22.433,-22.433,-22.433,-22.433 +07.10,-22.479,-22.479,-22.479,-22.479,-22.479,-22.479,-22.479 +07.15,-22.539,-22.539,-22.539,-22.539,-22.539,-22.539,-22.539 +07.20,-22.600,-22.600,-22.600,-22.600,-22.600,-22.600,-22.600 +07.25,-22.651,-22.651,-22.651,-22.651,-22.651,-22.651,-22.651 +07.30,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688 +07.35,-22.711,-22.711,-22.711,-22.711,-22.711,-22.711,-22.711 +07.40,-22.723,-22.723,-22.723,-22.723,-22.723,-22.723,-22.723 +07.45,-22.727,-22.727,-22.727,-22.727,-22.727,-22.727,-22.727 +07.50,-22.724,-22.724,-22.724,-22.724,-22.724,-22.724,-22.724 +07.55,-22.717,-22.717,-22.717,-22.717,-22.717,-22.717,-22.717 +07.60,-22.707,-22.707,-22.707,-22.707,-22.707,-22.707,-22.707 +07.65,-22.693,-22.693,-22.693,-22.693,-22.693,-22.693,-22.693 +07.70,-22.679,-22.679,-22.679,-22.679,-22.679,-22.679,-22.679 +07.75,-22.662,-22.662,-22.662,-22.662,-22.662,-22.662,-22.662 +07.80,-22.645,-22.645,-22.645,-22.645,-22.645,-22.645,-22.645 +07.85,-22.628,-22.628,-22.628,-22.628,-22.628,-22.628,-22.628 +07.90,-22.610,-22.610,-22.610,-22.610,-22.610,-22.610,-22.610 +07.95,-22.592,-22.592,-22.592,-22.592,-22.592,-22.592,-22.592 +08.00,-22.573,-22.573,-22.573,-22.573,-22.573,-22.573,-22.573 +08.05,-22.554,-22.554,-22.554,-22.554,-22.554,-22.554,-22.554 +08.10,-22.535,-22.535,-22.535,-22.535,-22.535,-22.535,-22.535 +08.15,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516 +08.20,-22.496,-22.496,-22.496,-22.496,-22.496,-22.496,-22.496 +08.25,-22.475,-22.475,-22.475,-22.475,-22.475,-22.475,-22.475 +08.30,-22.454,-22.454,-22.454,-22.454,-22.454,-22.454,-22.454 +08.35,-22.433,-22.433,-22.433,-22.433,-22.433,-22.433,-22.433 +08.40,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411 +08.45,-22.389,-22.389,-22.389,-22.389,-22.389,-22.389,-22.389 +08.50,-22.367,-22.367,-22.367,-22.367,-22.367,-22.367,-22.367 +08.55,-22.344,-22.344,-22.344,-22.344,-22.344,-22.344,-22.344 +08.60,-22.321,-22.321,-22.321,-22.321,-22.321,-22.321,-22.321 +08.65,-22.298,-22.298,-22.298,-22.298,-22.298,-22.298,-22.298 +08.70,-22.275,-22.275,-22.275,-22.275,-22.275,-22.275,-22.275 +08.75,-22.251,-22.251,-22.251,-22.251,-22.251,-22.251,-22.251 +08.80,-22.227,-22.227,-22.227,-22.227,-22.227,-22.227,-22.227 +08.85,-22.204,-22.204,-22.204,-22.204,-22.204,-22.204,-22.204 +08.90,-22.180,-22.180,-22.180,-22.180,-22.180,-22.180,-22.180 +08.95,-22.156,-22.156,-22.156,-22.156,-22.156,-22.156,-22.156 +09.00,-22.131,-22.131,-22.131,-22.131,-22.131,-22.131,-22.131 diff --git a/data/radiation/abund_20_rad_loss.dat b/data/radiation/abund_20_rad_loss.dat new file mode 100644 index 0000000..c0ab712 --- /dev/null +++ b/data/radiation/abund_20_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.199,-23.175,-23.169,-23.164,-23.165,-23.166,-23.167 +04.05,-22.741,-22.705,-22.696,-22.690,-22.689,-22.690,-22.690 +04.10,-22.318,-22.276,-22.267,-22.259,-22.259,-22.259,-22.259 +04.15,-21.974,-21.931,-21.923,-21.915,-21.915,-21.915,-21.915 +04.20,-21.783,-21.743,-21.735,-21.729,-21.729,-21.729,-21.729 +04.25,-21.761,-21.726,-21.719,-21.715,-21.715,-21.716,-21.716 +04.30,-21.832,-21.803,-21.798,-21.796,-21.797,-21.799,-21.800 +04.35,-21.923,-21.900,-21.897,-21.898,-21.900,-21.903,-21.904 +04.40,-21.988,-21.972,-21.970,-21.974,-21.977,-21.983,-21.985 +04.45,-22.009,-21.999,-21.999,-22.005,-22.009,-22.016,-22.019 +04.50,-21.998,-21.993,-21.994,-22.001,-22.005,-22.013,-22.017 +04.55,-21.969,-21.968,-21.969,-21.978,-21.982,-21.991,-21.995 +04.60,-21.925,-21.927,-21.930,-21.939,-21.944,-21.954,-21.958 +04.65,-21.863,-21.868,-21.872,-21.884,-21.889,-21.899,-21.904 +04.70,-21.783,-21.789,-21.793,-21.806,-21.812,-21.822,-21.827 +04.75,-21.688,-21.695,-21.699,-21.712,-21.717,-21.726,-21.730 +04.80,-21.589,-21.595,-21.599,-21.611,-21.615,-21.621,-21.624 +04.85,-21.503,-21.508,-21.512,-21.522,-21.525,-21.527,-21.529 +04.90,-21.447,-21.453,-21.456,-21.465,-21.467,-21.468,-21.469 +04.95,-21.415,-21.420,-21.423,-21.432,-21.434,-21.436,-21.437 +05.00,-21.397,-21.402,-21.405,-21.412,-21.414,-21.418,-21.421 +05.05,-21.408,-21.411,-21.413,-21.418,-21.420,-21.427,-21.430 +05.10,-21.426,-21.429,-21.430,-21.433,-21.435,-21.444,-21.450 +05.15,-21.420,-21.423,-21.424,-21.427,-21.429,-21.440,-21.447 +05.20,-21.398,-21.400,-21.401,-21.404,-21.407,-21.419,-21.428 +05.25,-21.378,-21.380,-21.381,-21.383,-21.386,-21.399,-21.408 +05.30,-21.368,-21.370,-21.371,-21.373,-21.376,-21.389,-21.398 +05.35,-21.362,-21.364,-21.365,-21.367,-21.369,-21.382,-21.391 +05.40,-21.367,-21.368,-21.369,-21.371,-21.373,-21.384,-21.393 +05.45,-21.414,-21.416,-21.416,-21.418,-21.420,-21.428,-21.435 +05.50,-21.517,-21.518,-21.519,-21.521,-21.522,-21.528,-21.532 +05.55,-21.616,-21.617,-21.618,-21.620,-21.621,-21.624,-21.627 +05.60,-21.661,-21.662,-21.663,-21.665,-21.666,-21.668,-21.670 +05.65,-21.662,-21.663,-21.664,-21.666,-21.667,-21.669,-21.670 +05.70,-21.645,-21.647,-21.648,-21.650,-21.651,-21.653,-21.653 +05.75,-21.632,-21.633,-21.634,-21.636,-21.638,-21.640,-21.640 +05.80,-21.627,-21.628,-21.629,-21.632,-21.634,-21.636,-21.637 +05.85,-21.621,-21.623,-21.624,-21.628,-21.630,-21.632,-21.633 +05.90,-21.610,-21.613,-21.614,-21.618,-21.620,-21.623,-21.624 +05.95,-21.600,-21.603,-21.604,-21.609,-21.611,-21.614,-21.615 +06.00,-21.600,-21.603,-21.604,-21.609,-21.611,-21.614,-21.616 +06.05,-21.613,-21.615,-21.617,-21.621,-21.623,-21.626,-21.627 +06.10,-21.633,-21.634,-21.635,-21.639,-21.641,-21.644,-21.645 +06.15,-21.656,-21.657,-21.658,-21.661,-21.662,-21.665,-21.665 +06.20,-21.688,-21.689,-21.690,-21.692,-21.693,-21.695,-21.696 +06.25,-21.739,-21.740,-21.741,-21.743,-21.744,-21.745,-21.746 +06.30,-21.821,-21.822,-21.823,-21.825,-21.826,-21.827,-21.827 +06.35,-21.934,-21.935,-21.936,-21.938,-21.939,-21.940,-21.940 +06.40,-22.056,-22.056,-22.057,-22.058,-22.059,-22.060,-22.060 +06.45,-22.162,-22.162,-22.163,-22.164,-22.164,-22.165,-22.166 +06.50,-22.244,-22.244,-22.244,-22.245,-22.245,-22.246,-22.246 +06.55,-22.298,-22.298,-22.298,-22.298,-22.299,-22.299,-22.299 +06.60,-22.327,-22.327,-22.327,-22.327,-22.327,-22.328,-22.328 +06.65,-22.337,-22.337,-22.337,-22.337,-22.338,-22.338,-22.338 +06.70,-22.335,-22.335,-22.335,-22.335,-22.335,-22.335,-22.335 +06.75,-22.325,-22.325,-22.325,-22.325,-22.325,-22.325,-22.325 +06.80,-22.313,-22.313,-22.313,-22.313,-22.313,-22.313,-22.313 +06.85,-22.301,-22.301,-22.301,-22.301,-22.301,-22.301,-22.301 +06.90,-22.295,-22.295,-22.295,-22.295,-22.295,-22.295,-22.295 +06.95,-22.297,-22.297,-22.297,-22.297,-22.297,-22.297,-22.297 +07.00,-22.312,-22.312,-22.312,-22.312,-22.312,-22.312,-22.312 +07.05,-22.343,-22.343,-22.343,-22.343,-22.343,-22.343,-22.343 +07.10,-22.394,-22.394,-22.394,-22.394,-22.394,-22.394,-22.394 +07.15,-22.461,-22.461,-22.461,-22.461,-22.461,-22.461,-22.461 +07.20,-22.530,-22.530,-22.530,-22.530,-22.530,-22.530,-22.530 +07.25,-22.590,-22.590,-22.590,-22.590,-22.590,-22.590,-22.590 +07.30,-22.635,-22.635,-22.635,-22.635,-22.635,-22.635,-22.635 +07.35,-22.666,-22.666,-22.666,-22.665,-22.665,-22.665,-22.665 +07.40,-22.683,-22.683,-22.683,-22.683,-22.683,-22.683,-22.683 +07.45,-22.691,-22.691,-22.691,-22.691,-22.691,-22.691,-22.691 +07.50,-22.692,-22.692,-22.692,-22.692,-22.692,-22.692,-22.692 +07.55,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688,-22.688 +07.60,-22.679,-22.679,-22.679,-22.679,-22.679,-22.679,-22.679 +07.65,-22.668,-22.668,-22.668,-22.668,-22.668,-22.668,-22.668 +07.70,-22.654,-22.654,-22.654,-22.654,-22.654,-22.654,-22.654 +07.75,-22.639,-22.639,-22.639,-22.639,-22.639,-22.639,-22.639 +07.80,-22.624,-22.624,-22.624,-22.624,-22.624,-22.624,-22.624 +07.85,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607 +07.90,-22.590,-22.590,-22.590,-22.590,-22.590,-22.590,-22.590 +07.95,-22.573,-22.573,-22.573,-22.573,-22.573,-22.573,-22.573 +08.00,-22.556,-22.556,-22.556,-22.556,-22.556,-22.556,-22.556 +08.05,-22.538,-22.538,-22.538,-22.538,-22.538,-22.538,-22.538 +08.10,-22.520,-22.520,-22.520,-22.520,-22.520,-22.520,-22.520 +08.15,-22.502,-22.502,-22.502,-22.502,-22.502,-22.502,-22.502 +08.20,-22.483,-22.483,-22.483,-22.483,-22.483,-22.483,-22.483 +08.25,-22.463,-22.463,-22.463,-22.463,-22.463,-22.463,-22.463 +08.30,-22.443,-22.443,-22.443,-22.443,-22.443,-22.443,-22.443 +08.35,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422,-22.422 +08.40,-22.401,-22.401,-22.401,-22.401,-22.401,-22.401,-22.401 +08.45,-22.380,-22.380,-22.380,-22.380,-22.380,-22.380,-22.380 +08.50,-22.358,-22.358,-22.358,-22.358,-22.358,-22.358,-22.358 +08.55,-22.336,-22.336,-22.336,-22.336,-22.335,-22.335,-22.335 +08.60,-22.313,-22.313,-22.313,-22.313,-22.313,-22.313,-22.313 +08.65,-22.290,-22.290,-22.290,-22.290,-22.290,-22.290,-22.290 +08.70,-22.267,-22.267,-22.267,-22.267,-22.267,-22.267,-22.267 +08.75,-22.244,-22.244,-22.244,-22.244,-22.244,-22.244,-22.244 +08.80,-22.220,-22.220,-22.220,-22.220,-22.220,-22.220,-22.220 +08.85,-22.197,-22.197,-22.197,-22.197,-22.197,-22.197,-22.197 +08.90,-22.173,-22.173,-22.173,-22.173,-22.173,-22.173,-22.173 +08.95,-22.149,-22.149,-22.149,-22.149,-22.149,-22.149,-22.149 +09.00,-22.125,-22.125,-22.125,-22.125,-22.125,-22.125,-22.125 diff --git a/data/radiation/abund_25_rad_loss.dat b/data/radiation/abund_25_rad_loss.dat new file mode 100644 index 0000000..6feeaec --- /dev/null +++ b/data/radiation/abund_25_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.152,-23.131,-23.126,-23.122,-23.123,-23.125,-23.125 +04.05,-22.708,-22.674,-22.667,-22.661,-22.660,-22.661,-22.662 +04.10,-22.296,-22.256,-22.247,-22.240,-22.240,-22.240,-22.241 +04.15,-21.958,-21.917,-21.909,-21.902,-21.901,-21.902,-21.902 +04.20,-21.769,-21.730,-21.723,-21.717,-21.717,-21.717,-21.717 +04.25,-21.746,-21.712,-21.706,-21.702,-21.702,-21.703,-21.704 +04.30,-21.815,-21.788,-21.783,-21.782,-21.782,-21.784,-21.785 +04.35,-21.903,-21.882,-21.879,-21.881,-21.882,-21.886,-21.888 +04.40,-21.965,-21.950,-21.949,-21.953,-21.956,-21.962,-21.965 +04.45,-21.982,-21.972,-21.973,-21.979,-21.982,-21.990,-21.994 +04.50,-21.966,-21.961,-21.962,-21.969,-21.973,-21.982,-21.986 +04.55,-21.933,-21.932,-21.933,-21.941,-21.945,-21.954,-21.959 +04.60,-21.887,-21.890,-21.892,-21.901,-21.905,-21.915,-21.920 +04.65,-21.827,-21.832,-21.835,-21.846,-21.851,-21.862,-21.867 +04.70,-21.748,-21.754,-21.758,-21.770,-21.776,-21.786,-21.791 +04.75,-21.656,-21.662,-21.666,-21.679,-21.684,-21.693,-21.697 +04.80,-21.561,-21.567,-21.571,-21.582,-21.586,-21.592,-21.595 +04.85,-21.480,-21.486,-21.489,-21.499,-21.501,-21.504,-21.506 +04.90,-21.433,-21.438,-21.441,-21.450,-21.452,-21.454,-21.455 +04.95,-21.407,-21.413,-21.416,-21.424,-21.426,-21.429,-21.430 +05.00,-21.393,-21.398,-21.401,-21.408,-21.410,-21.415,-21.417 +05.05,-21.405,-21.409,-21.411,-21.416,-21.418,-21.424,-21.428 +05.10,-21.423,-21.426,-21.427,-21.431,-21.433,-21.442,-21.447 +05.15,-21.417,-21.420,-21.421,-21.424,-21.426,-21.438,-21.445 +05.20,-21.394,-21.397,-21.398,-21.401,-21.403,-21.416,-21.424 +05.25,-21.373,-21.376,-21.376,-21.379,-21.382,-21.395,-21.404 +05.30,-21.363,-21.365,-21.366,-21.368,-21.371,-21.383,-21.393 +05.35,-21.355,-21.357,-21.358,-21.360,-21.362,-21.375,-21.384 +05.40,-21.357,-21.359,-21.360,-21.362,-21.364,-21.374,-21.383 +05.45,-21.400,-21.402,-21.402,-21.404,-21.406,-21.414,-21.421 +05.50,-21.494,-21.495,-21.496,-21.498,-21.499,-21.505,-21.509 +05.55,-21.579,-21.581,-21.582,-21.584,-21.585,-21.588,-21.590 +05.60,-21.611,-21.613,-21.614,-21.616,-21.617,-21.619,-21.620 +05.65,-21.602,-21.604,-21.605,-21.607,-21.608,-21.610,-21.611 +05.70,-21.578,-21.580,-21.580,-21.583,-21.584,-21.586,-21.587 +05.75,-21.557,-21.559,-21.559,-21.562,-21.563,-21.565,-21.566 +05.80,-21.545,-21.546,-21.547,-21.551,-21.552,-21.554,-21.555 +05.85,-21.533,-21.535,-21.536,-21.540,-21.542,-21.545,-21.545 +05.90,-21.519,-21.521,-21.522,-21.527,-21.529,-21.532,-21.533 +05.95,-21.507,-21.510,-21.511,-21.516,-21.518,-21.521,-21.522 +06.00,-21.506,-21.509,-21.510,-21.515,-21.517,-21.521,-21.522 +06.05,-21.519,-21.521,-21.523,-21.527,-21.529,-21.533,-21.534 +06.10,-21.539,-21.541,-21.542,-21.546,-21.547,-21.551,-21.552 +06.15,-21.563,-21.564,-21.566,-21.568,-21.570,-21.572,-21.573 +06.20,-21.596,-21.597,-21.598,-21.601,-21.602,-21.604,-21.604 +06.25,-21.649,-21.650,-21.651,-21.653,-21.654,-21.656,-21.656 +06.30,-21.734,-21.735,-21.735,-21.737,-21.738,-21.740,-21.740 +06.35,-21.851,-21.852,-21.852,-21.854,-21.855,-21.857,-21.857 +06.40,-21.977,-21.977,-21.978,-21.979,-21.980,-21.982,-21.982 +06.45,-22.087,-22.087,-22.087,-22.089,-22.089,-22.090,-22.091 +06.50,-22.170,-22.170,-22.171,-22.171,-22.172,-22.173,-22.173 +06.55,-22.225,-22.225,-22.225,-22.225,-22.226,-22.226,-22.226 +06.60,-22.253,-22.253,-22.253,-22.253,-22.253,-22.254,-22.254 +06.65,-22.262,-22.262,-22.262,-22.262,-22.262,-22.262,-22.262 +06.70,-22.257,-22.257,-22.257,-22.257,-22.258,-22.258,-22.258 +06.75,-22.246,-22.246,-22.246,-22.246,-22.246,-22.246,-22.246 +06.80,-22.232,-22.232,-22.232,-22.232,-22.232,-22.233,-22.233 +06.85,-22.220,-22.220,-22.220,-22.220,-22.220,-22.220,-22.220 +06.90,-22.213,-22.213,-22.213,-22.213,-22.213,-22.213,-22.213 +06.95,-22.215,-22.215,-22.215,-22.215,-22.215,-22.215,-22.215 +07.00,-22.231,-22.231,-22.231,-22.231,-22.231,-22.231,-22.231 +07.05,-22.264,-22.264,-22.264,-22.264,-22.264,-22.264,-22.264 +07.10,-22.318,-22.318,-22.318,-22.318,-22.318,-22.318,-22.318 +07.15,-22.390,-22.390,-22.390,-22.390,-22.390,-22.390,-22.390 +07.20,-22.466,-22.466,-22.466,-22.466,-22.466,-22.466,-22.466 +07.25,-22.533,-22.533,-22.533,-22.533,-22.533,-22.533,-22.533 +07.30,-22.585,-22.585,-22.585,-22.585,-22.585,-22.585,-22.585 +07.35,-22.621,-22.621,-22.621,-22.621,-22.621,-22.621,-22.621 +07.40,-22.644,-22.644,-22.644,-22.644,-22.644,-22.644,-22.644 +07.45,-22.656,-22.656,-22.656,-22.656,-22.656,-22.656,-22.656 +07.50,-22.660,-22.660,-22.660,-22.660,-22.660,-22.660,-22.660 +07.55,-22.658,-22.658,-22.658,-22.658,-22.658,-22.658,-22.658 +07.60,-22.652,-22.652,-22.652,-22.652,-22.652,-22.652,-22.652 +07.65,-22.642,-22.642,-22.642,-22.642,-22.642,-22.642,-22.642 +07.70,-22.630,-22.630,-22.630,-22.630,-22.630,-22.630,-22.630 +07.75,-22.616,-22.616,-22.616,-22.616,-22.616,-22.616,-22.616 +07.80,-22.601,-22.601,-22.601,-22.601,-22.601,-22.601,-22.601 +07.85,-22.586,-22.586,-22.586,-22.586,-22.586,-22.586,-22.586 +07.90,-22.571,-22.571,-22.571,-22.571,-22.571,-22.571,-22.570 +07.95,-22.555,-22.555,-22.555,-22.555,-22.555,-22.555,-22.555 +08.00,-22.538,-22.538,-22.538,-22.538,-22.538,-22.538,-22.538 +08.05,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522 +08.10,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505 +08.15,-22.487,-22.487,-22.487,-22.487,-22.487,-22.487,-22.487 +08.20,-22.469,-22.469,-22.469,-22.469,-22.469,-22.469,-22.469 +08.25,-22.451,-22.451,-22.451,-22.451,-22.451,-22.451,-22.451 +08.30,-22.431,-22.431,-22.431,-22.431,-22.431,-22.431,-22.431 +08.35,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411 +08.40,-22.391,-22.391,-22.391,-22.391,-22.391,-22.391,-22.391 +08.45,-22.370,-22.370,-22.370,-22.370,-22.370,-22.370,-22.370 +08.50,-22.348,-22.348,-22.348,-22.348,-22.348,-22.348,-22.348 +08.55,-22.327,-22.327,-22.327,-22.327,-22.327,-22.327,-22.327 +08.60,-22.304,-22.304,-22.304,-22.304,-22.304,-22.304,-22.304 +08.65,-22.282,-22.282,-22.282,-22.282,-22.282,-22.282,-22.282 +08.70,-22.259,-22.259,-22.259,-22.259,-22.259,-22.259,-22.259 +08.75,-22.236,-22.236,-22.236,-22.236,-22.236,-22.236,-22.236 +08.80,-22.213,-22.213,-22.213,-22.213,-22.213,-22.213,-22.213 +08.85,-22.190,-22.190,-22.190,-22.190,-22.190,-22.190,-22.190 +08.90,-22.166,-22.166,-22.166,-22.166,-22.166,-22.166,-22.166 +08.95,-22.142,-22.142,-22.142,-22.142,-22.142,-22.142,-22.142 +09.00,-22.119,-22.119,-22.119,-22.119,-22.119,-22.119,-22.119 diff --git a/data/radiation/abund_30_rad_loss.dat b/data/radiation/abund_30_rad_loss.dat new file mode 100644 index 0000000..71ada48 --- /dev/null +++ b/data/radiation/abund_30_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.110,-23.092,-23.087,-23.084,-23.085,-23.087,-23.088 +04.05,-22.677,-22.647,-22.639,-22.634,-22.634,-22.635,-22.636 +04.10,-22.275,-22.237,-22.229,-22.222,-22.222,-22.223,-22.223 +04.15,-21.943,-21.903,-21.895,-21.889,-21.888,-21.889,-21.889 +04.20,-21.756,-21.718,-21.711,-21.706,-21.705,-21.706,-21.706 +04.25,-21.732,-21.699,-21.693,-21.690,-21.690,-21.691,-21.691 +04.30,-21.800,-21.773,-21.769,-21.767,-21.768,-21.771,-21.772 +04.35,-21.885,-21.865,-21.862,-21.864,-21.866,-21.870,-21.872 +04.40,-21.944,-21.930,-21.929,-21.933,-21.936,-21.943,-21.946 +04.45,-21.956,-21.948,-21.948,-21.954,-21.958,-21.967,-21.971 +04.50,-21.936,-21.932,-21.933,-21.940,-21.943,-21.953,-21.958 +04.55,-21.900,-21.899,-21.901,-21.908,-21.911,-21.921,-21.926 +04.60,-21.854,-21.856,-21.858,-21.866,-21.870,-21.880,-21.886 +04.65,-21.794,-21.799,-21.802,-21.812,-21.817,-21.827,-21.833 +04.70,-21.717,-21.723,-21.726,-21.738,-21.743,-21.754,-21.759 +04.75,-21.627,-21.633,-21.637,-21.649,-21.653,-21.662,-21.667 +04.80,-21.535,-21.541,-21.544,-21.555,-21.559,-21.565,-21.568 +04.85,-21.459,-21.464,-21.468,-21.477,-21.480,-21.483,-21.485 +04.90,-21.419,-21.424,-21.428,-21.436,-21.438,-21.440,-21.441 +04.95,-21.400,-21.406,-21.409,-21.417,-21.419,-21.422,-21.423 +05.00,-21.390,-21.395,-21.397,-21.405,-21.407,-21.411,-21.414 +05.05,-21.402,-21.406,-21.408,-21.414,-21.416,-21.422,-21.426 +05.10,-21.420,-21.424,-21.425,-21.429,-21.431,-21.440,-21.445 +05.15,-21.414,-21.417,-21.418,-21.421,-21.424,-21.435,-21.442 +05.20,-21.391,-21.393,-21.394,-21.397,-21.400,-21.413,-21.421 +05.25,-21.369,-21.372,-21.372,-21.375,-21.378,-21.391,-21.400 +05.30,-21.357,-21.360,-21.360,-21.363,-21.365,-21.378,-21.387 +05.35,-21.348,-21.350,-21.351,-21.354,-21.356,-21.368,-21.377 +05.40,-21.348,-21.350,-21.351,-21.353,-21.355,-21.365,-21.374 +05.45,-21.386,-21.388,-21.389,-21.391,-21.393,-21.401,-21.408 +05.50,-21.472,-21.474,-21.475,-21.477,-21.478,-21.484,-21.488 +05.55,-21.546,-21.548,-21.549,-21.551,-21.552,-21.555,-21.557 +05.60,-21.567,-21.569,-21.570,-21.573,-21.574,-21.576,-21.577 +05.65,-21.550,-21.552,-21.553,-21.556,-21.557,-21.559,-21.559 +05.70,-21.521,-21.522,-21.523,-21.526,-21.527,-21.529,-21.530 +05.75,-21.494,-21.496,-21.497,-21.500,-21.501,-21.503,-21.504 +05.80,-21.477,-21.479,-21.480,-21.483,-21.484,-21.487,-21.487 +05.85,-21.461,-21.463,-21.464,-21.468,-21.470,-21.473,-21.473 +05.90,-21.444,-21.447,-21.448,-21.452,-21.454,-21.458,-21.458 +05.95,-21.432,-21.434,-21.436,-21.440,-21.442,-21.446,-21.447 +06.00,-21.431,-21.433,-21.435,-21.439,-21.441,-21.445,-21.446 +06.05,-21.443,-21.446,-21.447,-21.451,-21.453,-21.457,-21.458 +06.10,-21.464,-21.465,-21.467,-21.470,-21.472,-21.475,-21.476 +06.15,-21.488,-21.489,-21.490,-21.493,-21.495,-21.497,-21.498 +06.20,-21.522,-21.523,-21.524,-21.526,-21.527,-21.529,-21.530 +06.25,-21.576,-21.577,-21.577,-21.580,-21.581,-21.582,-21.583 +06.30,-21.662,-21.663,-21.664,-21.666,-21.667,-21.668,-21.669 +06.35,-21.782,-21.783,-21.784,-21.786,-21.787,-21.788,-21.788 +06.40,-21.911,-21.912,-21.912,-21.914,-21.915,-21.916,-21.916 +06.45,-22.024,-22.024,-22.024,-22.025,-22.026,-22.027,-22.028 +06.50,-22.108,-22.109,-22.109,-22.110,-22.110,-22.111,-22.111 +06.55,-22.163,-22.163,-22.163,-22.164,-22.164,-22.165,-22.165 +06.60,-22.191,-22.191,-22.191,-22.191,-22.191,-22.192,-22.192 +06.65,-22.198,-22.198,-22.198,-22.199,-22.199,-22.199,-22.199 +06.70,-22.193,-22.193,-22.193,-22.193,-22.193,-22.193,-22.193 +06.75,-22.180,-22.180,-22.180,-22.180,-22.180,-22.181,-22.181 +06.80,-22.166,-22.166,-22.166,-22.166,-22.166,-22.166,-22.166 +06.85,-22.152,-22.152,-22.152,-22.152,-22.152,-22.152,-22.152 +06.90,-22.145,-22.145,-22.145,-22.145,-22.145,-22.145,-22.145 +06.95,-22.147,-22.147,-22.147,-22.148,-22.148,-22.148,-22.148 +07.00,-22.164,-22.164,-22.164,-22.164,-22.164,-22.164,-22.164 +07.05,-22.198,-22.198,-22.198,-22.198,-22.198,-22.198,-22.198 +07.10,-22.255,-22.255,-22.255,-22.255,-22.255,-22.255,-22.255 +07.15,-22.330,-22.330,-22.330,-22.330,-22.330,-22.330,-22.330 +07.20,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411,-22.411 +07.25,-22.484,-22.484,-22.484,-22.484,-22.484,-22.484,-22.484 +07.30,-22.541,-22.541,-22.541,-22.541,-22.541,-22.541,-22.541 +07.35,-22.582,-22.582,-22.582,-22.582,-22.582,-22.582,-22.582 +07.40,-22.609,-22.609,-22.609,-22.609,-22.609,-22.609,-22.609 +07.45,-22.624,-22.624,-22.624,-22.624,-22.624,-22.624,-22.624 +07.50,-22.631,-22.631,-22.631,-22.631,-22.631,-22.631,-22.631 +07.55,-22.631,-22.631,-22.631,-22.631,-22.631,-22.631,-22.631 +07.60,-22.626,-22.626,-22.626,-22.626,-22.626,-22.626,-22.626 +07.65,-22.618,-22.618,-22.618,-22.618,-22.618,-22.618,-22.618 +07.70,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607 +07.75,-22.594,-22.594,-22.594,-22.594,-22.594,-22.594,-22.594 +07.80,-22.581,-22.581,-22.581,-22.581,-22.581,-22.581,-22.581 +07.85,-22.566,-22.566,-22.566,-22.566,-22.566,-22.566,-22.566 +07.90,-22.552,-22.552,-22.552,-22.552,-22.552,-22.552,-22.552 +07.95,-22.537,-22.537,-22.537,-22.537,-22.537,-22.537,-22.537 +08.00,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522 +08.05,-22.506,-22.506,-22.506,-22.506,-22.506,-22.506,-22.506 +08.10,-22.490,-22.490,-22.490,-22.490,-22.490,-22.490,-22.490 +08.15,-22.474,-22.474,-22.474,-22.474,-22.474,-22.474,-22.474 +08.20,-22.457,-22.457,-22.457,-22.457,-22.457,-22.457,-22.457 +08.25,-22.439,-22.439,-22.439,-22.439,-22.439,-22.439,-22.439 +08.30,-22.420,-22.420,-22.420,-22.420,-22.420,-22.420,-22.420 +08.35,-22.401,-22.401,-22.401,-22.401,-22.401,-22.401,-22.401 +08.40,-22.381,-22.381,-22.381,-22.381,-22.381,-22.381,-22.381 +08.45,-22.360,-22.360,-22.360,-22.360,-22.360,-22.360,-22.360 +08.50,-22.339,-22.339,-22.339,-22.339,-22.339,-22.339,-22.339 +08.55,-22.318,-22.318,-22.318,-22.318,-22.318,-22.318,-22.318 +08.60,-22.296,-22.296,-22.296,-22.296,-22.296,-22.296,-22.296 +08.65,-22.274,-22.274,-22.274,-22.274,-22.274,-22.274,-22.274 +08.70,-22.252,-22.252,-22.252,-22.252,-22.252,-22.252,-22.252 +08.75,-22.229,-22.229,-22.229,-22.229,-22.229,-22.229,-22.229 +08.80,-22.206,-22.206,-22.206,-22.206,-22.206,-22.206,-22.206 +08.85,-22.183,-22.183,-22.183,-22.183,-22.183,-22.183,-22.183 +08.90,-22.159,-22.159,-22.159,-22.159,-22.159,-22.159,-22.159 +08.95,-22.136,-22.136,-22.136,-22.136,-22.136,-22.136,-22.136 +09.00,-22.112,-22.112,-22.112,-22.112,-22.112,-22.112,-22.112 diff --git a/data/radiation/abund_35_rad_loss.dat b/data/radiation/abund_35_rad_loss.dat new file mode 100644 index 0000000..9172db8 --- /dev/null +++ b/data/radiation/abund_35_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.076,-23.060,-23.056,-23.054,-23.054,-23.057,-23.058 +04.05,-22.652,-22.623,-22.617,-22.612,-22.612,-22.613,-22.614 +04.10,-22.257,-22.221,-22.213,-22.207,-22.207,-22.207,-22.208 +04.15,-21.931,-21.892,-21.884,-21.878,-21.877,-21.878,-21.878 +04.20,-21.744,-21.708,-21.701,-21.696,-21.695,-21.696,-21.696 +04.25,-21.720,-21.688,-21.682,-21.679,-21.679,-21.680,-21.681 +04.30,-21.786,-21.760,-21.756,-21.755,-21.756,-21.759,-21.760 +04.35,-21.870,-21.851,-21.848,-21.850,-21.852,-21.857,-21.859 +04.40,-21.926,-21.913,-21.912,-21.917,-21.920,-21.927,-21.930 +04.45,-21.935,-21.927,-21.928,-21.934,-21.938,-21.947,-21.951 +04.50,-21.912,-21.908,-21.909,-21.915,-21.919,-21.929,-21.934 +04.55,-21.873,-21.873,-21.874,-21.880,-21.884,-21.893,-21.899 +04.60,-21.826,-21.828,-21.830,-21.838,-21.841,-21.852,-21.858 +04.65,-21.767,-21.772,-21.775,-21.784,-21.788,-21.799,-21.805 +04.70,-21.691,-21.697,-21.700,-21.711,-21.716,-21.727,-21.732 +04.75,-21.603,-21.609,-21.612,-21.624,-21.628,-21.637,-21.642 +04.80,-21.513,-21.519,-21.522,-21.533,-21.536,-21.543,-21.546 +04.85,-21.441,-21.446,-21.450,-21.459,-21.461,-21.465,-21.467 +04.90,-21.407,-21.412,-21.416,-21.424,-21.426,-21.428,-21.429 +04.95,-21.394,-21.399,-21.403,-21.411,-21.413,-21.415,-21.417 +05.00,-21.386,-21.392,-21.394,-21.402,-21.404,-21.408,-21.411 +05.05,-21.400,-21.404,-21.406,-21.411,-21.414,-21.420,-21.424 +05.10,-21.418,-21.422,-21.423,-21.427,-21.429,-21.438,-21.443 +05.15,-21.411,-21.415,-21.416,-21.419,-21.421,-21.432,-21.439 +05.20,-21.388,-21.390,-21.391,-21.395,-21.397,-21.410,-21.418 +05.25,-21.365,-21.368,-21.369,-21.372,-21.374,-21.387,-21.396 +05.30,-21.353,-21.355,-21.356,-21.359,-21.361,-21.373,-21.383 +05.35,-21.342,-21.345,-21.345,-21.348,-21.350,-21.362,-21.371 +05.40,-21.340,-21.342,-21.343,-21.345,-21.347,-21.357,-21.366 +05.45,-21.375,-21.377,-21.378,-21.380,-21.382,-21.390,-21.396 +05.50,-21.454,-21.456,-21.457,-21.459,-21.460,-21.466,-21.470 +05.55,-21.519,-21.521,-21.522,-21.524,-21.525,-21.528,-21.530 +05.60,-21.532,-21.534,-21.535,-21.538,-21.539,-21.541,-21.542 +05.65,-21.510,-21.511,-21.512,-21.515,-21.516,-21.518,-21.519 +05.70,-21.476,-21.477,-21.478,-21.481,-21.482,-21.484,-21.485 +05.75,-21.446,-21.448,-21.448,-21.451,-21.453,-21.455,-21.456 +05.80,-21.425,-21.426,-21.428,-21.431,-21.432,-21.435,-21.436 +05.85,-21.406,-21.408,-21.410,-21.413,-21.415,-21.418,-21.419 +05.90,-21.388,-21.390,-21.392,-21.396,-21.398,-21.401,-21.402 +05.95,-21.375,-21.377,-21.378,-21.383,-21.385,-21.389,-21.390 +06.00,-21.373,-21.376,-21.377,-21.382,-21.384,-21.388,-21.389 +06.05,-21.386,-21.388,-21.390,-21.394,-21.396,-21.400,-21.401 +06.10,-21.406,-21.408,-21.409,-21.413,-21.415,-21.418,-21.419 +06.15,-21.431,-21.432,-21.433,-21.437,-21.438,-21.440,-21.441 +06.20,-21.465,-21.466,-21.467,-21.470,-21.471,-21.473,-21.474 +06.25,-21.520,-21.521,-21.522,-21.524,-21.525,-21.527,-21.527 +06.30,-21.607,-21.608,-21.609,-21.611,-21.612,-21.614,-21.614 +06.35,-21.729,-21.730,-21.731,-21.733,-21.734,-21.735,-21.736 +06.40,-21.861,-21.861,-21.862,-21.863,-21.864,-21.866,-21.866 +06.45,-21.975,-21.975,-21.975,-21.977,-21.977,-21.979,-21.979 +06.50,-22.061,-22.061,-22.061,-22.062,-22.062,-22.063,-22.063 +06.55,-22.115,-22.115,-22.115,-22.116,-22.116,-22.117,-22.117 +06.60,-22.143,-22.143,-22.143,-22.143,-22.143,-22.144,-22.144 +06.65,-22.149,-22.149,-22.149,-22.150,-22.150,-22.150,-22.150 +06.70,-22.143,-22.143,-22.143,-22.143,-22.143,-22.143,-22.143 +06.75,-22.130,-22.130,-22.130,-22.130,-22.130,-22.130,-22.130 +06.80,-22.114,-22.114,-22.114,-22.114,-22.114,-22.114,-22.114 +06.85,-22.100,-22.100,-22.100,-22.101,-22.101,-22.101,-22.101 +06.90,-22.093,-22.093,-22.093,-22.093,-22.093,-22.093,-22.093 +06.95,-22.095,-22.095,-22.095,-22.096,-22.096,-22.096,-22.096 +07.00,-22.112,-22.112,-22.112,-22.112,-22.112,-22.112,-22.112 +07.05,-22.147,-22.148,-22.148,-22.148,-22.148,-22.148,-22.148 +07.10,-22.206,-22.206,-22.206,-22.206,-22.206,-22.206,-22.206 +07.15,-22.284,-22.284,-22.284,-22.284,-22.284,-22.284,-22.284 +07.20,-22.368,-22.368,-22.368,-22.368,-22.368,-22.368,-22.368 +07.25,-22.444,-22.444,-22.444,-22.444,-22.444,-22.444,-22.444 +07.30,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505 +07.35,-22.550,-22.550,-22.550,-22.550,-22.550,-22.550,-22.550 +07.40,-22.580,-22.580,-22.580,-22.580,-22.580,-22.580,-22.580 +07.45,-22.598,-22.598,-22.598,-22.598,-22.598,-22.598,-22.598 +07.50,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607,-22.607 +07.55,-22.608,-22.608,-22.608,-22.608,-22.608,-22.608,-22.608 +07.60,-22.605,-22.605,-22.605,-22.605,-22.605,-22.605,-22.605 +07.65,-22.598,-22.598,-22.598,-22.598,-22.598,-22.598,-22.598 +07.70,-22.588,-22.588,-22.588,-22.588,-22.588,-22.588,-22.588 +07.75,-22.576,-22.576,-22.576,-22.576,-22.576,-22.576,-22.576 +07.80,-22.563,-22.563,-22.563,-22.563,-22.563,-22.563,-22.563 +07.85,-22.550,-22.550,-22.550,-22.550,-22.550,-22.550,-22.550 +07.90,-22.536,-22.536,-22.536,-22.536,-22.536,-22.536,-22.536 +07.95,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522,-22.522 +08.00,-22.508,-22.508,-22.508,-22.508,-22.508,-22.508,-22.508 +08.05,-22.493,-22.493,-22.493,-22.493,-22.493,-22.493,-22.493 +08.10,-22.478,-22.478,-22.478,-22.478,-22.478,-22.478,-22.478 +08.15,-22.462,-22.462,-22.462,-22.462,-22.462,-22.462,-22.462 +08.20,-22.446,-22.446,-22.446,-22.446,-22.446,-22.446,-22.446 +08.25,-22.428,-22.428,-22.428,-22.428,-22.428,-22.428,-22.428 +08.30,-22.410,-22.410,-22.410,-22.410,-22.410,-22.410,-22.410 +08.35,-22.392,-22.392,-22.392,-22.392,-22.392,-22.392,-22.392 +08.40,-22.372,-22.372,-22.372,-22.372,-22.372,-22.372,-22.372 +08.45,-22.352,-22.352,-22.352,-22.352,-22.352,-22.352,-22.352 +08.50,-22.332,-22.332,-22.332,-22.332,-22.332,-22.332,-22.332 +08.55,-22.311,-22.311,-22.311,-22.311,-22.311,-22.311,-22.311 +08.60,-22.289,-22.289,-22.289,-22.289,-22.289,-22.289,-22.289 +08.65,-22.267,-22.267,-22.267,-22.267,-22.267,-22.267,-22.267 +08.70,-22.245,-22.245,-22.245,-22.245,-22.245,-22.245,-22.245 +08.75,-22.223,-22.223,-22.223,-22.223,-22.223,-22.223,-22.223 +08.80,-22.200,-22.200,-22.200,-22.200,-22.200,-22.200,-22.200 +08.85,-22.177,-22.177,-22.177,-22.177,-22.177,-22.177,-22.177 +08.90,-22.154,-22.154,-22.154,-22.154,-22.154,-22.154,-22.154 +08.95,-22.130,-22.130,-22.130,-22.130,-22.130,-22.130,-22.130 +09.00,-22.107,-22.107,-22.107,-22.107,-22.107,-22.107,-22.107 diff --git a/data/radiation/abund_40_rad_loss.dat b/data/radiation/abund_40_rad_loss.dat new file mode 100644 index 0000000..1b0a3a5 --- /dev/null +++ b/data/radiation/abund_40_rad_loss.dat @@ -0,0 +1,102 @@ +Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 +04.00,-23.040,-23.026,-23.023,-23.021,-23.021,-23.024,-23.025 +04.05,-22.625,-22.598,-22.592,-22.588,-22.588,-22.589,-22.590 +04.10,-22.238,-22.203,-22.196,-22.190,-22.190,-22.191,-22.191 +04.15,-21.916,-21.879,-21.871,-21.865,-21.865,-21.866,-21.866 +04.20,-21.732,-21.696,-21.690,-21.685,-21.684,-21.685,-21.685 +04.25,-21.706,-21.676,-21.670,-21.667,-21.667,-21.669,-21.669 +04.30,-21.771,-21.746,-21.743,-21.742,-21.743,-21.746,-21.747 +04.35,-21.853,-21.835,-21.832,-21.835,-21.837,-21.842,-21.844 +04.40,-21.906,-21.894,-21.893,-21.898,-21.901,-21.909,-21.913 +04.45,-21.912,-21.905,-21.905,-21.912,-21.915,-21.925,-21.930 +04.50,-21.885,-21.882,-21.883,-21.889,-21.893,-21.903,-21.908 +04.55,-21.844,-21.844,-21.845,-21.851,-21.854,-21.864,-21.870 +04.60,-21.796,-21.799,-21.801,-21.807,-21.811,-21.821,-21.827 +04.65,-21.738,-21.743,-21.745,-21.754,-21.758,-21.769,-21.776 +04.70,-21.663,-21.669,-21.672,-21.682,-21.687,-21.698,-21.704 +04.75,-21.576,-21.582,-21.586,-21.597,-21.601,-21.610,-21.615 +04.80,-21.489,-21.495,-21.498,-21.508,-21.512,-21.518,-21.522 +04.85,-21.421,-21.427,-21.430,-21.439,-21.441,-21.445,-21.447 +04.90,-21.393,-21.399,-21.402,-21.411,-21.413,-21.414,-21.416 +04.95,-21.387,-21.393,-21.396,-21.404,-21.406,-21.409,-21.410 +05.00,-21.383,-21.388,-21.391,-21.398,-21.400,-21.405,-21.407 +05.05,-21.397,-21.402,-21.404,-21.409,-21.411,-21.418,-21.421 +05.10,-21.416,-21.419,-21.421,-21.424,-21.427,-21.435,-21.441 +05.15,-21.408,-21.412,-21.413,-21.416,-21.419,-21.430,-21.437 +05.20,-21.384,-21.387,-21.388,-21.391,-21.394,-21.406,-21.414 +05.25,-21.361,-21.364,-21.365,-21.368,-21.370,-21.383,-21.392 +05.30,-21.347,-21.350,-21.351,-21.354,-21.356,-21.368,-21.377 +05.35,-21.335,-21.338,-21.339,-21.341,-21.343,-21.355,-21.364 +05.40,-21.331,-21.333,-21.334,-21.337,-21.338,-21.349,-21.357 +05.45,-21.362,-21.364,-21.365,-21.368,-21.369,-21.377,-21.383 +05.50,-21.434,-21.436,-21.437,-21.440,-21.441,-21.446,-21.450 +05.55,-21.489,-21.492,-21.493,-21.495,-21.496,-21.499,-21.501 +05.60,-21.495,-21.497,-21.498,-21.501,-21.502,-21.504,-21.505 +05.65,-21.467,-21.469,-21.470,-21.473,-21.474,-21.476,-21.476 +05.70,-21.429,-21.431,-21.432,-21.435,-21.436,-21.438,-21.439 +05.75,-21.396,-21.398,-21.399,-21.402,-21.403,-21.405,-21.406 +05.80,-21.372,-21.374,-21.375,-21.378,-21.380,-21.382,-21.383 +05.85,-21.351,-21.353,-21.354,-21.358,-21.360,-21.363,-21.364 +05.90,-21.331,-21.333,-21.335,-21.339,-21.341,-21.345,-21.345 +05.95,-21.317,-21.319,-21.321,-21.326,-21.328,-21.332,-21.333 +06.00,-21.315,-21.318,-21.320,-21.324,-21.326,-21.330,-21.332 +06.05,-21.328,-21.331,-21.332,-21.336,-21.338,-21.342,-21.343 +06.10,-21.349,-21.351,-21.352,-21.356,-21.357,-21.361,-21.362 +06.15,-21.374,-21.375,-21.376,-21.379,-21.381,-21.383,-21.384 +06.20,-21.408,-21.410,-21.410,-21.413,-21.414,-21.416,-21.417 +06.25,-21.463,-21.465,-21.465,-21.468,-21.469,-21.470,-21.471 +06.30,-21.552,-21.553,-21.554,-21.556,-21.557,-21.559,-21.559 +06.35,-21.676,-21.677,-21.677,-21.679,-21.680,-21.682,-21.682 +06.40,-21.809,-21.810,-21.810,-21.812,-21.813,-21.814,-21.814 +06.45,-21.925,-21.925,-21.925,-21.927,-21.927,-21.929,-21.929 +06.50,-22.011,-22.011,-22.012,-22.012,-22.013,-22.014,-22.014 +06.55,-22.066,-22.066,-22.066,-22.067,-22.067,-22.068,-22.068 +06.60,-22.093,-22.093,-22.093,-22.093,-22.094,-22.094,-22.094 +06.65,-22.099,-22.099,-22.099,-22.099,-22.099,-22.100,-22.100 +06.70,-22.092,-22.092,-22.092,-22.092,-22.092,-22.092,-22.092 +06.75,-22.078,-22.078,-22.078,-22.078,-22.078,-22.078,-22.078 +06.80,-22.062,-22.062,-22.062,-22.062,-22.062,-22.062,-22.062 +06.85,-22.048,-22.048,-22.048,-22.048,-22.048,-22.048,-22.048 +06.90,-22.040,-22.040,-22.040,-22.040,-22.040,-22.040,-22.040 +06.95,-22.042,-22.042,-22.043,-22.043,-22.043,-22.043,-22.043 +07.00,-22.059,-22.060,-22.060,-22.060,-22.060,-22.060,-22.060 +07.05,-22.096,-22.096,-22.096,-22.096,-22.096,-22.096,-22.096 +07.10,-22.155,-22.155,-22.155,-22.156,-22.156,-22.156,-22.156 +07.15,-22.236,-22.236,-22.236,-22.236,-22.236,-22.236,-22.236 +07.20,-22.323,-22.323,-22.323,-22.323,-22.323,-22.323,-22.323 +07.25,-22.403,-22.403,-22.403,-22.403,-22.403,-22.403,-22.403 +07.30,-22.468,-22.468,-22.468,-22.468,-22.468,-22.468,-22.468 +07.35,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516,-22.516 +07.40,-22.549,-22.549,-22.549,-22.549,-22.549,-22.549,-22.549 +07.45,-22.569,-22.569,-22.569,-22.569,-22.569,-22.569,-22.569 +07.50,-22.580,-22.580,-22.580,-22.580,-22.580,-22.580,-22.580 +07.55,-22.584,-22.584,-22.584,-22.584,-22.584,-22.584,-22.584 +07.60,-22.582,-22.582,-22.582,-22.582,-22.582,-22.582,-22.582 +07.65,-22.576,-22.576,-22.576,-22.576,-22.576,-22.576,-22.576 +07.70,-22.567,-22.567,-22.567,-22.567,-22.567,-22.567,-22.567 +07.75,-22.556,-22.556,-22.556,-22.556,-22.556,-22.556,-22.556 +07.80,-22.544,-22.544,-22.544,-22.544,-22.544,-22.544,-22.544 +07.85,-22.531,-22.531,-22.531,-22.531,-22.531,-22.531,-22.531 +07.90,-22.518,-22.518,-22.518,-22.518,-22.518,-22.518,-22.518 +07.95,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505,-22.505 +08.00,-22.492,-22.492,-22.492,-22.492,-22.492,-22.492,-22.492 +08.05,-22.478,-22.478,-22.478,-22.478,-22.478,-22.478,-22.478 +08.10,-22.464,-22.464,-22.464,-22.464,-22.464,-22.464,-22.464 +08.15,-22.449,-22.449,-22.449,-22.449,-22.449,-22.449,-22.449 +08.20,-22.433,-22.433,-22.433,-22.433,-22.433,-22.433,-22.433 +08.25,-22.417,-22.417,-22.417,-22.417,-22.417,-22.417,-22.417 +08.30,-22.400,-22.400,-22.400,-22.400,-22.400,-22.400,-22.400 +08.35,-22.381,-22.381,-22.381,-22.381,-22.381,-22.381,-22.381 +08.40,-22.363,-22.363,-22.363,-22.363,-22.363,-22.363,-22.363 +08.45,-22.343,-22.343,-22.343,-22.343,-22.343,-22.343,-22.343 +08.50,-22.323,-22.323,-22.323,-22.323,-22.323,-22.323,-22.323 +08.55,-22.302,-22.302,-22.302,-22.302,-22.302,-22.302,-22.302 +08.60,-22.281,-22.281,-22.281,-22.281,-22.281,-22.281,-22.281 +08.65,-22.260,-22.260,-22.260,-22.260,-22.260,-22.260,-22.260 +08.70,-22.238,-22.238,-22.238,-22.238,-22.238,-22.238,-22.238 +08.75,-22.215,-22.215,-22.215,-22.215,-22.215,-22.215,-22.215 +08.80,-22.193,-22.193,-22.193,-22.193,-22.193,-22.193,-22.193 +08.85,-22.170,-22.170,-22.170,-22.170,-22.170,-22.170,-22.170 +08.90,-22.147,-22.147,-22.147,-22.147,-22.147,-22.147,-22.147 +08.95,-22.124,-22.124,-22.124,-22.124,-22.124,-22.124,-22.124 +09.00,-22.100,-22.100,-22.100,-22.100,-22.100,-22.100,-22.100 From b5019ab583d6916ff397626a8f3dcb34c3c2b739 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 18 Dec 2023 14:54:21 -1000 Subject: [PATCH 05/24] add function to read abundance tables and calculate new losses --- source/loop.cpp | 48 +++++++++++++++++++++++++++++++++++++++------- source/loop.h | 11 ++++++++++- source/util/misc.h | 2 -- 3 files changed, 51 insertions(+), 10 deletions(-) diff --git a/source/loop.cpp b/source/loop.cpp index 8c8ab67..412077f 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -79,6 +79,11 @@ void Loop::Setup(void) // Calculate needed He abundance corrections CalculateAbundanceCorrection(parameters.helium_to_hydrogen_ratio); + if use_variable_abundances + { + ReadRadiativeLossData(); // Initialize the radiative loss arrays + } + //Reserve memory for results results.time.resize(parameters.N); results.heat.resize(parameters.N); @@ -370,22 +375,26 @@ double Loop::CalculateRadiativeLoss(double temperature) double Loop::CalculateRadiativeLoss(double temperature, double abundance_factor) { + double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for radiative loss function + int array_length = sizeof(abundance_array) / sizeof(abundance_array[0]); + + // Find the nearest value of AF to the values in our discrete list, then return that value + int abundance_index = find_closest(abundance_factor, abundance_array, array_length); + + double log_temperature = std::log10(temperature); + array_length = sizeof(log10_temperature_array) / sizeof(log10_temperature_array[0]); + int temperature_index = find_closest(log_temperature, log10_temperature_array, array_length); + return std::pow( 10.0, log10_loss_rate_array[temperature_index, abundance_index] ); } double Loop::CalculateAbundanceFactor(double density, double initial_density); { double initial_abundance_factor = 4.0; // Assumes initially coronal plasma - double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for radiative loss function - int array_length = sizeof(abundance_array) / sizeof(abundance_array[0]); // Calculate using a weighted average of the density // AF = 1.0 + (AF_0 - 1) * (n_0 / n) - double abundance_factor = 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); - - // Find the nearest value of AF to the values in our discrete list, then return that value - int index = find_closest(abundance_factor, abundance_array, array_length); - return abundance_array[index]; + return 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); } double Loop::CalculateCollisionFrequency(double temperature_e,double density) @@ -468,6 +477,31 @@ void Loop::CalculateAbundanceCorrection(double helium_to_hydrogen_ratio) parameters.ion_mass_correction = (1.0 + 4.0*helium_to_hydrogen_ratio)/(2.0 + 3.0*helium_to_hydrogen_ratio)*(1.0 + z_avg)/z_avg; } +void Loop::ReadRadiativeLossData() +{ + std::ifstream fin("data/radiation/abund_10_rad_loss.dat"); + + double number; + char comma; + + for( int i=0; i < 100; ++i) + { + for (int j=0; j < 8; ++j) + { + fin >> number; + if( j == 0 ) + { + log10_temperature_array[i] = number; + } + else + { + log10_loss_rate_array[i,j] = number; + } + fin >> comma; + } + } +} + double Loop::CalculateVelocity(double temperature_e, double temperature_i, double pressure_e) { double c4 = CalculateC4(); diff --git a/source/loop.h b/source/loop.h index aaced7f..8d1a2e5 100644 --- a/source/loop.h +++ b/source/loop.h @@ -29,6 +29,10 @@ class Loop { /* Current state of the system */ state_type __state; + /* Arrays variable abundance radiative losses */ + double log10_loss_rate_array[101,7]; + double log10_temperature_array[101]; + // Calculate c4 // @return ratio of average to base velocity // @@ -52,6 +56,11 @@ class Loop { // void CalculateAbundanceCorrection(double helium_to_hydrogen_ratio); + // Read the csv data file radiative loss rate as a function of temperature and + // abundance factor. + void ReadRadiativeLossData(); + + public: /* Instance of the object */ @@ -242,7 +251,7 @@ class Loop { // also assumed to be initially coronal (AF = 4.0) // // @return the abundance factor (unitless) - static double CalculateAbundanceFactor(double density, double initial_density); + double CalculateAbundanceFactor(double density, double initial_density); }; // Pointer to the class diff --git a/source/util/misc.h b/source/util/misc.h index 261b622..9766e00 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -8,5 +8,3 @@ // Find the closest value in an array int find_closest(double x, double array[], int array_length); -// Read csv data - From 3b089f1b90fcefa3d906399285808122a5a41c2e Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Mon, 18 Dec 2023 16:06:35 -1000 Subject: [PATCH 06/24] lots of debugging to get it to run --- config/ebtel.example.cfg.xml | 1 + data/radiation/abund_10_rad_loss.dat | 1 - source/helper.h | 5 ++++ source/loop.cpp | 38 ++++++++++++++++------------ source/loop.h | 8 ++---- source/util/misc.cpp | 4 +-- source/util/misc.h | 1 + 7 files changed, 33 insertions(+), 25 deletions(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 84fe89e..28c2e93 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -12,6 +12,7 @@ False False False + True ebtel++_results_file.txt 1e-6 0.5 diff --git a/data/radiation/abund_10_rad_loss.dat b/data/radiation/abund_10_rad_loss.dat index bc07666..e6c1cf9 100644 --- a/data/radiation/abund_10_rad_loss.dat +++ b/data/radiation/abund_10_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.309,-23.274,-23.265,-23.259,-23.259,-23.260,-23.260 04.05,-22.814,-22.770,-22.760,-22.751,-22.751,-22.751,-22.751 04.10,-22.365,-22.318,-22.307,-22.299,-22.298,-22.298,-22.298 diff --git a/source/helper.h b/source/helper.h index c04166a..c856a8b 100644 --- a/source/helper.h +++ b/source/helper.h @@ -68,6 +68,11 @@ struct Parameters { double surface_gravity; /* Number of grid points */ size_t N; + + /* Arrays for variable abundance radiative losses */ + double log10_loss_rate_array[101][7]; + double log10_temperature_array[101]; + double initial_density; }; // Structure to hold all results diff --git a/source/loop.cpp b/source/loop.cpp index 412077f..fd08c5e 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -42,6 +42,7 @@ Loop::Loop(char *config) parameters.use_flux_limiting = string2bool(get_element_text(root,"use_flux_limiting")); parameters.calculate_dem = string2bool(get_element_text(root,"calculate_dem")); parameters.use_adaptive_solver = string2bool(get_element_text(root,"use_adaptive_solver")); + parameters.use_variable_abundances = string2bool(get_element_text(root,"use_variable_abundances")); parameters.save_terms = string2bool(get_element_text(root,"save_terms")); //String parameters parameters.output_filename = get_element_text(root,"output_filename"); @@ -79,7 +80,7 @@ void Loop::Setup(void) // Calculate needed He abundance corrections CalculateAbundanceCorrection(parameters.helium_to_hydrogen_ratio); - if use_variable_abundances + if (parameters.use_variable_abundances) { ReadRadiativeLossData(); // Initialize the radiative loss arrays } @@ -140,6 +141,11 @@ state_type Loop::CalculateInitialConditions(void) temperature_old = temperature; density_old = density; } + + if( parameters.use_variable_abundances) + { + parameters.initial_density = density; + } // Set current state in order pressure_e, pressure_i, density pi_initial = parameters.boltzmann_correction*BOLTZMANN_CONSTANT*density*temperature; @@ -198,9 +204,9 @@ void Loop::CalculateDerivs(const state_type &state, state_type &derivs, double t double f_e = CalculateThermalConduction(state[3],state[2],"electron"); double f_i = CalculateThermalConduction(state[4],state[2],"ion"); double radiative_loss; - if use_variable_abundances + if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(state[2], results.density[0]); + double abundance_factor = CalculateAbundanceFactor(state[2]); radiative_loss = CalculateRadiativeLoss(state[3], abundance_factor); } else @@ -278,9 +284,9 @@ void Loop::SaveTerms(void) double f_i = CalculateThermalConduction(__state[4], __state[2], "ion"); double c1 = CalculateC1(__state[3], __state[4], __state[2]); double radiative_loss; - if use_variable_abundances + if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(__state[2], results.density[0]); + double abundance_factor = CalculateAbundanceFactor(__state[2]); radiative_loss = CalculateRadiativeLoss(__state[3], abundance_factor); } else @@ -382,19 +388,19 @@ double Loop::CalculateRadiativeLoss(double temperature, double abundance_factor) int abundance_index = find_closest(abundance_factor, abundance_array, array_length); double log_temperature = std::log10(temperature); - array_length = sizeof(log10_temperature_array) / sizeof(log10_temperature_array[0]); - int temperature_index = find_closest(log_temperature, log10_temperature_array, array_length); + array_length = sizeof(parameters.log10_temperature_array) / sizeof(parameters.log10_temperature_array[0]); + int temperature_index = find_closest(log_temperature, parameters.log10_temperature_array, array_length); - return std::pow( 10.0, log10_loss_rate_array[temperature_index, abundance_index] ); + return std::pow( 10.0, parameters.log10_loss_rate_array[temperature_index][abundance_index] ); } -double Loop::CalculateAbundanceFactor(double density, double initial_density); +double Loop::CalculateAbundanceFactor(double density) { double initial_abundance_factor = 4.0; // Assumes initially coronal plasma // Calculate using a weighted average of the density // AF = 1.0 + (AF_0 - 1) * (n_0 / n) - return 1.0 + (initial_abundance_factor - 1.0) * (initial_density / density); + return 1.0 + (initial_abundance_factor - 1.0) * (parameters.initial_density / density); } double Loop::CalculateCollisionFrequency(double temperature_e,double density) @@ -415,9 +421,9 @@ double Loop::CalculateC1(double temperature_e, double temperature_i, double dens double loss_correction = 1.0; double scale_height = CalculateScaleHeight(temperature_e,temperature_i); double radiative_loss; - if use_variable_abundances + if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(density, results.density[0]); + double abundance_factor = CalculateAbundanceFactor(density); radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); } else @@ -491,11 +497,11 @@ void Loop::ReadRadiativeLossData() fin >> number; if( j == 0 ) { - log10_temperature_array[i] = number; + parameters.log10_temperature_array[i] = number; } else { - log10_loss_rate_array[i,j] = number; + parameters.log10_loss_rate_array[i][j] = number; } fin >> comma; } @@ -508,9 +514,9 @@ double Loop::CalculateVelocity(double temperature_e, double temperature_i, doubl double density = pressure_e/(BOLTZMANN_CONSTANT*temperature_e); double c1 = CalculateC1(temperature_e,temperature_i,density); double radiative_loss; - if use_variable_abundances + if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(density, results.density[0]); + double abundance_factor = CalculateAbundanceFactor(density); radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); } else diff --git a/source/loop.h b/source/loop.h index 8d1a2e5..c62db2c 100644 --- a/source/loop.h +++ b/source/loop.h @@ -29,10 +29,6 @@ class Loop { /* Current state of the system */ state_type __state; - /* Arrays variable abundance radiative losses */ - double log10_loss_rate_array[101,7]; - double log10_temperature_array[101]; - // Calculate c4 // @return ratio of average to base velocity // @@ -71,7 +67,7 @@ class Loop { /* Terms structure */ static Terms terms; - + // Constructor // @config main configuration file // @@ -251,7 +247,7 @@ class Loop { // also assumed to be initially coronal (AF = 4.0) // // @return the abundance factor (unitless) - double CalculateAbundanceFactor(double density, double initial_density); + static double CalculateAbundanceFactor(double density); }; // Pointer to the class diff --git a/source/util/misc.cpp b/source/util/misc.cpp index 97b78f6..7f33df4 100644 --- a/source/util/misc.cpp +++ b/source/util/misc.cpp @@ -9,14 +9,14 @@ Miscellaneous functions for manipulating data int find_closest(double x, double array[], int array_length) { /* Traverses the whole array, so O(N) search. - * A binary search would be more efficient for long arrays O(ln N). + * A binary search would be more efficient for long arrays: O(ln N). */ int index = 0; double closest_value = array[0]; for( int i=1; i < array_length; ++i ) { - if( abs(x - closest_value) > abs(x - array[i]) ) + if( std::abs(x - closest_value) > std::abs(x - array[i]) ) { closest_value = array[i]; index = i; diff --git a/source/util/misc.h b/source/util/misc.h index 9766e00..fcc4290 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -4,6 +4,7 @@ // * // **** +#include // Find the closest value in an array int find_closest(double x, double array[], int array_length); From 9efa88a9c48746b38e89ad6c97f8b972c04d2e40 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 19 Dec 2023 13:29:52 -1000 Subject: [PATCH 07/24] fixed format and input of radiative loss files --- data/radiation/abund_15_rad_loss.dat | 1 - data/radiation/abund_20_rad_loss.dat | 1 - data/radiation/abund_25_rad_loss.dat | 1 - data/radiation/abund_30_rad_loss.dat | 1 - data/radiation/abund_35_rad_loss.dat | 1 - data/radiation/abund_40_rad_loss.dat | 1 - source/helper.h | 2 +- source/loop.cpp | 64 +++++++++++++++++----------- source/loop.h | 2 +- source/util/misc.h | 4 +- 10 files changed, 45 insertions(+), 33 deletions(-) diff --git a/data/radiation/abund_15_rad_loss.dat b/data/radiation/abund_15_rad_loss.dat index a3e9125..8dc5265 100644 --- a/data/radiation/abund_15_rad_loss.dat +++ b/data/radiation/abund_15_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.249,-23.220,-23.213,-23.208,-23.208,-23.209,-23.210 04.05,-22.775,-22.735,-22.726,-22.718,-22.718,-22.718,-22.719 04.10,-22.340,-22.296,-22.286,-22.278,-22.277,-22.277,-22.278 diff --git a/data/radiation/abund_20_rad_loss.dat b/data/radiation/abund_20_rad_loss.dat index c0ab712..70a8bb7 100644 --- a/data/radiation/abund_20_rad_loss.dat +++ b/data/radiation/abund_20_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.199,-23.175,-23.169,-23.164,-23.165,-23.166,-23.167 04.05,-22.741,-22.705,-22.696,-22.690,-22.689,-22.690,-22.690 04.10,-22.318,-22.276,-22.267,-22.259,-22.259,-22.259,-22.259 diff --git a/data/radiation/abund_25_rad_loss.dat b/data/radiation/abund_25_rad_loss.dat index 6feeaec..fec8450 100644 --- a/data/radiation/abund_25_rad_loss.dat +++ b/data/radiation/abund_25_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.152,-23.131,-23.126,-23.122,-23.123,-23.125,-23.125 04.05,-22.708,-22.674,-22.667,-22.661,-22.660,-22.661,-22.662 04.10,-22.296,-22.256,-22.247,-22.240,-22.240,-22.240,-22.241 diff --git a/data/radiation/abund_30_rad_loss.dat b/data/radiation/abund_30_rad_loss.dat index 71ada48..cd5f663 100644 --- a/data/radiation/abund_30_rad_loss.dat +++ b/data/radiation/abund_30_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.110,-23.092,-23.087,-23.084,-23.085,-23.087,-23.088 04.05,-22.677,-22.647,-22.639,-22.634,-22.634,-22.635,-22.636 04.10,-22.275,-22.237,-22.229,-22.222,-22.222,-22.223,-22.223 diff --git a/data/radiation/abund_35_rad_loss.dat b/data/radiation/abund_35_rad_loss.dat index 9172db8..6654221 100644 --- a/data/radiation/abund_35_rad_loss.dat +++ b/data/radiation/abund_35_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.076,-23.060,-23.056,-23.054,-23.054,-23.057,-23.058 04.05,-22.652,-22.623,-22.617,-22.612,-22.612,-22.613,-22.614 04.10,-22.257,-22.221,-22.213,-22.207,-22.207,-22.207,-22.208 diff --git a/data/radiation/abund_40_rad_loss.dat b/data/radiation/abund_40_rad_loss.dat index 1b0a3a5..8a1f015 100644 --- a/data/radiation/abund_40_rad_loss.dat +++ b/data/radiation/abund_40_rad_loss.dat @@ -1,4 +1,3 @@ -Log_t, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 04.00,-23.040,-23.026,-23.023,-23.021,-23.021,-23.024,-23.025 04.05,-22.625,-22.598,-22.592,-22.588,-22.588,-22.589,-22.590 04.10,-22.238,-22.203,-22.196,-22.190,-22.190,-22.191,-22.191 diff --git a/source/helper.h b/source/helper.h index c856a8b..7cb05c7 100644 --- a/source/helper.h +++ b/source/helper.h @@ -70,7 +70,7 @@ struct Parameters { size_t N; /* Arrays for variable abundance radiative losses */ - double log10_loss_rate_array[101][7]; + double log10_loss_rate_array[101][7][7]; double log10_temperature_array[101]; double initial_density; }; diff --git a/source/loop.cpp b/source/loop.cpp index fd08c5e..aa02557 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -206,8 +206,7 @@ void Loop::CalculateDerivs(const state_type &state, state_type &derivs, double t double radiative_loss; if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(state[2]); - radiative_loss = CalculateRadiativeLoss(state[3], abundance_factor); + radiative_loss = CalculateRadiativeLoss(state[3], state[2]); } else { @@ -286,8 +285,7 @@ void Loop::SaveTerms(void) double radiative_loss; if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(__state[2]); - radiative_loss = CalculateRadiativeLoss(__state[3], abundance_factor); + radiative_loss = CalculateRadiativeLoss(__state[3], __state[2]); } else { @@ -379,19 +377,25 @@ double Loop::CalculateRadiativeLoss(double temperature) return chi * std::pow( 10.0, (alpha*log_temperature) ); } -double Loop::CalculateRadiativeLoss(double temperature, double abundance_factor) +double Loop::CalculateRadiativeLoss(double temperature, double density) { - double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for radiative loss function + double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for abundance factors int array_length = sizeof(abundance_array) / sizeof(abundance_array[0]); + double abundance_factor = CalculateAbundanceFactor(density); // Find the nearest value of AF to the values in our discrete list, then return that value int abundance_index = find_closest(abundance_factor, abundance_array, array_length); + + double log10_density_array[] = {8.0, 8.69897, 9.0, 9.69897, 10.0, 10.69897, 11.0}; // Discrete values available for density + double log10_density = std::log10(density); + array_length = sizeof(log10_density_array) / sizeof(log10_density_array[0]); + int density_index = find_closest(log10_density, log10_density_array, array_length); - double log_temperature = std::log10(temperature); + double log10_temperature = std::log10(temperature); array_length = sizeof(parameters.log10_temperature_array) / sizeof(parameters.log10_temperature_array[0]); - int temperature_index = find_closest(log_temperature, parameters.log10_temperature_array, array_length); + int temperature_index = find_closest(log10_temperature, parameters.log10_temperature_array, array_length); - return std::pow( 10.0, parameters.log10_loss_rate_array[temperature_index][abundance_index] ); + return std::pow( 10.0, parameters.log10_loss_rate_array[temperature_index][density_index][abundance_index] ); } double Loop::CalculateAbundanceFactor(double density) @@ -423,8 +427,7 @@ double Loop::CalculateC1(double temperature_e, double temperature_i, double dens double radiative_loss; if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(density); - radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); + radiative_loss = CalculateRadiativeLoss(temperature_e, density); } else { @@ -485,26 +488,40 @@ void Loop::CalculateAbundanceCorrection(double helium_to_hydrogen_ratio) void Loop::ReadRadiativeLossData() { - std::ifstream fin("data/radiation/abund_10_rad_loss.dat"); + // Reads in the radiative loss files. Only need to do once during the setup. + std::string path = "data/radiation/"; + std::string filenames[] = {"abund_10_rad_loss.dat","abund_15_rad_loss.dat","abund_20_rad_loss.dat","abund_25_rad_loss.dat", + "abund_30_rad_loss.dat","abund_35_rad_loss.dat","abund_40_rad_loss.dat"}; + std::ifstream fin; + int n_files = sizeof(filenames)/sizeof(filenames[0]); double number; char comma; - for( int i=0; i < 100; ++i) + for (int i=0; i < n_files; ++i) // Loop over files for different abundances { - for (int j=0; j < 8; ++j) + fin.open(path+filenames[i]); + + for (int j=0; j < 100; ++j) // Loop over temperatures (rows in files) { - fin >> number; - if( j == 0 ) - { - parameters.log10_temperature_array[i] = number; - } - else + + for (int k=0; k < 8; ++k) // Loop over densities (columns in files) { - parameters.log10_loss_rate_array[i][j] = number; + fin >> number; + if ( k == 0 ) + { + parameters.log10_temperature_array[j] = number; + } + else + { + parameters.log10_loss_rate_array[j][k][i] = number; + //[temperature_index][density_index][abundance_index] + } + fin >> comma; } - fin >> comma; } + fin.close(); + fin.clear(); } } @@ -516,8 +533,7 @@ double Loop::CalculateVelocity(double temperature_e, double temperature_i, doubl double radiative_loss; if (parameters.use_variable_abundances) { - double abundance_factor = CalculateAbundanceFactor(density); - radiative_loss = CalculateRadiativeLoss(temperature_e, abundance_factor); + radiative_loss = CalculateRadiativeLoss(temperature_e, density); } else { diff --git a/source/loop.h b/source/loop.h index c62db2c..3d59edc 100644 --- a/source/loop.h +++ b/source/loop.h @@ -215,7 +215,7 @@ class Loop { // @return radiative loss function (in erg cm$^3$ s$^{-1}$) // static double CalculateRadiativeLoss(double temperature); - static double CalculateRadiativeLoss(double temperature, double abundance_factor); + static double CalculateRadiativeLoss(double temperature, double density); // Calculate derivatives of EBTEL equations // @state current state of the loop diff --git a/source/util/misc.h b/source/util/misc.h index fcc4290..db813b0 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -6,6 +6,8 @@ #include -// Find the closest value in an array +// Find the closest value in an array to the input x +// Uses a simple search, O(N), so it's best for short arrays +// Returns the index of the array that's closest to value x int find_closest(double x, double array[], int array_length); From f3ca13e6dd8868ac5f2d63be05fe3d416022697a Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Wed, 10 Jan 2024 14:45:10 -1000 Subject: [PATCH 08/24] add correction for draining and later upflow events --- source/helper.h | 4 ++++ source/loop.cpp | 38 +++++++++++++++++++++++++++++++++++--- 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/source/helper.h b/source/helper.h index 7cb05c7..30e4e81 100644 --- a/source/helper.h +++ b/source/helper.h @@ -73,6 +73,10 @@ struct Parameters { double log10_loss_rate_array[101][7][7]; double log10_temperature_array[101]; double initial_density; + double previous_density; + double initial_abundance_factor; + double previous_abundance_factor; + bool upflowing; }; // Structure to hold all results diff --git a/source/loop.cpp b/source/loop.cpp index aa02557..654e281 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -145,6 +145,10 @@ state_type Loop::CalculateInitialConditions(void) if( parameters.use_variable_abundances) { parameters.initial_density = density; + parameters.previous_density = density; + parameters.initial_abundance_factor = 4.0; // Assumes initially coronal plasma + parameters.previous_abundance_factor = parameters.initial_abundance_factor; + parameters.upflowing = false; } // Set current state in order pressure_e, pressure_i, density @@ -274,6 +278,11 @@ void Loop::SaveResults(int i,double time) results.density[i] = __state[2]; results.velocity[i] = velocity; } + + if( parameters.use_variable_abundances ) + { + parameters.previous_density = __state[2]; + } } void Loop::SaveTerms(void) @@ -400,11 +409,34 @@ double Loop::CalculateRadiativeLoss(double temperature, double density) double Loop::CalculateAbundanceFactor(double density) { - double initial_abundance_factor = 4.0; // Assumes initially coronal plasma - // Calculate using a weighted average of the density // AF = 1.0 + (AF_0 - 1) * (n_0 / n) - return 1.0 + (initial_abundance_factor - 1.0) * (parameters.initial_density / density); + double abundance_factor; + + if( density > parameters.previous_density && !parameters.upflowing ) + { + // When the plasma starts to upflow, store the coronal density as this + // is the "initial" density needed for the calculation of AF + parameters.upflowing = true; + parameters.initial_density = parameters.previous_density; + parameters.initial_abundance_factor = parameters.previous_abundance_factor; + } + + if( density <= parameters.previous_density ) + { + // If the density has not increased, it is no longer upflowing, and so the + // AF does not change since elements will not preferentially drain + parameters.upflowing = false; + abundance_factor = parameters.previous_abundance_factor; + } + + if( parameters.upflowing ) + { + abundance_factor = 1.0 + (parameters.initial_abundance_factor - 1.0) * (parameters.initial_density / density); + } + + parameters.previous_abundance_factor = abundance_factor; + return abundance_factor; } double Loop::CalculateCollisionFrequency(double temperature_e,double density) From a10895dc330efd4ea3ba7aa09111cfd5ec8315a0 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Thu, 11 Jan 2024 13:03:08 -1000 Subject: [PATCH 09/24] use original radiation for initial conditions --- source/helper.h | 2 +- source/loop.cpp | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/source/helper.h b/source/helper.h index 30e4e81..7289a1f 100644 --- a/source/helper.h +++ b/source/helper.h @@ -76,7 +76,7 @@ struct Parameters { double previous_density; double initial_abundance_factor; double previous_abundance_factor; - bool upflowing; + bool upflowing, initial_radiation; }; // Structure to hold all results diff --git a/source/loop.cpp b/source/loop.cpp index 654e281..3bb4966 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -122,6 +122,11 @@ state_type Loop::CalculateInitialConditions(void) double heat = heater->Get_Heating(0.0); state_type state; + if( parameters.use_variable_abundances ) + { + parameters.initial_radiation = true; + } + while(i 0) @@ -149,6 +154,7 @@ state_type Loop::CalculateInitialConditions(void) parameters.initial_abundance_factor = 4.0; // Assumes initially coronal plasma parameters.previous_abundance_factor = parameters.initial_abundance_factor; parameters.upflowing = false; + parameters.initial_radiation = false; } // Set current state in order pressure_e, pressure_i, density @@ -457,7 +463,7 @@ double Loop::CalculateC1(double temperature_e, double temperature_i, double dens double loss_correction = 1.0; double scale_height = CalculateScaleHeight(temperature_e,temperature_i); double radiative_loss; - if (parameters.use_variable_abundances) + if (parameters.use_variable_abundances && !parameters.initial_radiation) { radiative_loss = CalculateRadiativeLoss(temperature_e, density); } From 5ebfa610521f6278f09fb3575654abc5256e51bf Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Thu, 11 Jan 2024 16:15:22 -1000 Subject: [PATCH 10/24] add include guard for misc.h --- source/util/misc.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/source/util/misc.h b/source/util/misc.h index db813b0..4fac1cc 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -4,6 +4,9 @@ // * // **** +#ifndef MISC_H +#define MISC_H + #include // Find the closest value in an array to the input x @@ -11,3 +14,4 @@ // Returns the index of the array that's closest to value x int find_closest(double x, double array[], int array_length); +#endif \ No newline at end of file From 37127408e53835b6e060875d23e4e1006db78430 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Sat, 11 May 2024 15:30:12 -0400 Subject: [PATCH 11/24] change bool use_variable_abundances to string radiation input --- config/ebtel.example.cfg.xml | 2 +- source/helper.h | 7 +++++++ source/loop.cpp | 24 +++++++++++++++++++++++- 3 files changed, 31 insertions(+), 2 deletions(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 28c2e93..254e071 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -12,7 +12,7 @@ False False False - True + variable ebtel++_results_file.txt 1e-6 0.5 diff --git a/source/helper.h b/source/helper.h index 7289a1f..afa257c 100644 --- a/source/helper.h +++ b/source/helper.h @@ -52,6 +52,13 @@ struct Parameters { bool save_terms; /* Switch for using the adaptive solver option */ bool use_adaptive_solver; + /* What radiative losses to assume: + power_law: use the default power-law fit to radiative losses (Klimchuk et al 2008) + variable: use a look-up table with time-variable abundance factor f + photospheric: use a look-up table with abundance factor f = 1 + coronal: use a look-up table with abundance factor f = 4 + */ + std::string radiation; /* Switch for using time-variable abundances */ bool use_variable_abundances; /* Path to output file */ diff --git a/source/loop.cpp b/source/loop.cpp index 3bb4966..7acf36f 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -42,7 +42,24 @@ Loop::Loop(char *config) parameters.use_flux_limiting = string2bool(get_element_text(root,"use_flux_limiting")); parameters.calculate_dem = string2bool(get_element_text(root,"calculate_dem")); parameters.use_adaptive_solver = string2bool(get_element_text(root,"use_adaptive_solver")); - parameters.use_variable_abundances = string2bool(get_element_text(root,"use_variable_abundances")); + parameters.radiation = get_element_text(root,"radiation"); + if (parameters.radiation == "power_law") + { + parameters.use_variable_abundances = false; + } + else if( (parameters.radiation == "variable") || + (parameters.radiation == "photospheric") || + (parameters.radiation == "coronal") ) + { + parameters.use_variable_abundances = true; + } + else + { + std::string filename(config); + std::string error_message = "Invalid option for radiation in "+filename+ + ".\n Valid options are power_law, variable, photospheric, or coronal."; + throw std::runtime_error(error_message); + } parameters.save_terms = string2bool(get_element_text(root,"save_terms")); //String parameters parameters.output_filename = get_element_text(root,"output_filename"); @@ -415,6 +432,11 @@ double Loop::CalculateRadiativeLoss(double temperature, double density) double Loop::CalculateAbundanceFactor(double density) { + if (parameters.radiation == "photospheric") + return 1.0; + if (parameters.radiation == "coronal") + return 4.0; + // Calculate using a weighted average of the density // AF = 1.0 + (AF_0 - 1) * (n_0 / n) double abundance_factor; From f2117d0af41d2eae798fd5f8edc07ff7548d8543 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Sun, 2 Jun 2024 22:16:03 -1000 Subject: [PATCH 12/24] add radiation to test configurations --- SConstruct | 8 ++++++-- config/ebtel.example.cfg.xml | 10 ++++------ tests/test_compare_hydrad.py | 1 + tests/test_compare_idl.py | 1 + tests/test_single_fluid.py | 1 + tests/test_solver.py | 1 + 6 files changed, 14 insertions(+), 8 deletions(-) diff --git a/SConstruct b/SConstruct index 72cded6..2d9a0b4 100644 --- a/SConstruct +++ b/SConstruct @@ -31,9 +31,13 @@ env = Environment(CXX=CXX, CXXFLAGS=cxx_flags) if 'darwin' in sys.platform: print("Using Mac OS X compile options.") - env.Append(CPPPATH=['/opt/local/include', '/usr/include/malloc']) + if 'HOMEBREW_PREFIX' in os.environ: + env.Append(CPPPATH=['/opt/local/include','/usr/include/malloc','/opt/homebrew/include']) + env.Append(LIBPATH=['/opt/homebrew/lib']) + else: + env.Append(CPPPATH=['/opt/local/include', '/usr/include/malloc']) + env.Append(LIBPATH=['/opt/local/lib']) env.Append(LIBS=['boost_program_options-mt']) - env.Append(LIBPATH=['/opt/local/lib']) elif 'linux' in sys.platform: print("Using Linux compile options.") env.Append(CPPPATH=['/usr/include']) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 254e071..52d3112 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -3,7 +3,7 @@ 5000.0 1.0 1e+300 - 40.0e+8 + 80e+8 1.0 False True @@ -12,8 +12,8 @@ False False False - variable - ebtel++_results_file.txt + coronal + /Users/reep/Documents/Forks/ebtel_abundances/src/data/cor_L80_H0.01_t20.txt 1e-6 0.5 2.0 @@ -28,9 +28,7 @@ 3.5e-5 1.0 - - - + diff --git a/tests/test_compare_hydrad.py b/tests/test_compare_hydrad.py index 4e9b1fc..7bf3ac2 100644 --- a/tests/test_compare_hydrad.py +++ b/tests/test_compare_hydrad.py @@ -26,6 +26,7 @@ def base_config(): 'calculate_dem': False, 'save_terms': False, 'use_adaptive_solver': True, + 'radiation': 'power_law', 'adaptive_solver_error': 1e-8, 'adaptive_solver_safety': 0.5, 'c1_cond0': 6.0, diff --git a/tests/test_compare_idl.py b/tests/test_compare_idl.py index 1dea552..5938a59 100644 --- a/tests/test_compare_idl.py +++ b/tests/test_compare_idl.py @@ -23,6 +23,7 @@ def base_config(): 'use_flux_limiting': False, 'use_adaptive_solver': False, 'calculate_dem': False, + 'radiation': 'power_law', 'save_terms': False, 'adaptive_solver_error': 1e-6, 'adaptive_solver_safety': 0.5, diff --git a/tests/test_single_fluid.py b/tests/test_single_fluid.py index 2938963..198615c 100644 --- a/tests/test_single_fluid.py +++ b/tests/test_single_fluid.py @@ -25,6 +25,7 @@ def base_config(): 'calculate_dem': False, 'save_terms': False, 'use_adaptive_solver': True, + 'radiation': 'power_law', 'adaptive_solver_error': 1e-6, 'adaptive_solver_safety': 0.5, 'c1_cond0': 2.0, diff --git a/tests/test_solver.py b/tests/test_solver.py index 8d82620..e6d833d 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -23,6 +23,7 @@ def base_config(): 'use_flux_limiting': True, 'calculate_dem': True, 'save_terms': False, + 'radiation': 'power_law', 'adaptive_solver_error': 1e-6, 'adaptive_solver_safety': 0.5, 'c1_cond0': 2.0, From 2878be1d3a0129434751aca4fe778b64cc92d107 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Sun, 2 Jun 2024 22:35:30 -1000 Subject: [PATCH 13/24] revert config file to original values --- config/ebtel.example.cfg.xml | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 52d3112..3ec97cf 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -3,7 +3,7 @@ 5000.0 1.0 1e+300 - 80e+8 + 40.0e+8 1.0 False True @@ -13,7 +13,7 @@ False False coronal - /Users/reep/Documents/Forks/ebtel_abundances/src/data/cor_L80_H0.01_t20.txt + ebtel++_results_file.txt 1e-6 0.5 2.0 @@ -28,7 +28,9 @@ 3.5e-5 1.0 - + + + From 076e95eb621b612631f350667910a72c526668e4 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Sun, 2 Jun 2024 22:36:42 -1000 Subject: [PATCH 14/24] use the power law by default in config --- config/ebtel.example.cfg.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 3ec97cf..b60cd8b 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -12,7 +12,7 @@ False False False - coronal + power_law ebtel++_results_file.txt 1e-6 0.5 From c9fe089a391c34544e6af6670d6c87ac89326b98 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 14:53:57 -1000 Subject: [PATCH 15/24] add descriptions of variables in helper.h --- config/ebtel.example.cfg.xml | 10 ++++------ source/helper.h | 21 ++++++++++++++++++--- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index b60cd8b..52d3112 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -3,7 +3,7 @@ 5000.0 1.0 1e+300 - 40.0e+8 + 80e+8 1.0 False True @@ -12,8 +12,8 @@ False False False - power_law - ebtel++_results_file.txt + coronal + /Users/reep/Documents/Forks/ebtel_abundances/src/data/cor_L80_H0.01_t20.txt 1e-6 0.5 2.0 @@ -28,9 +28,7 @@ 3.5e-5 1.0 - - - + diff --git a/source/helper.h b/source/helper.h index afa257c..a4a9ca7 100644 --- a/source/helper.h +++ b/source/helper.h @@ -76,14 +76,29 @@ struct Parameters { /* Number of grid points */ size_t N; - /* Arrays for variable abundance radiative losses */ - double log10_loss_rate_array[101][7][7]; + /* Variables and arrays used for variable abundance radiative losses */ + /* The temperature values in the look-up table for radiative losses */ double log10_temperature_array[101]; + /* The look-up table's radiative loss rate as a + * function of [temperature][density][abundance] */ + double log10_loss_rate_array[101][7][7]; + /* The density in the corona before upflows begin, used to calculate + * the change in abundance factor */ double initial_density; + /* The density at the previous time step*/ double previous_density; + /* The abundance factor in the corona before upflows begin */ double initial_abundance_factor; + /* The abundance factor at the previous time step */ double previous_abundance_factor; - bool upflowing, initial_radiation; + /* Whether the flows are upflowing (into the corona) or not, + * which is used to determine whether the abundance factor changes + * with time. That is, flows out of the corona do not affect + * the abundance factor. */ + bool upflowing; + /* The power law losses are used to calculate the initial conditions, so + * this bool is used to tell the code that. */ + bool initial_radiation; }; // Structure to hold all results From 615e5213adbb417551a0cffb37dd04071832a08c Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 14:59:38 -1000 Subject: [PATCH 16/24] revert the config file --- config/ebtel.example.cfg.xml | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/config/ebtel.example.cfg.xml b/config/ebtel.example.cfg.xml index 52d3112..eb175bb 100644 --- a/config/ebtel.example.cfg.xml +++ b/config/ebtel.example.cfg.xml @@ -3,7 +3,7 @@ 5000.0 1.0 1e+300 - 80e+8 + 40e+8 1.0 False True @@ -12,8 +12,8 @@ False False False - coronal - /Users/reep/Documents/Forks/ebtel_abundances/src/data/cor_L80_H0.01_t20.txt + power_law + ebtel++_results_file.txt 1e-6 0.5 2.0 @@ -28,7 +28,9 @@ 3.5e-5 1.0 - + + + From 1bf1db71896e4b27b498b29d5ab0e9cfda73205d Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 15:26:23 -1000 Subject: [PATCH 17/24] rename function --- source/loop.cpp | 4 ++-- source/loop.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/source/loop.cpp b/source/loop.cpp index 7acf36f..518895a 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -95,7 +95,7 @@ Loop::~Loop(void) void Loop::Setup(void) { // Calculate needed He abundance corrections - CalculateAbundanceCorrection(parameters.helium_to_hydrogen_ratio); + CalculateIonMassCorrection(parameters.helium_to_hydrogen_ratio); if (parameters.use_variable_abundances) { @@ -539,7 +539,7 @@ double Loop::CalculateScaleHeight(double temperature_e,double temperature_i) return BOLTZMANN_CONSTANT*(temperature_e + parameters.boltzmann_correction*temperature_i)/(parameters.ion_mass_correction*PROTON_MASS)/(parameters.surface_gravity * (double)SOLAR_SURFACE_GRAVITY); } -void Loop::CalculateAbundanceCorrection(double helium_to_hydrogen_ratio) +void Loop::CalculateIonMassCorrection(double helium_to_hydrogen_ratio) { double z_avg = (1.0 + 2.0*helium_to_hydrogen_ratio)/(1.0 + helium_to_hydrogen_ratio); parameters.boltzmann_correction = 1.0/z_avg; diff --git a/source/loop.h b/source/loop.h index 3d59edc..410cbca 100644 --- a/source/loop.h +++ b/source/loop.h @@ -48,9 +48,9 @@ class Loop { // static double CalculateCollisionFrequency(double temperature_e,double density); - // Calculate correction for He abundance + // Calculate correction to the ion mass for He abundance // - void CalculateAbundanceCorrection(double helium_to_hydrogen_ratio); + void CalculateIonMassCorrection(double helium_to_hydrogen_ratio); // Read the csv data file radiative loss rate as a function of temperature and // abundance factor. From 9f3d71b7308fb85ce012db1de01851126687eb0f Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 16:43:48 -1000 Subject: [PATCH 18/24] remove hardcoding of radiation filenames and number of files -- requires > C++ 17 --- SConstruct | 2 +- source/helper.h | 5 +++++ source/loop.cpp | 25 ++++++++++++++++++++----- 3 files changed, 26 insertions(+), 6 deletions(-) diff --git a/SConstruct b/SConstruct index 2d9a0b4..17f2c11 100644 --- a/SConstruct +++ b/SConstruct @@ -18,7 +18,7 @@ AddOption('--linkflags', dest='linkflags', type='string', nargs=1, action='store subdirs = ['source'] -cxx_flags = ['-std=c++14',] +cxx_flags = ['-std=c++17',] try: CXX = os.environ['CXX'] except KeyError: diff --git a/source/helper.h b/source/helper.h index a4a9ca7..9101713 100644 --- a/source/helper.h +++ b/source/helper.h @@ -15,9 +15,14 @@ General purpose includes to be used everywhere #include #include #include +#include +#include #include "boost/array.hpp" #include "util/xmlreader.h" +namespace fs = std::filesystem; + + // Structure to hold all input parameters struct Parameters { /* Total simulation time (in s) */ diff --git a/source/loop.cpp b/source/loop.cpp index 518895a..402e029 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -550,17 +550,32 @@ void Loop::ReadRadiativeLossData() { // Reads in the radiative loss files. Only need to do once during the setup. std::string path = "data/radiation/"; - std::string filenames[] = {"abund_10_rad_loss.dat","abund_15_rad_loss.dat","abund_20_rad_loss.dat","abund_25_rad_loss.dat", - "abund_30_rad_loss.dat","abund_35_rad_loss.dat","abund_40_rad_loss.dat"}; + std::vector filenames; std::ifstream fin; - int n_files = sizeof(filenames)/sizeof(filenames[0]); + std::string filename; + /* We use a set here to read in the filenames because each filename is unique, + * and this will therefore sort automatically. */ + std::set file_set; + + // Read in the filenames of each file in the radiation directory: + for (const auto & entry : fs::directory_iterator(path)) + { + file_set.insert(entry.path()); + } + for (auto &file : file_set) + { + filenames.push_back(file.c_str()); + } + int n_abund = filenames.size(); + double number; char comma; - for (int i=0; i < n_files; ++i) // Loop over files for different abundances + for (int i=0; i < n_abund; ++i) // Loop over files for different abundances { - fin.open(path+filenames[i]); + //fin.open(path+filenames[i]); + fin.open(filenames[i]); for (int j=0; j < 100; ++j) // Loop over temperatures (rows in files) { From 4cbd02f8299fc1f4a9ce7467d84da43ad340eb39 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 16:58:03 -1000 Subject: [PATCH 19/24] use_variable_abundances -> use_lookup_table_losses --- source/helper.h | 4 ++-- source/loop.cpp | 28 ++++++++++++++-------------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/source/helper.h b/source/helper.h index 9101713..f7254e0 100644 --- a/source/helper.h +++ b/source/helper.h @@ -64,8 +64,8 @@ struct Parameters { coronal: use a look-up table with abundance factor f = 4 */ std::string radiation; - /* Switch for using time-variable abundances */ - bool use_variable_abundances; + /* Switch for using look-up tables for the radiative loss calculation */ + bool use_lookup_table_losses; /* Path to output file */ std::string output_filename; /* XML node holding DEM calculation parameters */ diff --git a/source/loop.cpp b/source/loop.cpp index 402e029..1aaccce 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -45,13 +45,13 @@ Loop::Loop(char *config) parameters.radiation = get_element_text(root,"radiation"); if (parameters.radiation == "power_law") { - parameters.use_variable_abundances = false; + parameters.use_lookup_table_losses = false; } else if( (parameters.radiation == "variable") || (parameters.radiation == "photospheric") || (parameters.radiation == "coronal") ) { - parameters.use_variable_abundances = true; + parameters.use_lookup_table_losses = true; } else { @@ -97,7 +97,7 @@ void Loop::Setup(void) // Calculate needed He abundance corrections CalculateIonMassCorrection(parameters.helium_to_hydrogen_ratio); - if (parameters.use_variable_abundances) + if (parameters.use_lookup_table_losses) { ReadRadiativeLossData(); // Initialize the radiative loss arrays } @@ -139,7 +139,7 @@ state_type Loop::CalculateInitialConditions(void) double heat = heater->Get_Heating(0.0); state_type state; - if( parameters.use_variable_abundances ) + if( parameters.use_lookup_table_losses ) { parameters.initial_radiation = true; } @@ -164,7 +164,7 @@ state_type Loop::CalculateInitialConditions(void) density_old = density; } - if( parameters.use_variable_abundances) + if( parameters.use_lookup_table_losses ) { parameters.initial_density = density; parameters.previous_density = density; @@ -231,7 +231,7 @@ void Loop::CalculateDerivs(const state_type &state, state_type &derivs, double t double f_e = CalculateThermalConduction(state[3],state[2],"electron"); double f_i = CalculateThermalConduction(state[4],state[2],"ion"); double radiative_loss; - if (parameters.use_variable_abundances) + if (parameters.use_lookup_table_losses) { radiative_loss = CalculateRadiativeLoss(state[3], state[2]); } @@ -302,7 +302,7 @@ void Loop::SaveResults(int i,double time) results.velocity[i] = velocity; } - if( parameters.use_variable_abundances ) + if( parameters.use_lookup_table_losses ) { parameters.previous_density = __state[2]; } @@ -315,7 +315,7 @@ void Loop::SaveTerms(void) double f_i = CalculateThermalConduction(__state[4], __state[2], "ion"); double c1 = CalculateC1(__state[3], __state[4], __state[2]); double radiative_loss; - if (parameters.use_variable_abundances) + if (parameters.use_lookup_table_losses) { radiative_loss = CalculateRadiativeLoss(__state[3], __state[2]); } @@ -485,7 +485,7 @@ double Loop::CalculateC1(double temperature_e, double temperature_i, double dens double loss_correction = 1.0; double scale_height = CalculateScaleHeight(temperature_e,temperature_i); double radiative_loss; - if (parameters.use_variable_abundances && !parameters.initial_radiation) + if (parameters.use_lookup_table_losses && !parameters.initial_radiation) { radiative_loss = CalculateRadiativeLoss(temperature_e, density); } @@ -549,7 +549,7 @@ void Loop::CalculateIonMassCorrection(double helium_to_hydrogen_ratio) void Loop::ReadRadiativeLossData() { // Reads in the radiative loss files. Only need to do once during the setup. - std::string path = "data/radiation/"; + std::string data_path = "data/radiation/"; std::vector filenames; std::ifstream fin; std::string filename; @@ -558,8 +558,9 @@ void Loop::ReadRadiativeLossData() * and this will therefore sort automatically. */ std::set file_set; - // Read in the filenames of each file in the radiation directory: - for (const auto & entry : fs::directory_iterator(path)) + /* Read in the filenames of each file in the radiation directory. + * fs::directory_iterator requires C++ 17 or newer. */ + for (const auto & entry : fs::directory_iterator(data_path)) { file_set.insert(entry.path()); } @@ -574,7 +575,6 @@ void Loop::ReadRadiativeLossData() for (int i=0; i < n_abund; ++i) // Loop over files for different abundances { - //fin.open(path+filenames[i]); fin.open(filenames[i]); for (int j=0; j < 100; ++j) // Loop over temperatures (rows in files) @@ -606,7 +606,7 @@ double Loop::CalculateVelocity(double temperature_e, double temperature_i, doubl double density = pressure_e/(BOLTZMANN_CONSTANT*temperature_e); double c1 = CalculateC1(temperature_e,temperature_i,density); double radiative_loss; - if (parameters.use_variable_abundances) + if (parameters.use_lookup_table_losses) { radiative_loss = CalculateRadiativeLoss(temperature_e, density); } From c5a31075088128b01d47534a68b0781c30be2d85 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Tue, 2 Jul 2024 21:24:33 -1000 Subject: [PATCH 20/24] remove hardcoding of indices --- data/radiation/abund_10_rad_loss.dat | 1 + data/radiation/abund_15_rad_loss.dat | 1 + data/radiation/abund_20_rad_loss.dat | 1 + data/radiation/abund_25_rad_loss.dat | 1 + data/radiation/abund_30_rad_loss.dat | 1 + data/radiation/abund_35_rad_loss.dat | 1 + data/radiation/abund_40_rad_loss.dat | 1 + source/helper.h | 8 ++- source/loop.cpp | 80 ++++++++++++++++++++-------- source/util/misc.cpp | 22 +++++++- source/util/misc.h | 6 ++- 11 files changed, 95 insertions(+), 28 deletions(-) diff --git a/data/radiation/abund_10_rad_loss.dat b/data/radiation/abund_10_rad_loss.dat index e6c1cf9..91545c2 100644 --- a/data/radiation/abund_10_rad_loss.dat +++ b/data/radiation/abund_10_rad_loss.dat @@ -1,3 +1,4 @@ +1.0,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.309,-23.274,-23.265,-23.259,-23.259,-23.260,-23.260 04.05,-22.814,-22.770,-22.760,-22.751,-22.751,-22.751,-22.751 04.10,-22.365,-22.318,-22.307,-22.299,-22.298,-22.298,-22.298 diff --git a/data/radiation/abund_15_rad_loss.dat b/data/radiation/abund_15_rad_loss.dat index 8dc5265..287f789 100644 --- a/data/radiation/abund_15_rad_loss.dat +++ b/data/radiation/abund_15_rad_loss.dat @@ -1,3 +1,4 @@ +1.5,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.249,-23.220,-23.213,-23.208,-23.208,-23.209,-23.210 04.05,-22.775,-22.735,-22.726,-22.718,-22.718,-22.718,-22.719 04.10,-22.340,-22.296,-22.286,-22.278,-22.277,-22.277,-22.278 diff --git a/data/radiation/abund_20_rad_loss.dat b/data/radiation/abund_20_rad_loss.dat index 70a8bb7..89173cf 100644 --- a/data/radiation/abund_20_rad_loss.dat +++ b/data/radiation/abund_20_rad_loss.dat @@ -1,3 +1,4 @@ +2.0,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.199,-23.175,-23.169,-23.164,-23.165,-23.166,-23.167 04.05,-22.741,-22.705,-22.696,-22.690,-22.689,-22.690,-22.690 04.10,-22.318,-22.276,-22.267,-22.259,-22.259,-22.259,-22.259 diff --git a/data/radiation/abund_25_rad_loss.dat b/data/radiation/abund_25_rad_loss.dat index fec8450..70b1fb5 100644 --- a/data/radiation/abund_25_rad_loss.dat +++ b/data/radiation/abund_25_rad_loss.dat @@ -1,3 +1,4 @@ +2.5,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.152,-23.131,-23.126,-23.122,-23.123,-23.125,-23.125 04.05,-22.708,-22.674,-22.667,-22.661,-22.660,-22.661,-22.662 04.10,-22.296,-22.256,-22.247,-22.240,-22.240,-22.240,-22.241 diff --git a/data/radiation/abund_30_rad_loss.dat b/data/radiation/abund_30_rad_loss.dat index cd5f663..3d7326e 100644 --- a/data/radiation/abund_30_rad_loss.dat +++ b/data/radiation/abund_30_rad_loss.dat @@ -1,3 +1,4 @@ +3.0,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.110,-23.092,-23.087,-23.084,-23.085,-23.087,-23.088 04.05,-22.677,-22.647,-22.639,-22.634,-22.634,-22.635,-22.636 04.10,-22.275,-22.237,-22.229,-22.222,-22.222,-22.223,-22.223 diff --git a/data/radiation/abund_35_rad_loss.dat b/data/radiation/abund_35_rad_loss.dat index 6654221..b683255 100644 --- a/data/radiation/abund_35_rad_loss.dat +++ b/data/radiation/abund_35_rad_loss.dat @@ -1,3 +1,4 @@ +3.5,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.076,-23.060,-23.056,-23.054,-23.054,-23.057,-23.058 04.05,-22.652,-22.623,-22.617,-22.612,-22.612,-22.613,-22.614 04.10,-22.257,-22.221,-22.213,-22.207,-22.207,-22.207,-22.208 diff --git a/data/radiation/abund_40_rad_loss.dat b/data/radiation/abund_40_rad_loss.dat index 8a1f015..f58ddbc 100644 --- a/data/radiation/abund_40_rad_loss.dat +++ b/data/radiation/abund_40_rad_loss.dat @@ -1,3 +1,4 @@ +4.0,8.0,8.69897,9.0,9.69897,10.0,10.69897,11.0 04.00,-23.040,-23.026,-23.023,-23.021,-23.021,-23.024,-23.025 04.05,-22.625,-22.598,-22.592,-22.588,-22.588,-22.589,-22.590 04.10,-22.238,-22.203,-22.196,-22.190,-22.190,-22.191,-22.191 diff --git a/source/helper.h b/source/helper.h index f7254e0..9738e3d 100644 --- a/source/helper.h +++ b/source/helper.h @@ -83,10 +83,14 @@ struct Parameters { /* Variables and arrays used for variable abundance radiative losses */ /* The temperature values in the look-up table for radiative losses */ - double log10_temperature_array[101]; + //double log10_temperature_array[101]; + std::vector log10_temperature_array; + std::vector log10_density_array; + std::vector abundance_array; /* The look-up table's radiative loss rate as a * function of [temperature][density][abundance] */ - double log10_loss_rate_array[101][7][7]; + std::vector > > log10_loss_rate_array; + //double log10_loss_rate_array[101][7][7]; /* The density in the corona before upflows begin, used to calculate * the change in abundance factor */ double initial_density; diff --git a/source/loop.cpp b/source/loop.cpp index 1aaccce..16109ea 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -411,23 +411,21 @@ double Loop::CalculateRadiativeLoss(double temperature) double Loop::CalculateRadiativeLoss(double temperature, double density) { - double abundance_array[] = {1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0}; // Discrete values available for abundance factors - int array_length = sizeof(abundance_array) / sizeof(abundance_array[0]); + int array_length = parameters.abundance_array.size(); double abundance_factor = CalculateAbundanceFactor(density); // Find the nearest value of AF to the values in our discrete list, then return that value - int abundance_index = find_closest(abundance_factor, abundance_array, array_length); + int abundance_index = find_closest(abundance_factor, parameters.abundance_array, array_length); - double log10_density_array[] = {8.0, 8.69897, 9.0, 9.69897, 10.0, 10.69897, 11.0}; // Discrete values available for density double log10_density = std::log10(density); - array_length = sizeof(log10_density_array) / sizeof(log10_density_array[0]); - int density_index = find_closest(log10_density, log10_density_array, array_length); + array_length = parameters.log10_density_array.size(); + int density_index = find_closest(log10_density, parameters.log10_density_array, array_length); double log10_temperature = std::log10(temperature); - array_length = sizeof(parameters.log10_temperature_array) / sizeof(parameters.log10_temperature_array[0]); + array_length = parameters.log10_temperature_array.size(); int temperature_index = find_closest(log10_temperature, parameters.log10_temperature_array, array_length); - return std::pow( 10.0, parameters.log10_loss_rate_array[temperature_index][density_index][abundance_index] ); + return std::pow( 10.0, parameters.log10_loss_rate_array[abundance_index][temperature_index][density_index] ); } double Loop::CalculateAbundanceFactor(double density) @@ -551,9 +549,14 @@ void Loop::ReadRadiativeLossData() // Reads in the radiative loss files. Only need to do once during the setup. std::string data_path = "data/radiation/"; std::vector filenames; + std::vector loss_rate_1D; + std::vector > loss_rate_2D; std::ifstream fin; - std::string filename; - + std::string filename, line; + int i, j, k; + double number; + char comma; + /* We use a set here to read in the filenames because each filename is unique, * and this will therefore sort automatically. */ std::set file_set; @@ -569,35 +572,68 @@ void Loop::ReadRadiativeLossData() filenames.push_back(file.c_str()); } int n_abund = filenames.size(); - - double number; - char comma; - for (int i=0; i < n_abund; ++i) // Loop over files for different abundances + /* Find the number of rows and columns in the file to determine + * the number of density and temperature points in the look-up tables. */ + int n_temperature = 0; + int n_density = 0; + fin.open(filenames[0]); + getline(fin, line); + for( i=0; i < line.size(); ++i ) + { + if( line[i] == ',' ) n_density++; + } + while( getline(fin, line) ) n_temperature++; + fin.close(); + fin.clear(); + + for (i=0; i < n_abund; ++i) // Loop over files for different abundances { fin.open(filenames[i]); - for (int j=0; j < 100; ++j) // Loop over temperatures (rows in files) + for(k=0; k < n_density+1; ++k) // Read the first row to get the abundance factor + { + fin >> number; + fin >> comma; + if ( k == 0 ) + { + parameters.abundance_array.push_back(number); + } + else if ( i == 0 && k > 0 ) + { + parameters.log10_density_array.push_back(number); + } + } + + for (j=0; j < n_temperature; ++j) // Loop over temperatures (rows in files) { - for (int k=0; k < 8; ++k) // Loop over densities (columns in files) + for (k=0; k < n_density+1; ++k) // Loop over densities (columns in files) { fin >> number; - if ( k == 0 ) + fin >> comma; + if ( k == 0 && i == 0) { - parameters.log10_temperature_array[j] = number; + /* Since all files use the same temperature array, + * we only store the values from the first file*/ + parameters.log10_temperature_array.push_back(number); } - else + else if ( k > 0 ) { - parameters.log10_loss_rate_array[j][k][i] = number; - //[temperature_index][density_index][abundance_index] + loss_rate_1D.push_back(number); + // [abundance_index][temperature_index][density_index] } - fin >> comma; } + loss_rate_2D.push_back(loss_rate_1D); + loss_rate_1D.clear(); } fin.close(); fin.clear(); + + parameters.log10_loss_rate_array.push_back(loss_rate_2D); + loss_rate_2D.clear(); } + } double Loop::CalculateVelocity(double temperature_e, double temperature_i, double pressure_e) diff --git a/source/util/misc.cpp b/source/util/misc.cpp index 7f33df4..b404907 100644 --- a/source/util/misc.cpp +++ b/source/util/misc.cpp @@ -9,8 +9,7 @@ Miscellaneous functions for manipulating data int find_closest(double x, double array[], int array_length) { /* Traverses the whole array, so O(N) search. - * A binary search would be more efficient for long arrays: O(ln N). - */ + * A binary search would be more efficient for long arrays: O(ln N). */ int index = 0; double closest_value = array[0]; @@ -25,3 +24,22 @@ int find_closest(double x, double array[], int array_length) return index; } + +int find_closest(double x, std::vector array, int array_length) +{ + /* Traverses the whole array, so O(N) search. + * A binary search would be more efficient for long arrays: O(ln N). */ + int index = 0; + double closest_value = array[0]; + + for( int i=1; i < array_length; ++i ) + { + if( std::abs(x - closest_value) > std::abs(x - array[i]) ) + { + closest_value = array[i]; + index = i; + } + } + + return index; +} \ No newline at end of file diff --git a/source/util/misc.h b/source/util/misc.h index 4fac1cc..ab2fa05 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -8,10 +8,12 @@ #define MISC_H #include +#include -// Find the closest value in an array to the input x +// Find the closest value in an array (or vector) to the input x // Uses a simple search, O(N), so it's best for short arrays -// Returns the index of the array that's closest to value x +// Returns the index of the array (vector) that's closest to value x int find_closest(double x, double array[], int array_length); +int find_closest(double x, std::vector array, int array_length); #endif \ No newline at end of file From d4628f5360a020037c855211c35f9dfbbb550f7b Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Wed, 3 Jul 2024 14:55:16 -1000 Subject: [PATCH 21/24] add comment about initial conditions radiation --- source/helper.h | 3 +-- source/loop.cpp | 4 ++++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/source/helper.h b/source/helper.h index 9738e3d..308347d 100644 --- a/source/helper.h +++ b/source/helper.h @@ -88,9 +88,8 @@ struct Parameters { std::vector log10_density_array; std::vector abundance_array; /* The look-up table's radiative loss rate as a - * function of [temperature][density][abundance] */ + * function of [abundance][temperature][density] */ std::vector > > log10_loss_rate_array; - //double log10_loss_rate_array[101][7][7]; /* The density in the corona before upflows begin, used to calculate * the change in abundance factor */ double initial_density; diff --git a/source/loop.cpp b/source/loop.cpp index 16109ea..bbc8c91 100644 --- a/source/loop.cpp +++ b/source/loop.cpp @@ -141,6 +141,10 @@ state_type Loop::CalculateInitialConditions(void) if( parameters.use_lookup_table_losses ) { + /* The electron density has not been determined yet, so the look-up table + * radiative losses cannot be used for the initial conditions calculation. + * This check tells the call to CalculateC1 to use power-law losses for the + * initial conditions, instead. */ parameters.initial_radiation = true; } From 195ecdffcd6dc3cbda271fa48680a1e083e0dc38 Mon Sep 17 00:00:00 2001 From: Jeffrey Reep Date: Thu, 4 Jul 2024 11:15:12 -1000 Subject: [PATCH 22/24] remove overloaded find_closest --- source/util/misc.cpp | 20 -------------------- source/util/misc.h | 5 ++--- 2 files changed, 2 insertions(+), 23 deletions(-) diff --git a/source/util/misc.cpp b/source/util/misc.cpp index b404907..82b4aa5 100644 --- a/source/util/misc.cpp +++ b/source/util/misc.cpp @@ -5,26 +5,6 @@ Miscellaneous functions for manipulating data #include "misc.h" - -int find_closest(double x, double array[], int array_length) -{ - /* Traverses the whole array, so O(N) search. - * A binary search would be more efficient for long arrays: O(ln N). */ - int index = 0; - double closest_value = array[0]; - - for( int i=1; i < array_length; ++i ) - { - if( std::abs(x - closest_value) > std::abs(x - array[i]) ) - { - closest_value = array[i]; - index = i; - } - } - - return index; -} - int find_closest(double x, std::vector array, int array_length) { /* Traverses the whole array, so O(N) search. diff --git a/source/util/misc.h b/source/util/misc.h index ab2fa05..268622a 100644 --- a/source/util/misc.h +++ b/source/util/misc.h @@ -10,10 +10,9 @@ #include #include -// Find the closest value in an array (or vector) to the input x +// Find the closest value in a vector to the input x // Uses a simple search, O(N), so it's best for short arrays -// Returns the index of the array (vector) that's closest to value x -int find_closest(double x, double array[], int array_length); +// Returns the index of the vector that's closest to value x int find_closest(double x, std::vector array, int array_length); #endif \ No newline at end of file From 5ddf87afe2c906fd7859b7f022fbadb638fcc41b Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Thu, 4 Jul 2024 18:28:17 -0400 Subject: [PATCH 23/24] add tests for new radiative loss option --- tests/test_rad_loss.py | 57 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 tests/test_rad_loss.py diff --git a/tests/test_rad_loss.py b/tests/test_rad_loss.py new file mode 100644 index 0000000..8502208 --- /dev/null +++ b/tests/test_rad_loss.py @@ -0,0 +1,57 @@ +""" +Tests for the different kinds of radiative loss functions +""" +from collections import OrderedDict + +import pytest + +from .helpers import run_ebtelplusplus + + +@pytest.fixture(scope='module') +def base_config(): + base_config = { + 'total_time': 5e3, + 'tau': 1.0, + 'tau_max': 10.0, + 'loop_length': 4e9, + 'saturation_limit': 1/6, + 'force_single_fluid': False, + 'use_c1_loss_correction': True, + 'use_c1_grav_correction': True, + 'use_flux_limiting': True, + 'calculate_dem': False, + 'save_terms': False, + 'use_adaptive_solver': True, + 'adaptive_solver_error': 1e-6, + 'adaptive_solver_safety': 0.5, + 'c1_cond0': 2.0, + 'c1_rad0': 0.6, + 'helium_to_hydrogen_ratio': 0.075, + 'surface_gravity': 1.0, + 'heating': OrderedDict({ + 'partition': 1.0, + 'background': 1e-6, + 'events': [ + {'event': {'rise_start': 0.0, 'rise_end': 100.0, 'decay_start': 100.0, + 'decay_end': 200.0, 'magnitude': 0.1}} + ] + }), + } + return base_config + +@pytest.mark.parametrize('radiation', ['power_law', 'variable', 'coronal', 'photospheric']) +def test_rad_loss_options(base_config, radiation): + # Just a smoke test to make sure the new radiative loss options work + base_config['radiation'] = radiation + results = run_ebtelplusplus(base_config) + quantities = [ + 'electron_temperature', + 'ion_temperature', + 'density', + 'electron_pressure', + 'ion_pressure', + 'velocity' + ] + for q in quantities: + assert q in results From 488c154a71a4e08ca3631f1c97e0a5851d4c1835 Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Thu, 4 Jul 2024 18:28:37 -0400 Subject: [PATCH 24/24] document radiative loss options --- docs/mkdocs/configuration.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/mkdocs/configuration.md b/docs/mkdocs/configuration.md index 46cf077..31bbd9b 100644 --- a/docs/mkdocs/configuration.md +++ b/docs/mkdocs/configuration.md @@ -21,6 +21,7 @@ An ebtel++ run is configured by a single XML configuration file. The table below | **c1_rad0** | `float` | Nominal value of $c_1$ during radiative phase; see Eq. 16 of [Cargill et al. (2012a)][cargill_2012a] | | **helium_to_hydrogen_ratio** | `float` | Ratio of helium to hydrogen abundance; used in correction to ion mass, ion equation of state | | **surface_gravity** | `float` | Surface gravity in units of solar surface gravity; should be set to 1.0 unless using for extra-solar cases | +| **radiation** | `string` | The kind of radiative loss function to use. Must be either "power_law" (to use radiative losses of [Klimchuk et al. (2008)][klimchuk_2008]), "coronal" (to use radiative losses computed with coronal abundances), "photospheric" (to use radiative losses computed with photospheric abundances), or "variable" (to vary the radiative loss function from coronal to photospheric as a function of density and temperature) | ### Heating The time dependent heating is configured in a separate node. It includes the following parameters,