From 154004c8ca23fee0edb26eabf0a4d959c800cac7 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 5 Mar 2025 16:59:04 -0500 Subject: [PATCH] feat(eos): added helmholtz eos as module Aaron Dotter implimented a C++ version of Frank Timmes' fortran code helmholtz.f90. I have taken that and refactored it to work in the 4DSEE code style. This has mostly involved some light moving of stuff around. The biggest change is removing all globals, and reorienting memory to be heap allocated and contiguous. This is because there was too much memory being stack allocated. --- src/eos/meson.build | 23 + src/eos/private/helm.cpp | 823 ++++++++++++++++++++++++++++++++ src/eos/private/helmholtz.cpp | 865 ---------------------------------- src/eos/public/helm.h | 359 ++++++++++++++ 4 files changed, 1205 insertions(+), 865 deletions(-) create mode 100644 src/eos/meson.build create mode 100644 src/eos/private/helm.cpp delete mode 100644 src/eos/private/helmholtz.cpp create mode 100644 src/eos/public/helm.h diff --git a/src/eos/meson.build b/src/eos/meson.build new file mode 100644 index 0000000..ca3e18f --- /dev/null +++ b/src/eos/meson.build @@ -0,0 +1,23 @@ +# Define the library +eos_sources = files( + 'private/helm.cpp', +) + +eos_headers = files( + 'public/helm.h' +) + +# Define the libconst library so it can be linked against by other parts of the build system +libeos = static_library('eos', + eos_sources, + include_directories: include_directories('public'), + cpp_args: ['-fvisibility=default'], + dependencies: [const_dep, quill_dep, probe_dep, config_dep, mfem_dep], + install : true) + +eos_dep = declare_dependency( + include_directories: include_directories('public'), + link_with: libeos, +) +# Make headers accessible +install_headers(eos_headers, subdir : '4DSSE/eos') \ No newline at end of file diff --git a/src/eos/private/helm.cpp b/src/eos/private/helm.cpp new file mode 100644 index 0000000..faaed84 --- /dev/null +++ b/src/eos/private/helm.cpp @@ -0,0 +1,823 @@ +// this file contains a C++ conversion of Frank Timmes' fortran code +// helmholtz.f90, which implements the Helmholtz Free Energy EOS described +// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025 + +#include +#include +#include +#include +#include +#include +#include +#include + + +#include "helm.h" +#include "const.h" +#include "probe.h" +#include "config.h" +#include "quill/LogMacros.h" + +using namespace std; + +double** heap_allocate_contiguous_2D_memory(int rows, int cols) { + double **array = new double*[rows]; + array[0] = new double[rows * cols]; + + for (int i = 1; i < rows; i++) { + array[i] = array[0] + i * cols; + } + + return array; +} + +void heap_deallocate_contiguous_2D_memory(double **array) { + delete[] array[0]; + delete[] array; +} +// interpolating polynomila function definitions +namespace helmholtz { + double psi0(double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; } + + double dpsi0(double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); } + + double ddpsi0(double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); } + + double psi1(double z) { return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); } + + double dpsi1(double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; } + + double ddpsi1(double z) { return z * (z * (-60.0*z + 96.0) - 36.0); } + + double psi2(double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); } + + double dpsi2(double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); } + + double ddpsi2(double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); } + + double xpsi0(double z) { return z * z * (2.0 * z - 3.0) + 1.0; } + + double xdpsi0(double z) { return z * (6.0*z - 6.0); } + + double xpsi1(double z) { return z * ( z * (z - 2.0) + 1.0); } + + double xdpsi1(double z) { return z * (3.0 * z - 4.0) + 1.0; } + + double h3(std::array fi, double w0t, double w1t, double w0mt, double w1mt, + double w0d, double w1d, double w0md, double w1md) { + return fi[0] * w0d * w0t + fi[1] * w0md * w0t + + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt + + fi[4] * w0d * w1t + fi[5] * w0md * w1t + + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + + fi[8] * w1d * w0t + fi[9] * w1md * w0t + + fi[10] * w1d * w0mt + fi[11] * w1md * w0mt + + fi[12] * w1d * w1t + fi[13] * w1md * w1t + + fi[14] * w1d * w1mt + fi[15] * w1md * w1mt; + } + + double h5(std::array fi, double w0t, double w1t, double w2t,double w0mt, + double w1mt, double w2mt, double w0d, double w1d, double w2d, + double w0md, double w1md, double w2md) { + return fi[0] * w0d * w0t + fi[1] * w0md * w0t + + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt + + fi[4] * w0d * w1t + fi[5] * w0md * w1t + + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + + fi[8] * w0d * w2t + fi[ 9] * w0md * w2t + + fi[10] * w0d * w2mt + fi[11] * w0md * w2mt + + fi[12] * w1d * w0t + fi[13] * w1md * w0t + + fi[14] * w1d * w0mt + fi[15] * w1md * w0mt + + fi[16] * w2d * w0t + fi[17] * w2md * w0t + + fi[18] * w2d * w0mt + fi[19] * w2md * w0mt + + fi[20] * w1d * w1t + fi[21] * w1md * w1t + + fi[22] * w1d * w1mt + fi[23] * w1md * w1mt + + fi[24] * w2d * w1t + fi[25] * w2md * w1t + + fi[26] * w2d * w1mt + fi[27] * w2md * w1mt + + fi[28] * w1d * w2t + fi[29] * w1md * w2t + + fi[30] * w1d * w2mt + fi[31] * w1md * w2mt + + fi[32] * w2d * w2t + fi[33] * w2md * w2t + + fi[34] * w2d * w2mt + fi[35] * w2md * w2mt; + } + + // this function reads in the HELM table and stores in the above arrays + HELMTable read_helm_table(const std::string filename) { + Config& config = Config::getInstance(); + std::string logFile = config.get("EOS:Helm:LogFile", "log"); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger = logManager.getLogger(logFile); + LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename); + + HELMTable table; + string data; + int i, j; + + //set T and Rho (d) arrays + for (j=0; j> table.f[i][j]; + id >> table.fd[i][j]; + id >> table.ft[i][j]; + id >> table.fdd[i][j]; + id >> table.ftt[i][j]; + id >> table.fdt[i][j]; + id >> table.fddt[i][j]; + id >> table.fdtt[i][j]; + id >> table.fddtt[i][j]; + } + } + + //read the pressure derivative with density + for (j=0; j> table.dpdf[i][j]; + id >> table.dpdfd[i][j]; + id >> table.dpdft[i][j]; + id >> table.dpdfdt[i][j]; + } + } + + //read the electron chemical potential + for (j=0; j> table.ef[i][j]; + id >> table.efd[i][j]; + id >> table.eft[i][j]; + id >> table.efdt[i][j]; + } + } + + //read the number density + for (j=0; j> table.xf[i][j]; + id >> table.xfd[i][j]; + id >> table.xft[i][j]; + id >> table.xfdt[i][j]; + } + } + + helm_table.close(); //done reading + + // construct the temperature and density deltas and their inverses + for (j=0; j("EOS:Helm:LogFile", "log"); + Probe::LogManager& logManager = Probe::LogManager::getInstance(); + quill::Logger* logger = logManager.getLogger(logFile); + + Constants& constants = Constants::getInstance(); + const double pi = std::numbers::pi; + const double amu = constants.get("u").value; + const double h = constants.get("h").value; + const double avo = constants.get("N_a").value; + const double kerg = constants.get("kB").value; + const double clight = constants.get("c").value; + + const double kergavo = kerg*avo; + const double clight2 = clight * clight; + const double ssol = 5.6704e-5; + const double asol = 4 * ssol / clight; + const double asoli3 = asol / 3; + const double sioncon = (2.0 * pi * amu * kerg)/(h*h); + const double qe = 4.8032042712e-10; + const double esqu = qe * qe; + const double third = 1./3.; + std::array fi; + double T, rho, s, free, abar, zbar, ytot1, ye, ptot; + double prad, srad, erad, etot, stot; + double dpraddd, dpraddt, dpradda, dpraddz; + double deraddd, deraddt, deradda, deraddz; + double dsraddd, dsraddt, dsradda, dsraddz; + double dpiondd, dpiondt, dpionda, dpiondz; + double deiondd, deiondt, deionda, deiondz; + double dsiondd, dsiondt, dsionda, dsiondz; + double dpcouldd, dpcouldt, dpcoulda, dpcouldz; + double decouldd, decouldt, decoulda, decouldz; + double dscouldd, dscouldt, dscoulda, dscouldz; + double dpgasdd, dpgasdt, dpgasda, dpgasdz; + double degasdd, degasdt, degasda, degasdz; + double dsgasdd, dsgasdt, dsgasda, dsgasdz; + double dpepdd, dpepdt, dpepda, dpepdz; + double dsepdd, dsepdt, dsepda, dsepdz; + double deepdd, deepdt, deepda, deepdz; + double dxnidd, dxnida; + double dpresdd, dpresdt, dpresda, dpresdz; + double denerdd, denerdt, denerda, denerdz; + double dentrdd, dentrdt, dentrda, dentrdz; + double xni, sion, pion, eion, x, y, z; + double pgas, egas, sgas, pele, sele, eele, pcoul, scoul, ecoul; + int jat, iat; + + T = q.T; rho = q.rho; abar = q.abar; zbar = q.zbar; + + ytot1 = 1/abar; + ye = zbar * ytot1; + + double rhoi = 1/rho; + double Ti = 1/T; + double kT= kerg*T; + double kTinv = 1/kT; + + /*** radiation ***/ + prad = asoli3 * pow(T,4); + dpraddd = 0.; + dpraddt = 4 * prad * Ti; + dpradda = 0.; + dpraddz = 0.; + + erad = 3 * prad * rhoi; + deraddd = -erad * rhoi; + deraddt = 3 * dpraddt * rhoi; + deradda = 0.; + deraddz = 0.; + + srad = (prad*rhoi + erad)*Ti; + dsraddd = (dpraddd*rhoi - prad*rhoi*rhoi + deraddd)*Ti; + dsraddt = (dpraddt*rhoi + deraddt - srad)*Ti; + dsradda = 0.; + dsraddz = 0.; + + + /*** ions ***/ + xni = avo * ytot1 * rho; + dxnidd = avo * ytot1; + dxnida = -xni * ytot1; + + pion = xni * kT; + dpiondd = dxnidd * kT; + dpiondt = xni * kerg; + dpionda = dxnida * kT; + dpiondz = 0.; + + eion = 1.5 * pion * rhoi; + deiondd = (1.5 * dpiondd - eion)*rhoi; + deiondt = 1.5 * dpiondt*rhoi; + deionda = 1.5 * dpionda*rhoi; + deiondz = 0.; + + x = abar * abar * sqrt(abar) * rhoi/avo; + s = sioncon * T; + z = x * s * sqrt(s); + y = log(z); + sion = (pion * rhoi + eion) * Ti + kergavo * ytot1 * y; + + dsiondd = (dpiondd*rhoi - pion*rhoi*rhoi + deiondd)*Ti + - kergavo * rhoi * ytot1; + dsiondt = (dpiondt*rhoi + deiondt)*Ti - (pion*rhoi + eion) * Ti*Ti + + 1.5 * kergavo * Ti*ytot1; + x = avo*kerg/abar; + dsionda = (dpionda*rhoi + deionda)*Ti + kergavo*ytot1*ytot1* (2.5 - y); + dsiondz = 0.; + + + // electron-positron + // double xnem = xni * zbar; + double din = ye*rho; + + // hash the table location in t, rho: + jat = (log10(T) - tlo)*tstpi; + jat = max(0, min(jat, table.jmax-1)); // bounds + iat = max(0, min(iat, table.imax-1)); // bounds + iat = (log10(din)-dlo)*dstpi; + + //using the indices, access the table: + fi[0] = table.f[iat][jat]; + fi[1] = table.f[iat+1][jat]; + fi[2] = table.f[iat][jat+1]; + fi[3] = table.f[iat+1][jat+1]; + fi[4] = table.ft[iat][jat]; + fi[5] = table.ft[iat+1][jat]; + fi[6] = table.ft[iat][jat+1]; + fi[7] = table.ft[iat+1][jat+1]; + fi[8] = table.ftt[iat][jat]; + fi[9] = table.ftt[iat+1][jat]; + fi[10] = table.ftt[iat][jat+1]; + fi[11] = table.ftt[iat+1][jat+1]; + fi[12] = table.fd[iat][jat]; + fi[13] = table.fd[iat+1][jat]; + fi[14] = table.fd[iat][jat+1]; + fi[15] = table.fd[iat+1][jat+1]; + fi[16] = table.fdd[iat][jat]; + fi[17] = table.fdd[iat+1][jat]; + fi[18] = table.fdd[iat][jat+1]; + fi[19] = table.fdd[iat+1][jat+1]; + fi[20] = table.fdt[iat][jat]; + fi[21] = table.fdt[iat+1][jat]; + fi[22] = table.fdt[iat][jat+1]; + fi[23] = table.fdt[iat+1][jat+1]; + fi[24] = table.fddt[iat][jat]; + fi[25] = table.fddt[iat+1][jat]; + fi[26] = table.fddt[iat][jat+1]; + fi[27] = table.fddt[iat+1][jat+1]; + fi[28] = table.fdtt[iat][jat]; + fi[29] = table.fdtt[iat+1][jat]; + fi[30] = table.fdtt[iat][jat+1]; + fi[31] = table.fdtt[iat+1][jat+1]; + fi[32] = table.fddtt[iat][jat]; + fi[33] = table.fddtt[iat+1][jat]; + fi[34] = table.fddtt[iat][jat+1]; + fi[35] = table.fddtt[iat+1][jat+1]; + + double xt = (T - table.t[jat]) * table.dti_sav[jat]; + double xd = (din - table.d[iat]) * table.ddi_sav[iat]; + double mxt = 1-xt; + double mxd = 1-xd; + + //density and temperature basis functions + double si0t = psi0(xt); + double si1t = psi1(xt)*table.dt_sav[jat]; + double si2t = psi2(xt)*table.dt2_sav[jat]; + + double si0mt = psi0(mxt); + double si1mt = -psi1(mxt)*table.dt_sav[jat]; + double si2mt = psi2(mxt)*table.dt2_sav[jat]; + + double si0d = psi0(xd); + double si1d = psi1(xd)*table.dd_sav[iat]; + double si2d = psi2(xd)*table.dd2_sav[iat]; + + double si0md = psi0(mxd); + double si1md = -psi1(mxd)*table.dd_sav[iat]; + double si2md = psi2(mxd)*table.dd2_sav[iat]; + + // derivatives of the weight functions + double dsi0t = dpsi0(xt)*table.dti_sav[jat]; + double dsi1t = dpsi1(xt); + double dsi2t = dpsi2(xt)*table.dt_sav[jat]; + + double dsi0mt = -dpsi0(mxt)*table.dti_sav[jat]; + double dsi1mt = dpsi1(mxt); + double dsi2mt = -dpsi2(mxt)*table.dt_sav[jat]; + + double dsi0d = dpsi0(xd)*table.ddi_sav[iat]; + double dsi1d = dpsi1(xd); + double dsi2d = dpsi2(xd)*table.dd_sav[iat]; + + double dsi0md = -dpsi0(mxd)*table.ddi_sav[iat]; + double dsi1md = dpsi1(mxd); + double dsi2md = -dpsi2(mxd)*table.dd_sav[iat]; + + // second derivatives of the weight functions + double ddsi0t = ddpsi0(xt)*table.dt2i_sav[jat]; + double ddsi1t = ddpsi1(xt)*table.dti_sav[jat]; + double ddsi2t = ddpsi2(xt); + + double ddsi0mt = ddpsi0(mxt)*table.dt2i_sav[jat]; + double ddsi1mt = -ddpsi1(mxt)*table.dti_sav[jat]; + double ddsi2mt = ddpsi2(mxt); + + // Refactor these to use LOG_DEBUG(logger, message) + LOG_DEBUG(logger, " xt, mxt = {} {}", xt, mxt); + LOG_DEBUG(logger, " dti_sav[jat] = {}", table.dti_sav[jat]); + LOG_DEBUG(logger, " dt2i_sav[jat] = {}", table.dt2i_sav[jat]); + LOG_DEBUG(logger, " ddpsi0(xt) = {}", ddpsi0(xt)); + LOG_DEBUG(logger, " ddpsi1(xt) = {}", ddpsi1(xt)); + LOG_DEBUG(logger, " ddpsi2(xt) = {}", ddpsi2(xt)); + LOG_DEBUG(logger, " ddpsi0(mxt) = {}", ddpsi0(mxt)); + LOG_DEBUG(logger, " ddpsi1(mxt) = {}", ddpsi1(mxt)); + LOG_DEBUG(logger, " ddpsi2(mxt) = {}", ddpsi2(mxt)); + LOG_DEBUG(logger, " ddsi0t = {}", ddsi0t); + LOG_DEBUG(logger, " ddsi1t = {}", ddsi1t); + LOG_DEBUG(logger, " ddsi2t = {}", ddsi2t); + LOG_DEBUG(logger, " ddsi0mt = {}", ddsi0mt); + LOG_DEBUG(logger, " ddsi1mt = {}", ddsi1mt); + LOG_DEBUG(logger, " ddsi2mt = {}", ddsi2mt); + + + ddsi0t = 0.0; ddsi1t = 0.0; ddsi2t = 1.0; + ddsi0mt = 0.0; ddsi1mt = 0.0; ddsi2mt = 0.0; + + // free energy + free = h5(fi,si0t, si1t, si2t, si0mt, si1mt, si2mt, + si0d, si1d, si2d, si0md, si1md, si2md); + + // derivative with respect to density + double df_d = h5(fi, si0t, si1t, si2t, si0mt, si1mt, si2mt, + dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md); + + // derivative with respect to temperature + double df_t = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, + si0d, si1d, si2d, si0md, si1md, si2md); + + // second derivative with respect to temperature + double df_tt = h5(fi, ddsi0t, ddsi1t, ddsi2t, ddsi0mt, ddsi1mt, ddsi2mt, + si0d, si1d, si2d, si0md, si1md, si2md); + + LOG_DEBUG(logger, "df_tt = {}", df_tt); + + + // derivative with respect to temperature and density + double df_dt = h5(fi, dsi0t, dsi1t, dsi2t, dsi0mt, dsi1mt, dsi2mt, + dsi0d, dsi1d, dsi2d, dsi0md, dsi1md, dsi2md); + + + // pressure derivatives + si0t = xpsi0(xt); + si1t = xpsi1(xt)*table.dt_sav[jat]; + + si0mt = xpsi0(mxt); + si1mt = -xpsi1(mxt)*table.dt_sav[jat]; + + si0d = xpsi0(xd); + si1d = xpsi1(xd)*table.dd_sav[iat]; + + si0md = xpsi0(mxd); + si1md = -xpsi1(mxd)*table.dd_sav[iat]; + + + // derivatives of weight functions + dsi0t = xdpsi0(xt)*table.dti_sav[jat]; + dsi1t = xdpsi1(xt); + + dsi0mt = -xdpsi0(mxt)*table.dti_sav[jat]; + dsi1mt = xdpsi1(mxt); + + dsi0d = xdpsi0(xd)*table.ddi_sav[iat]; + dsi1d = xdpsi1(xd); + + dsi0md = -xdpsi0(mxd)*table.ddi_sav[iat]; + dsi1md = xdpsi1(mxd); + + // pressure derivative + fi[0] = table.dpdf[iat][jat]; + fi[1] = table.dpdf[iat+1][jat]; + fi[2] = table.dpdf[iat][jat+1]; + fi[3] = table.dpdf[iat+1][jat+1]; + fi[4] = table.dpdft[iat][jat]; + fi[5] = table.dpdft[iat+1][jat]; + fi[6] = table.dpdft[iat][jat+1]; + fi[7] = table.dpdft[iat+1][jat+1]; + fi[8] = table.dpdfd[iat][jat]; + fi[9] = table.dpdfd[iat+1][jat]; + fi[10] = table.dpdfd[iat][jat+1]; + fi[11] = table.dpdfd[iat+1][jat+1]; + fi[12] = table.dpdfdt[iat][jat]; + fi[13] = table.dpdfdt[iat+1][jat]; + fi[14] = table.dpdfdt[iat][jat+1]; + fi[15] = table.dpdfdt[iat+1][jat+1]; + + + // pressure derivative with density + dpepdd = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md); + dpepdd = max(ye * dpepdd,1e-30); + + + // look in the electron chemical potential table only once + fi[0] = table.ef[iat][jat]; + fi[1] = table.ef[iat+1][jat]; + fi[2] = table.ef[iat][jat+1]; + fi[3] = table.ef[iat+1][jat+1]; + fi[4] = table.eft[iat][jat]; + fi[5] = table.eft[iat+1][jat]; + fi[6] = table.eft[iat][jat+1]; + fi[7] = table.eft[iat+1][jat+1]; + fi[8] = table.efd[iat][jat]; + fi[9] = table.efd[iat+1][jat]; + fi[10] = table.efd[iat][jat+1]; + fi[11] = table.efd[iat+1][jat+1]; + fi[12] = table.efdt[iat][jat]; + fi[13] = table.efdt[iat+1][jat]; + fi[14] = table.efdt[iat][jat+1]; + fi[15] = table.efdt[iat+1][jat+1]; + + // electron chemical potential etaele + double etaele = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md); + + // derivative with respect to density + x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md); + // double detadd = ye * x; + + // derivative with respect to temperature + // double detadt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, + // si1d, si0md, si1md); + + // derivative with respect to abar and zbar + // double detada = -x * din * ytot1; + // double detadz = x * rho * ytot1; + + // look in the number density table only once + fi[0] = table.xf[iat][jat]; + fi[1] = table.xf[iat+1][jat]; + fi[2] = table.xf[iat][jat+1]; + fi[3] = table.xf[iat+1][jat+1]; + fi[4] = table.xft[iat][jat]; + fi[5] = table.xft[iat+1][jat]; + fi[6] = table.xft[iat][jat+1]; + fi[7] = table.xft[iat+1][jat+1]; + fi[8] = table.xfd[iat][jat]; + fi[9] = table.xfd[iat+1][jat]; + fi[10] = table.xfd[iat][jat+1]; + fi[11] = table.xfd[iat+1][jat+1]; + fi[12] = table.xfdt[iat][jat]; + fi[13] = table.xfdt[iat+1][jat]; + fi[14] = table.xfdt[iat][jat+1]; + fi[15] = table.xfdt[iat+1][jat+1]; + + // electron chemical potential etaele + double xnefer = h3(fi, si0t, si1t, si0mt, si1mt, si0d, si1d, si0md, si1md); + + // derivative with respect to density + x = h3(fi, si0t, si1t, si0mt, si1mt, dsi0d, dsi1d, dsi0md, dsi1md); + x = max(x,1e-30); + // double dxnedd = ye * x; + + // derivative with respect to temperature + // double dxnedt = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, si1d, si0md, si1md); + + // derivative with respect to abar and zbar + // double dxneda = -x * din * ytot1; + // double dxnedz = x * rho * ytot1; + + + // the desired electron-positron thermodynamic quantities + x = din * din; + pele = x * df_d; + dpepdt = x * df_dt; + //dpepdd is set above + s = dpepdd/ye - 2 * din * df_d; + dpepda = -ytot1 * (2 * pele + s * din); + dpepdz = rho * ytot1 * (2 * din * df_d + s); + + x = ye * ye; + sele = -df_t * ye; + dsepdt = -df_tt * ye; + dsepdd = -df_dt * x; + dsepda = ytot1 * (ye * df_dt * din - sele); + dsepdz = -ytot1 * (ye * df_dt * rho + df_t); + + eele = ye*free + T * sele; + deepdt = T * dsepdt; + deepdd = x * df_d + T * dsepdd; + deepda = -ye * ytot1 * (free + df_d * din) + T * dsepda; + deepdz = ytot1 * (free + ye * df_d * rho) + T * dsepdz; + + // coulomb + + const double a1 = -0.898004; + const double b1 = 0.96786; + const double c1 = 0.220703; + const double d1 = -0.86097; + const double e1 = 2.5269; + const double a2 = 0.29561; + const double b2 = 1.9885; + const double c2 = 0.288675; + + z = 4 * pi * third; + s = z * xni; + double dsdd = z * dxnidd; + double dsda = z * dxnida; + + double lami = pow(s,-1./3.); + double inv_lami = 1.0/lami; + z = -lami/3.; + double lamidd = z * dsdd/s; + double lamida = z * dsda/s; + + double plasg = zbar * zbar * esqu * kTinv * inv_lami; + z = -plasg * inv_lami; + double plasgdd = z * lamidd; + double plasgda = z * lamida; + double plasgdt = -plasg*kTinv * kerg; + double plasgdz = 2 * plasg/zbar; + + // yakovlev & shalybkov 1989 equations 82, 85, 86, 87 + if (plasg >= 1.0) { + x = pow(plasg,0.25); + y = avo * ytot1 * kerg; + ecoul = y * T * (a1*plasg + b1*x + c1/x + d1); + pcoul = third * rho * ecoul; + scoul = -y * (3*b1*x - 5*c1/x + d1 * (log(plasg) - 1) - e1); + + y = avo * ytot1 * kT * (a1 + 0.25/plasg*(b1*x - c1/x)); + decouldd = y * plasgdd; + decouldt = y * plasgdt + ecoul/T; + decoulda = y * plasgda - ecoul/abar; + decouldz = y * plasgdz; + + y = third * rho; + dpcouldd = third * ecoul + y * decouldd; + dpcouldt = y * decouldt; + dpcoulda = y * decoulda; + dpcouldz = y * decouldz; + + + y = -avo * kerg / (abar*plasg) * (0.75*b1*x+1.25*c1/x+d1); + dscouldd = y * plasgdd; + dscouldt = y * plasgdt; + dscoulda = y * plasgda - scoul/abar; + dscouldz = y * plasgdz; + } + + // yakovlev & shalybkov 1989 equations 102, 103, 104 + else { + x = plasg*sqrt(plasg); + y = pow(plasg,b2); + z = c2 * x - third * a2 * y; + pcoul = -pion * z; + ecoul = 3 * pcoul/rho; + scoul = -(avo/abar) * kerg * (c2*x -a2*(b2-1)/b2*y); + + s = 1.5*c2*x/plasg - third*a2*b2*y/plasg; + dpcouldd = -dpiondd*z - pion*s*plasgdd; + dpcouldt = -dpiondt*z - pion*s*plasgdt; + dpcoulda = -dpionda*z - pion*s*plasgda; + dpcouldz = -dpiondz*z - pion*s*plasgdz; + + s = 3/rho; + decouldd = s * dpcouldd - ecoul/rho; + decouldt = s * dpcouldt; + decoulda = s * dpcoulda; + decouldz = s * dpcouldz; + + s = -avo*kerg/(abar*plasg)*(1.5*c2*x - a2*(b2-1)*y); + dscouldd = s * plasgdd; + dscouldt = s * plasgdt; + dscoulda = s * plasgda - scoul/abar; + dscouldz = s * plasgdz; + } + + // "bomb proof" + x = prad + pion + pele + pcoul; + y = erad + eion + eele + ecoul; + z = srad + sion + sele + scoul; + + if(x <= 0 || y <= 0) { + // zero out coulomb terms + pcoul = 0; dpcouldd=0; dpcouldt=0; dpcoulda=0; dpcouldz=0; + ecoul = 0; decouldd=0; decouldt=0; decoulda=0; decouldz=0; + scoul = 0; dscouldd=0; dscouldt=0; dscoulda=0; dscouldz=0; + } + + + // gas subtotals + pgas = pion + pele + pcoul; + egas = eion + eele + ecoul; + sgas = sion + sele + scoul; + + LOG_DEBUG(logger, "pgas = {}", pgas); + LOG_DEBUG(logger, "prad = {}", prad); + LOG_DEBUG(logger, "pion = {}", pion); + LOG_DEBUG(logger, "pele = {}", pele); + LOG_DEBUG(logger, "pcoul = {}", pcoul); + + dpgasdd = dpiondd + dpepdd + dpcouldd; + dpgasdt = dpiondt + dpepdt + dpcouldt; + dpgasda = dpionda + dpepda + dpcoulda; + dpgasdz = dpiondz + dpepdz + dpcouldz; + + LOG_DEBUG(logger, "dpgasdd = {}", dpgasdd); + LOG_DEBUG(logger, "dpraddd = {}", dpraddd); + LOG_DEBUG(logger, "dpiondd = {}", dpiondd); + LOG_DEBUG(logger, "dpeledd = {}", dpepdd); + LOG_DEBUG(logger, "dpcouldd = {}", dpcouldd); + + degasdd = deiondd + deepdd + decouldd; + degasdt = deiondt + deepdt + decouldt; + degasda = deionda + deepda + decoulda; + degasdz = deiondz + deepdz + decouldz; + + LOG_DEBUG(logger, "degasdd = {}", degasdd); + LOG_DEBUG(logger, "deraddd = {}", deraddd); + LOG_DEBUG(logger, "deiondd = {}", deiondd); + LOG_DEBUG(logger, "deeledd = {}", deepdd); + LOG_DEBUG(logger, "decouldd = {}", decouldd); + + LOG_DEBUG(logger, "degasdt = {}", degasdt); + LOG_DEBUG(logger, "deraddt = {}", deraddt); + LOG_DEBUG(logger, "deiondt = {}", deiondt); + LOG_DEBUG(logger, "deeledt = {}", deepdt); + LOG_DEBUG(logger, "decouldt = {}", decouldt); + + dsgasdd = dsiondd + dsepdd + dscouldd; + dsgasdt = dsiondt + dsepdt + dscouldt; + dsgasda = dsionda + dsepda + dscoulda; + dsgasdz = dsiondz + dsepdz + dscouldz; + + + LOG_DEBUG(logger, "sgas = {}", sgas); + LOG_DEBUG(logger, "srad = {}", srad); + LOG_DEBUG(logger, "sion = {}", sion); + LOG_DEBUG(logger, "sele = {}", sele); + LOG_DEBUG(logger, "scoul = {}", scoul); + + LOG_DEBUG(logger, "dsgasdd = {}", dsgasdd); + LOG_DEBUG(logger, "dsraddd = {}", dsraddd); + LOG_DEBUG(logger, "dsiondd = {}", dsiondd); + LOG_DEBUG(logger, "dseledd = {}", dsepdd); + LOG_DEBUG(logger, "dscouldd = {}", dscouldd); + + LOG_DEBUG(logger, "dsgasdt = {}", dsgasdt); + LOG_DEBUG(logger, "dsraddt = {}", dsraddt); + LOG_DEBUG(logger, "dsiondt = {}", dsiondt); + LOG_DEBUG(logger, "dseledt = {}", dsepdt); + LOG_DEBUG(logger, "dscouldt = {}", dscouldt); + + //total + ptot = prad + pgas; + stot = srad + sgas; + etot = erad + egas; + + dpresdd = dpraddd + dpgasdd; + dpresdt = dpraddt + dpgasdt; + dpresda = dpradda + dpgasda; + dpresdz = dpraddz + dpgasdz; + + denerdd = deraddd + degasdd; + denerdt = deraddt + degasdt; + denerda = deradda + degasda; + denerdz = deraddz + degasdz; + + dentrdd = dsraddd + dsgasdd; + dentrdt = dsraddt + dsgasdt; + dentrda = dsradda + dsgasda; + dentrdz = dsraddz + dsgasdz; + + // derived quantities + double zz = ptot * rhoi; + double zzi = rho/ptot; + double chit = (T/ptot)*dpresdt; + double chirho = zzi*dpresdd; + double cv = denerdt; + x = zz * chit/(T*cv); + double gamma3 = x + 1.0; + double gamma1 = chit * x + chirho; + double grad_ad = x/gamma1; + double gamma2 = 1.0/(1.0 - grad_ad); + double cp = cv * gamma1/chirho; + z = 1.0 + (etot + clight2)*zzi; + double csound = clight * sqrt(gamma1/z); + + // package in q: + EOS eos; + eos.ptot = ptot; eos.etot = etot; eos.stot = stot; + eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas; + eos.prad = prad; eos.erad = erad; eos.srad = srad; + eos.chiT = chit, eos.chiRho = chirho, eos.csound = csound; + eos.gamma1 = gamma1, eos.gamma2 = gamma2; eos.gamma3 = gamma3; + eos.cV = cv, eos.cP = cp, eos.grad_ad = grad_ad; + eos.etaele = etaele, eos.xnefer = xnefer, eos.ye = ye; + + eos.dpresdd = dpresdd; eos.dentrdd = dentrdd; eos.denerdd = denerdd; + eos.dpresdt = dpresdt; eos.dentrdt = dentrdt; eos.denerdt = denerdt; + eos.dpresda = dpresda; eos.dentrda = dentrda; eos.denerda = denerda; + eos.dpresdz = dpresdz; eos.dentrdz = dentrdz; eos.denerdz = denerdz; + + // maxwell relations + x = rho * rho; + eos.dse = T * dentrdt/denerdt - 1.0; + eos.dpe = (denerdd*x + T*dpresdt)/ptot - 1.0; + eos.dsp = -dentrdd*x/dpresdt - 1.0; + + return eos; + // TODO: In future this might be made more preformant by avoiding the copy. + } +} \ No newline at end of file diff --git a/src/eos/private/helmholtz.cpp b/src/eos/private/helmholtz.cpp deleted file mode 100644 index 54dbf65..0000000 --- a/src/eos/private/helmholtz.cpp +++ /dev/null @@ -1,865 +0,0 @@ -// this file contains a C++ conversion of Frank Timmes' fortran code -// helmholtz.f90, which implements the Helmholtz Free Energy EOS described -// by Timmes & Swesty (2000) doi:10.1086/313304 -- Aaron Dotter 2025 - -#include -#include -#include -#include -#include - -using namespace std; - -// interpolating polynomila function definitions -double psi0(double z) {return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0;} - -double dpsi0(double z) {return z*z * ( z * (-30.0 * z + 60.0) - 30.0);} - -double ddpsi0(double z) {return z * ( z * (-120.0 * z + 180.0) - 60.0);} - -double psi1(double z) {return z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0);} - -double dpsi1(double z) {return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0;} - -double ddpsi1(double z) {return z * (z * (-60.0*z + 96.0) - 36.0);} - -double psi2(double z) {return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0);} - -double dpsi2(double z) {return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0);} - -double ddpsi2(double z) {return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.);} - -double xpsi0(double z) {return z * z * (2.0 * z - 3.0) + 1.0;} - -double xdpsi0(double z) {return z * (6.0*z - 6.0);} - -double xpsi1(double z) {return z * ( z * (z - 2.0) + 1.0);} - -double xdpsi1(double z) {return z * (3.0 * z - 4.0) + 1.0;} - -double h3(double fi[36], double w0t, double w1t, double w0mt, double w1mt, - double w0d, double w1d, double w0md, double w1md) { - return fi[0] * w0d * w0t + fi[1] * w0md * w0t - + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt - + fi[4] * w0d * w1t + fi[5] * w0md * w1t - + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt - + fi[8] * w1d * w0t + fi[9] * w1md * w0t - + fi[10] * w1d * w0mt + fi[11] * w1md * w0mt - + fi[12] * w1d * w1t + fi[13] * w1md * w1t - + fi[14] * w1d * w1mt + fi[15] * w1md * w1mt; - } - -double h5(double fi[36], double w0t, double w1t, double w2t,double w0mt, - double w1mt, double w2mt, double w0d, double w1d, double w2d, - double w0md, double w1md, double w2md) { - return fi[0] * w0d * w0t + fi[1] * w0md * w0t - + fi[2] * w0d * w0mt + fi[3] * w0md * w0mt - + fi[4] * w0d * w1t + fi[5] * w0md * w1t - + fi[6] * w0d * w1mt + fi[7] * w0md * w1mt - + fi[8] * w0d * w2t + fi[ 9] * w0md * w2t - + fi[10] * w0d * w2mt + fi[11] * w0md * w2mt - + fi[12] * w1d * w0t + fi[13] * w1md * w0t - + fi[14] * w1d * w0mt + fi[15] * w1md * w0mt - + fi[16] * w2d * w0t + fi[17] * w2md * w0t - + fi[18] * w2d * w0mt + fi[19] * w2md * w0mt - + fi[20] * w1d * w1t + fi[21] * w1md * w1t - + fi[22] * w1d * w1mt + fi[23] * w1md * w1mt - + fi[24] * w2d * w1t + fi[25] * w2md * w1t - + fi[26] * w2d * w1mt + fi[27] * w2md * w1mt - + fi[28] * w1d * w2t + fi[29] * w1md * w2t - + fi[30] * w1d * w2mt + fi[31] * w1md * w2mt - + fi[32] * w2d * w2t + fi[33] * w2md * w2t - + fi[34] * w2d * w2mt + fi[35] * w2md * w2mt; - } - -// some global quantities including table size, limits, and deltas -const int IMAX=541, JMAX=201; -const double tlo = 3.0, thi = 13.0; -const double dlo = -12.0, dhi = 15.0; -const double tstp = (thi - tlo)/(JMAX-1), tstpi=1/tstp; -const double dstp = (dhi-dlo)/(IMAX-1), dstpi = 1/dstp; -double t[JMAX], d[IMAX]; -double f[IMAX][JMAX], fd[IMAX][JMAX], ft[IMAX][JMAX]; -double fdd[IMAX][JMAX], ftt[IMAX][JMAX], fdt[IMAX][JMAX]; -double fddt[IMAX][JMAX], fdtt[IMAX][JMAX], fddtt[IMAX][JMAX]; -double dpdf[IMAX][JMAX], dpdfd[IMAX][JMAX], dpdft[IMAX][JMAX]; -double dpdfdt[IMAX][JMAX], ef[IMAX][JMAX], efd[IMAX][JMAX]; -double eft[IMAX][JMAX], efdt[IMAX][JMAX], xf[IMAX][JMAX]; -double xfd[IMAX][JMAX], xft[IMAX][JMAX], xfdt[IMAX][JMAX]; -double dt_sav[JMAX], dt2_sav[JMAX], dti_sav[JMAX], dt2i_sav[JMAX]; -double dd_sav[IMAX], dd2_sav[IMAX], ddi_sav[IMAX]; - -// this function reads in the HELM table and stores in the above arrays -void read_helm_table() { - string data; - int i, j; - - //set T and Rho (d) arrays - for (j=0; j> f[i][j]; - id >> fd[i][j]; - id >> ft[i][j]; - id >> fdd[i][j]; - id >> ftt[i][j]; - id >> fdt[i][j]; - id >> fddt[i][j]; - id >> fdtt[i][j]; - id >> fddtt[i][j]; - } - } - - //read the pressure derivative with density - for (j=0; j> dpdf[i][j]; - id >> dpdfd[i][j]; - id >> dpdft[i][j]; - id >> dpdfdt[i][j]; - } - } - - //read the electron chemical potential - for (j=0; j> ef[i][j]; - id >> efd[i][j]; - id >> eft[i][j]; - id >> efdt[i][j]; - } - } - - //read the number density - for (j=0; j> xf[i][j]; - id >> xfd[i][j]; - id >> xft[i][j]; - id >> xfdt[i][j]; - } - } - - helm_table.close(); //done reading - - // construct the temperature and density deltas and their inverses - for (j=0; j +#include +#include +#include +#include + +/** + * @brief 2D array template alias. + * @tparam T Type of the array elements. + * @tparam J Number of columns. + * @tparam I Number of rows. + */ +template +using array2D = std::array, I>; + +double** heap_allocate_contiguous_2D_memory(int rows, int cols); + +void heap_deallocate_contiguous_2D_memory(double **array); + +namespace helmholtz +{ + const double tlo = 3.0, thi = 13.0; + const double dlo = -12.0, dhi = 15.0; + const double tstp = (thi - tlo) / (JMAX - 1), tstpi = 1 / tstp; + const double dstp = (dhi - dlo) / (IMAX - 1), dstpi = 1 / dstp; + + /** + * @struct HELMTable + * @brief Structure to hold the Helmholtz EOS table data. + */ + struct HELMTable + { + bool loaded = false; + const int imax = IMAX, jmax = JMAX; + double t[JMAX], d[IMAX]; + double** f; + double** fd; + double** ft; + double** fdd; + double** ftt; + double** fdt; + double** fddt; + double** fdtt; + double** fddtt; + double** dpdf; + double** dpdfd; + double** dpdft; + double** dpdfdt; + double** ef; + double** efd; + double** eft; + double** efdt; + double** xf; + double** xfd; + double** xft; + double** xfdt; + + double dt_sav[JMAX], dt2_sav[JMAX], dti_sav[JMAX], dt2i_sav[JMAX]; + double dd_sav[IMAX], dd2_sav[IMAX], ddi_sav[IMAX]; + + // Constructor + HELMTable() { + f = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fd = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + ft = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fdd = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + ftt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fddt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fdtt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + fddtt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + dpdf = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + dpdfd = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + dpdft = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + dpdfdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + ef = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + efd = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + eft = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + efdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + xf = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + xfd = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + xft = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + xfdt = heap_allocate_contiguous_2D_memory(IMAX, JMAX); + } + + // Destructor + ~HELMTable() { + heap_deallocate_contiguous_2D_memory(f); + heap_deallocate_contiguous_2D_memory(fd); + heap_deallocate_contiguous_2D_memory(ft); + heap_deallocate_contiguous_2D_memory(fdd); + heap_deallocate_contiguous_2D_memory(ftt); + heap_deallocate_contiguous_2D_memory(fdt); + heap_deallocate_contiguous_2D_memory(fddt); + heap_deallocate_contiguous_2D_memory(fdtt); + heap_deallocate_contiguous_2D_memory(fddtt); + heap_deallocate_contiguous_2D_memory(dpdf); + heap_deallocate_contiguous_2D_memory(dpdfd); + heap_deallocate_contiguous_2D_memory(dpdft); + heap_deallocate_contiguous_2D_memory(dpdfdt); + heap_deallocate_contiguous_2D_memory(ef); + heap_deallocate_contiguous_2D_memory(efd); + heap_deallocate_contiguous_2D_memory(eft); + heap_deallocate_contiguous_2D_memory(efdt); + heap_deallocate_contiguous_2D_memory(xf); + heap_deallocate_contiguous_2D_memory(xfd); + heap_deallocate_contiguous_2D_memory(xft); + heap_deallocate_contiguous_2D_memory(xfdt); + } + + friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMTable& table) { + if (!table.loaded) { + os << "HELMTable not loaded\n"; + return os; + } + double tMin = std::numeric_limits::max(), tMax = std::numeric_limits::lowest(); + double dMin = std::numeric_limits::max(), dMax = std::numeric_limits::lowest(); + + for (const auto& temp : table.t) { + if (temp < tMin) tMin = temp; + if (temp > tMax) tMax = temp; + } + + for (const auto& dens : table.d) { + if (dens < dMin) dMin = dens; + if (dens > dMax) dMax = dens; + } + + os << "HELMTable Data:\n"; + os << " imax: " << table.imax << ", jmax: " << table.jmax << "\n"; + os << " Temperature Range: [" << tMin << ", " << tMax << "]\n"; + os << " Density Range: [" << dMin << ", " << dMax << "]\n"; + + return os; + } + + + }; + + /** + * @struct EOSInput + * @brief Structure to hold the input parameters for the EOS calculation. + */ + struct EOSInput + { + double T; ///< Temperature. + double rho; ///< Density. + double abar; ///< Mean atomic mass. + double zbar; ///< Mean atomic number. + + friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOSInput& eosInput) { + os << "EOSInput Data:\n"; + os << " Temperature: " << eosInput.T << "\n"; + os << " Density: " << eosInput.rho << "\n"; + os << " Mean Atomic Mass: " << eosInput.abar << "\n"; + os << " Mean Atomic Number: " << eosInput.zbar << "\n"; + + return os; + } + }; + + /** + * @struct EOS + * @brief Structure to hold the output parameters and derivatives of the EOS calculation. + */ + struct EOS + { + // output + double ye, etaele, xnefer; // + double ptot, pgas, prad; // pressure + double etot, egas, erad; // energy + double stot, sgas, srad; // entropy + + // derivatives + double dpresdd, dpresdt, dpresda, dpresdz; + double dentrdd, dentrdt, dentrda, dentrdz; + double denerdd, denerdt, denerda, denerdz; + + // these could become functions: + double chiT, chiRho, csound, grad_ad; // derived quantities + double gamma1, gamma2, gamma3, cV, cP; // derived quantities + double dse, dpe, dsp; // Maxwell relations + + friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOS& eos) { + os << "EOS Data:\n" << std::setw(20) << std::left; + os << " Electron Fraction: " << eos.ye << "\n"; + os << " Electron Chemical Potential: " << eos.etaele << "\n"; + os << " Electron Number Density: " << eos.xnefer << "\n"; + os << " Total Pressure: " << eos.ptot << "\n"; + os << " Gas Pressure: " << eos.pgas << "\n"; + os << " Radiation Pressure: " << eos.prad << "\n"; + os << " Total Energy: " << eos.etot << "\n"; + os << " Gas Energy: " << eos.egas << "\n"; + os << " Radiation Energy: " << eos.erad << "\n"; + os << " Total Entropy: " << eos.stot << "\n"; + os << " Gas Entropy: " << eos.sgas << "\n"; + os << " Radiation Entropy: " << eos.srad; + return os; + } + }; + + // interpolating polynomial function definitions + + /** + * @brief Interpolating polynomial function psi0. + * @param z Input value. + * @return Result of the polynomial. + */ + double psi0(double z); + + /** + * @brief Derivative of the interpolating polynomial function psi0. + * @param z Input value. + * @return Derivative of the polynomial. + */ + double dpsi0(double z); + + /** + * @brief Second derivative of the interpolating polynomial function psi0. + * @param z Input value. + * @return Second derivative of the polynomial. + */ + double ddpsi0(double z); + + /** + * @brief Interpolating polynomial function psi1. + * @param z Input value. + * @return Result of the polynomial. + */ + double psi1(double z); + + /** + * @brief Derivative of the interpolating polynomial function psi1. + * @param z Input value. + * @return Derivative of the polynomial. + */ + double dpsi1(double z); + + /** + * @brief Second derivative of the interpolating polynomial function psi1. + * @param z Input value. + * @return Second derivative of the polynomial. + */ + double ddpsi1(double z); + + /** + * @brief Interpolating polynomial function psi2. + * @param z Input value. + * @return Result of the polynomial. + */ + double psi2(double z); + + /** + * @brief Derivative of the interpolating polynomial function psi2. + * @param z Input value. + * @return Derivative of the polynomial. + */ + double dpsi2(double z); + + /** + * @brief Second derivative of the interpolating polynomial function psi2. + * @param z Input value. + * @return Second derivative of the polynomial. + */ + double ddpsi2(double z); + + /** + * @brief Interpolating polynomial function xpsi0. + * @param z Input value. + * @return Result of the polynomial. + */ + double xpsi0(double z); + + /** + * @brief Derivative of the interpolating polynomial function xpsi0. + * @param z Input value. + * @return Derivative of the polynomial. + */ + double xdpsi0(double z); + + /** + * @brief Interpolating polynomial function xpsi1. + * @param z Input value. + * @return Result of the polynomial. + */ + double xpsi1(double z); + + /** + * @brief Derivative of the interpolating polynomial function xpsi1. + * @param z Input value. + * @return Derivative of the polynomial. + */ + double xdpsi1(double z); + + /** + * @brief Interpolating polynomial function h3. + * @param fi Array of coefficients. + * @param w0t Weight 0 for temperature. + * @param w1t Weight 1 for temperature. + * @param w0mt Weight 0 for temperature (minus). + * @param w1mt Weight 1 for temperature (minus). + * @param w0d Weight 0 for density. + * @param w1d Weight 1 for density. + * @param w0md Weight 0 for density (minus). + * @param w1md Weight 1 for density (minus). + * @return Result of the polynomial. + */ + double h3(double fi[36], double w0t, double w1t, double w0mt, double w1mt, + double w0d, double w1d, double w0md, double w1md); + + /** + * @brief Interpolating polynomial function h5. + * @param fi Array of coefficients. + * @param w0t Weight 0 for temperature. + * @param w1t Weight 1 for temperature. + * @param w2t Weight 2 for temperature. + * @param w0mt Weight 0 for temperature (minus). + * @param w1mt Weight 1 for temperature (minus). + * @param w2mt Weight 2 for temperature (minus). + * @param w0d Weight 0 for density. + * @param w1d Weight 1 for density. + * @param w2d Weight 2 for density. + * @param w0md Weight 0 for density (minus). + * @param w1md Weight 1 for density (minus). + * @param w2md Weight 2 for density (minus). + * @return Result of the polynomial. + */ + double h5(double fi[36], double w0t, double w1t, double w2t, double w0mt, + double w1mt, double w2mt, double w0d, double w1d, double w2d, + double w0md, double w1md, double w2md); + + /** + * @brief Read the Helmholtz EOS table from a file. + * @param filename Path to the file containing the table. + * @return HELMTable structure containing the table data. + */ + HELMTable read_helm_table(const std::string filename); + + /** + * @brief Calculate the Helmholtz EOS components. + * @param q EOSInput parameters for the EOS calculation. + * @param table HELMTable structure containing the table data. + * @return EOS structure containing the calculated quantities. + */ + EOS get_helm_EOS(EOSInput &q, const HELMTable &table); + +} + +#endif // HELM_H \ No newline at end of file