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 <iostream>
+#include <fstream>
+#include <sstream>
+#include <cmath>
+#include <string>
+#include <stdexcept>
+#include <array>
+#include <numbers>
+
+
+#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<double, 36> 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<double, 36> 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<std::string>("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.jmax; j++) { table.t[j] = exp10(tlo + tstp*j); }
+
+    for (i=0; i<table.imax; i++) { table.d[i] = exp10(dlo + dstp*i); }
+    
+    ifstream helm_table(filename);
+    if (!helm_table) {
+      LOG_ERROR(logger, "read_helm_table : Error opening file {}", filename);
+      throw std::runtime_error("Error opening file " + filename);
+    }
+    //read the Helmholtz free energy and its derivatives
+    for (j=0; j<table.jmax; j++) {
+      for (i=0; i<table.imax; i++){
+        getline(helm_table, data);
+        stringstream id(data);
+        id >> 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.jmax; j++) {
+      for (i=0; i<table.imax; i++){
+        getline(helm_table, data);
+        stringstream id(data);
+        id >> 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.jmax; j++) {
+      for (i=0; i<table.imax; i++){
+        getline(helm_table, data);
+        stringstream id(data);
+        id >> 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.jmax; j++) {
+      for (i=0; i<table.imax; i++){
+        getline(helm_table, data);
+        stringstream id(data);
+        id >> 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<table.jmax; j++) {
+      double dth         = table.t[j+1] - table.t[j];
+      double dt2         = dth * dth;
+      double dti         = 1.0/dth;
+      double dt2i        = 1.0/dt2;
+      table.dt_sav[j]   = dth;
+      table.dt2_sav[j]  = dt2;
+      table.dti_sav[j]  = dti;
+      table.dt2i_sav[j] = dt2i;
+    }
+
+
+    for (i=0; i<table.imax; i++) {
+      double dd          = table.d[i+1] - table.d[i];
+      double dd2         = dd * dd;
+      double ddi         = 1.0/dd;
+      table.dd_sav[i]   = dd;
+      table.dd2_sav[i]  = dd2;
+      table.ddi_sav[i]  = ddi;
+    }
+
+    table.loaded = true;
+
+    return table;
+  }
+
+
+  /***
+      this function evaluates the EOS components:
+        ion, radiation, electron-positron and Coulomb interaction
+      and returns the calculated quantities in the input 
+  ***/
+  EOS get_helm_EOS(EOSInput &q, const HELMTable &table) {
+    Config& config = Config::getInstance();
+    std::string logFile = config.get<std::string>("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<double, 36> 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 <iostream>
-#include <fstream>
-#include <sstream>
-#include <cmath>
-#include <format>
-
-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<JMAX; j++) { t[j] = exp10(tlo + tstp*j); }
-
-  for (i=0; i<IMAX; i++) { d[i] = exp10(dlo + dstp*i); }
-  
-  ifstream helm_table("helm_table.dat");
-  //read the Helmholtz free energy and its derivatives
-  for (j=0; j<JMAX; j++) {
-    for (i=0; i<IMAX; i++){
-      getline(helm_table, data);
-      stringstream id(data);
-      id >> 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<JMAX; j++) {
-    for (i=0; i<IMAX; i++){
-      getline(helm_table, data);
-      stringstream id(data);
-      id >> dpdf[i][j];
-      id >> dpdfd[i][j];
-      id >> dpdft[i][j];
-      id >> dpdfdt[i][j];
-    }
-  }
-  
-  //read the electron chemical potential
-  for (j=0; j<JMAX; j++) {
-    for (i=0; i<IMAX; i++){
-      getline(helm_table, data);
-      stringstream id(data);
-      id >> ef[i][j];
-      id >> efd[i][j];
-      id >> eft[i][j];
-      id >> efdt[i][j];
-    }
-  }
-  
-  //read the number density
-  for (j=0; j<JMAX; j++) {
-    for (i=0; i<IMAX; i++){
-      getline(helm_table, data);
-      stringstream id(data);
-      id >> 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<JMAX; j++) {
-    double dth          = t[j+1] - t[j];
-    double dt2         = dth * dth;
-    double dti         = 1.0/dth;
-    double dt2i        = 1.0/dt2;
-    dt_sav[j]   = dth;
-    dt2_sav[j]  = dt2;
-    dti_sav[j]  = dti;
-    dt2i_sav[j] = dt2i;
-  }
-
-
-  for (i=0; i<IMAX; i++) {
-    double dd          = d[i+1] - d[i];
-    double dd2         = dd * dd;
-    double ddi         = 1.0/dd;
-    dd_sav[i]   = dd;
-    dd2_sav[i]  = dd2;
-    ddi_sav[i]  = ddi;
-  }
-}
-
-
-//this class stores the information for one EOS call
-class EOS {
-public:
-  //input
-  double T, rho, abar, zbar;
-
-  //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
-  
-};
-
-/***
-    this function evaluates the EOS components:
-      ion, radiation, electron-positron and Coulomb interaction
-    and returns the calculated quantities in the input 
-***/
-void helmEOS(EOS &q) {
-  const double pi = 3.1415926535897932384;
-  const double amu = 1.66053878283e-24;
-  const double h = 6.6260689633e-27;
-  const double avo  = 6.0221417930e23;
-  const double kerg = 1.380650424e-16;
-  const double kergavo = kerg*avo;
-  const double clight  = 2.99792458e10;
-  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.;
-  double T, rho, s, fi[36], 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 dxnedd, dxnedt, dxneda, dxnedz;
-  double dxnidd, dxnida, detadz, detada;
-  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, JMAX-1)); // bounds
-  iat = max(0, min(iat, IMAX-1)); // bounds
-  iat = (log10(din)-dlo)*dstpi;
-
-  //using the indices, access the table:
-  fi[0]  = f[iat][jat];
-  fi[1]  = f[iat+1][jat];
-  fi[2]  = f[iat][jat+1];
-  fi[3]  = f[iat+1][jat+1];
-  fi[4]  = ft[iat][jat];
-  fi[5]  = ft[iat+1][jat];
-  fi[6]  = ft[iat][jat+1];
-  fi[7]  = ft[iat+1][jat+1];
-  fi[8]  = ftt[iat][jat];
-  fi[9]  = ftt[iat+1][jat];
-  fi[10] = ftt[iat][jat+1];
-  fi[11] = ftt[iat+1][jat+1];
-  fi[12] = fd[iat][jat];
-  fi[13] = fd[iat+1][jat];
-  fi[14] = fd[iat][jat+1];
-  fi[15] = fd[iat+1][jat+1];
-  fi[16] = fdd[iat][jat];
-  fi[17] = fdd[iat+1][jat];
-  fi[18] = fdd[iat][jat+1];
-  fi[19] = fdd[iat+1][jat+1];
-  fi[20] = fdt[iat][jat];
-  fi[21] = fdt[iat+1][jat];
-  fi[22] = fdt[iat][jat+1];
-  fi[23] = fdt[iat+1][jat+1];
-  fi[24] = fddt[iat][jat];
-  fi[25] = fddt[iat+1][jat]; 
-  fi[26] = fddt[iat][jat+1];
-  fi[27] = fddt[iat+1][jat+1];
-  fi[28] = fdtt[iat][jat];
-  fi[29] = fdtt[iat+1][jat];
-  fi[30] = fdtt[iat][jat+1];
-  fi[31] = fdtt[iat+1][jat+1];
-  fi[32] = fddtt[iat][jat];
-  fi[33] = fddtt[iat+1][jat];
-  fi[34] = fddtt[iat][jat+1];
-  fi[35] = fddtt[iat+1][jat+1];
-
-  double xt = (T - t[jat]) * dti_sav[jat];
-  double xd = (din - d[iat]) * ddi_sav[iat];
-  double mxt = 1-xt;
-  double mxd = 1-xd;
-
-  //density and temperature basis functions
-  double si0t =   psi0(xt);
-  double si1t =   psi1(xt)*dt_sav[jat];
-  double si2t =   psi2(xt)*dt2_sav[jat];
-
-  double si0mt =  psi0(mxt);
-  double si1mt = -psi1(mxt)*dt_sav[jat];
-  double si2mt =  psi2(mxt)*dt2_sav[jat];
-
-  double si0d =   psi0(xd);
-  double si1d =   psi1(xd)*dd_sav[iat];
-  double si2d =   psi2(xd)*dd2_sav[iat];
-
-  double si0md =  psi0(mxd);
-  double si1md = -psi1(mxd)*dd_sav[iat];
-  double si2md =  psi2(mxd)*dd2_sav[iat];
-
-  // derivatives of the weight functions
-  double dsi0t =   dpsi0(xt)*dti_sav[jat];
-  double dsi1t =   dpsi1(xt);
-  double dsi2t =   dpsi2(xt)*dt_sav[jat];
-
-  double dsi0mt = -dpsi0(mxt)*dti_sav[jat];
-  double dsi1mt =  dpsi1(mxt);
-  double dsi2mt = -dpsi2(mxt)*dt_sav[jat];
-
-  double dsi0d =   dpsi0(xd)*ddi_sav[iat];
-  double dsi1d =   dpsi1(xd);
-  double dsi2d =   dpsi2(xd)*dd_sav[iat];
-
-  double dsi0md = -dpsi0(mxd)*ddi_sav[iat];
-  double dsi1md =  dpsi1(mxd);
-  double dsi2md = -dpsi2(mxd)*dd_sav[iat];
-
-  // second derivatives of the weight functions
-  double ddsi0t =   ddpsi0(xt)*dt2i_sav[jat];
-  double ddsi1t =   ddpsi1(xt)*dti_sav[jat];
-  double ddsi2t =   ddpsi2(xt);
-
-  double ddsi0mt =  ddpsi0(mxt)*dt2i_sav[jat];
-  double ddsi1mt = -ddpsi1(mxt)*dti_sav[jat];
-  double ddsi2mt =  ddpsi2(mxt);
-
-  cout << " xt, mxt  = " << xt << " " << mxt << endl;  
-  cout << " dti_sav[jat] = " << dti_sav[jat] << endl;
-  cout << " dt2i_sav[jat] = " << dt2i_sav[jat] << endl;
-
-  cout << " ddpsi0(xt) = " << ddpsi0(xt) << endl;
-  cout << " ddpsi1(xt) = " << ddpsi1(xt) << endl;
-  cout << " ddpsi2(xt) = " << ddpsi2(xt) << endl;
-
-  cout << " ddpsi0(mxt) = " << ddpsi0(mxt) << endl;
-  cout << " ddpsi1(mxt) = " << ddpsi1(mxt) << endl;
-  cout << " ddpsi2(mxt) = " << ddpsi2(mxt) << endl;
-  
-  cout << " ddsi0t = " << ddsi0t << endl;
-  cout << " ddsi1t = " << ddsi1t << endl;
-  cout << " ddsi2t = " << ddsi2t << endl;
-
-  cout << " ddsi0mt = " << ddsi0mt << endl;
-  cout << " ddsi1mt = " << ddsi1mt << endl;
-  cout << " ddsi2mt = " << ddsi2mt << endl;
-
-
-  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);
-
-  cout << "df_tt = " << df_tt << endl;
-
-  
-  // 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)*dt_sav[jat];
-
-  si0mt  =  xpsi0(mxt);
-  si1mt  =  -xpsi1(mxt)*dt_sav[jat];
-
-  si0d   =  xpsi0(xd);
-  si1d   =  xpsi1(xd)*dd_sav[iat];
-
-  si0md  =  xpsi0(mxd);
-  si1md  =  -xpsi1(mxd)*dd_sav[iat];
-
-  
-  // derivatives of weight functions
-  dsi0t  = xdpsi0(xt)*dti_sav[jat];
-  dsi1t  = xdpsi1(xt);
-    
-  dsi0mt = -xdpsi0(mxt)*dti_sav[jat];
-  dsi1mt = xdpsi1(mxt);
-
-  dsi0d  = xdpsi0(xd)*ddi_sav[iat];
-  dsi1d  = xdpsi1(xd);
-
-  dsi0md = -xdpsi0(mxd)*ddi_sav[iat];
-  dsi1md = xdpsi1(mxd);
-
-  // pressure derivative
-  fi[0]  = dpdf[iat][jat];
-  fi[1]  = dpdf[iat+1][jat];
-  fi[2]  = dpdf[iat][jat+1];
-  fi[3]  = dpdf[iat+1][jat+1];
-  fi[4]  = dpdft[iat][jat];
-  fi[5]  = dpdft[iat+1][jat];
-  fi[6]  = dpdft[iat][jat+1];
-  fi[7]  = dpdft[iat+1][jat+1];
-  fi[8]  = dpdfd[iat][jat];
-  fi[9]  = dpdfd[iat+1][jat];
-  fi[10] = dpdfd[iat][jat+1];
-  fi[11] = dpdfd[iat+1][jat+1];
-  fi[12] = dpdfdt[iat][jat];
-  fi[13] = dpdfdt[iat+1][jat];
-  fi[14] = dpdfdt[iat][jat+1];
-  fi[15] = 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]  = ef[iat][jat];
-  fi[1]  = ef[iat+1][jat];
-  fi[2]  = ef[iat][jat+1];
-  fi[3]  = ef[iat+1][jat+1];
-  fi[4]  = eft[iat][jat];
-  fi[5]  = eft[iat+1][jat];
-  fi[6]  = eft[iat][jat+1];
-  fi[7]  = eft[iat+1][jat+1];
-  fi[8]  = efd[iat][jat];
-  fi[9] = efd[iat+1][jat];
-  fi[10] = efd[iat][jat+1];
-  fi[11] = efd[iat+1][jat+1];
-  fi[12] = efdt[iat][jat];
-  fi[13] = efdt[iat+1][jat];
-  fi[14] = efdt[iat][jat+1];
-  fi[15] = 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
-  detada = -x * din * ytot1;
-  detadz =  x * rho * ytot1;
-
-  // look in the number density table only once
-  fi[0]  = xf[iat][jat];
-  fi[1]  = xf[iat+1][jat];
-  fi[2]  = xf[iat][jat+1];
-  fi[3]  = xf[iat+1][jat+1];
-  fi[4]  = xft[iat][jat];
-  fi[5]  = xft[iat+1][jat];
-  fi[6]  = xft[iat][jat+1];
-  fi[7]  = xft[iat+1][jat+1];
-  fi[8]  = xfd[iat][jat];
-  fi[9] = xfd[iat+1][jat];
-  fi[10] = xfd[iat][jat+1];
-  fi[11] = xfd[iat+1][jat+1];
-  fi[12] = xfdt[iat][jat];
-  fi[13] = xfdt[iat+1][jat];
-  fi[14] = xfdt[iat][jat+1];
-  fi[15] = 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);
-  dxnedd = ye * x;
-
-  // derivative with respect to temperature
-  dxnedt   = h3(fi, dsi0t, dsi1t, dsi0mt, dsi1mt, si0d, si1d, si0md, si1md);
-
-  // derivative with respect to abar and zbar
-  dxneda = -x * din * ytot1;
-  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;
-
-  cout << "pgas = " << pgas << endl;
-  cout << "prad = " << prad << endl;
-  cout << "pion = " << pion << endl;
-  cout << "pele = " << pele << endl;
-  cout << "pcoul= " << pcoul << endl;
-
-  dpgasdd = dpiondd + dpepdd + dpcouldd;
-  dpgasdt = dpiondt + dpepdt + dpcouldt;
-  dpgasda = dpionda + dpepda + dpcoulda;
-  dpgasdz = dpiondz + dpepdz + dpcouldz;
-
-  cout << "dpgasdd = " << dpgasdd << endl;
-  cout << "dpraddd = " << dpraddd << endl;
-  cout << "dpiondd = " << dpiondd << endl;
-  cout << "dpeledd = " << dpepdd << endl;
-  cout << "dpcouldd= " << dpcouldd << endl;
-
-  degasdd = deiondd + deepdd + decouldd;
-  degasdt = deiondt + deepdt + decouldt;
-  degasda = deionda + deepda + decoulda;
-  degasdz = deiondz + deepdz + decouldz;
-
-  cout << "degasdd = " << degasdd << endl;
-  cout << "deraddd = " << deraddd << endl;
-  cout << "deiondd = " << deiondd << endl;
-  cout << "deeledd = " << deepdd << endl;
-  cout << "decouldd= " << decouldd << endl;
-
-  cout << "degasdt = " << degasdt << endl;
-  cout << "deraddt = " << deraddt << endl;
-  cout << "deiondt = " << deiondt << endl;
-  cout << "deeledt = " << deepdt << endl;
-  cout << "decouldt= " << decouldt << endl;
-
-  dsgasdd = dsiondd + dsepdd + dscouldd;
-  dsgasdt = dsiondt + dsepdt + dscouldt;
-  dsgasda = dsionda + dsepda + dscoulda;
-  dsgasdz = dsiondz + dsepdz + dscouldz;
-
-  cout << "sgas = " << sgas << endl;
-  cout << "srad = " << srad << endl;
-  cout << "sion = " << sion << endl;
-  cout << "sele = " << sele << endl;
-  cout << "scoul= " << scoul << endl;
-
-  cout << "dsgasdd = " << dsgasdd << endl;
-  cout << "dsraddd = " << dsraddd << endl;
-  cout << "dsiondd = " << dsiondd << endl;
-  cout << "dseledd = " << dsepdd << endl;
-  cout << "dscouldd= " << dscouldd << endl;
-
-  cout << "dsgasdt = " << dsgasdt << endl;
-  cout << "dsraddt = " << dsraddt << endl;
-  cout << "dsiondt = " << dsiondt << endl;
-  cout << "dseledt = " << dsepdt << endl;
-  cout << "dscouldt= " << dscouldt << endl;
-  
-  //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:
-  q.ptot = ptot; q.etot = etot; q.stot = stot;
-  q.pgas = pgas; q.egas = egas; q.sgas = sgas;
-  q.prad = prad; q.erad = erad; q.srad = srad;
-  q.chiT = chiT, q.chiRho = chiRho, q.csound = csound;
-  q.gamma1 = gamma1, q.gamma2 = gamma2; q.gamma3 = gamma3;
-  q.cV = cV, q.cP = cP, q.grad_ad = grad_ad;
-  q.etaele = etaele, q.xnefer = xnefer, q.ye = ye;
-
-  q.dpresdd = dpresdd;  q.dentrdd = dentrdd;  q.denerdd = denerdd;
-  q.dpresdt = dpresdt;  q.dentrdt = dentrdt;  q.denerdt = denerdt;
-  q.dpresda = dpresda;  q.dentrda = dentrda;  q.denerda = denerda;
-  q.dpresdz = dpresdz;  q.dentrdz = dentrdz;  q.denerdz = denerdz;
-    
-  // Maxwell relations
-  x = rho * rho;
-  q.dse = T * dentrdt/denerdt - 1.0;
-  q.dpe = (denerdd*x + T*dpresdt)/ptot - 1.0;
-  q.dsp = -dentrdd*x/dpresdt - 1.0;
-  
-}
-
-int main() {
-  int i;
-  double asum, zsum;
-  const int nel=3;
-  double xmass[nel], aion[nel], zion[nel];
-  EOS eos1;
-  
-  xmass[0] = 0.75; aion[0] =  1.0; zion[0] = 1.0;
-  xmass[1] = 0.23; aion[1] =  4.0; zion[1] = 2.0; 
-  xmass[2] = 0.02; aion[2] = 12.0; zion[2] = 6.0;
-
-  eos1.T = 1.0e8;
-  eos1.rho = 1.0e6;
-  
-  asum = 0.0;
-  zsum = 0.0;
-  for (i=0; i<nel; i++) {
-    asum += xmass[i]/aion[i];
-    zsum += xmass[i]*zion[i]/aion[i];
-  }
-  eos1.abar = 1.0/asum;
-  eos1.zbar = eos1.abar*zsum;
-  
-  read_helm_table();
-
-  helmEOS(eos1);
-
-  cout << "    T=" << eos1.T << " rho = " << eos1.rho << " ye = " << eos1.ye << endl;
-  cout << "helm: " << "       value " << "     d/drho" << "        d/dT " << "        d/da" << "        d/dz" << std::scientific << endl;
-  cout << "p tot=" << " " << eos1.ptot << " " << eos1.dpresdd << " " << eos1.dpresdt << " " << eos1.dpresda << " " << eos1.dpresdz << endl;
-  cout << "e tot=" << " " << eos1.etot << " " << eos1.denerdd << " " << eos1.denerdt << " " << eos1.denerda << " " << eos1.denerdz << endl;
-  cout << "s tot=" << " " << eos1.stot << " " << eos1.dentrdd << " " << eos1.dentrdt << " " << eos1.dentrda << " " << eos1.dentrdz << endl;
-
-  cout <<  " chiT   = " << format("{0:14.8e}", eos1.chiT) << endl;
-  cout <<  " chiRho = " << format("{0:14.8e}", eos1.chiRho) << endl;
-  cout <<  " etaele = " << format("{0:14.8e}", eos1.etaele) << endl;
-  cout <<  " xnefer = " << format("{0:14.8e}", eos1.xnefer) << endl;
-  cout <<  "    cV  = " << format("{0:14.8e}", eos1.cV) << endl;
-  cout <<  "    cP  = " << format("{0:14.8e}", eos1.cP) << endl;
-  cout <<  " gamma1 = " << format("{0:14.8e}", eos1.gamma1) << endl;
-  cout <<  " gamma2 = " << format("{0:14.8e}", eos1.gamma2) << endl;
-  cout <<  " gamma3 = " << format("{0:14.8e}", eos1.gamma3) << endl;
-  cout <<  " csound = " << format("{0:14.8e}", eos1.csound) << endl;
-
-  cout << " Maxwell Relations " << endl;
-  cout << " dse = " << format("{0:14.8e}", eos1.dse) << endl;
-  cout << " dpe = " << format("{0:14.8e}", eos1.dpe) << endl;
-  cout << " dsp = " << format("{0:14.8e}", eos1.dsp) << endl;
-  return 0;
-}
diff --git a/src/eos/public/helm.h b/src/eos/public/helm.h
new file mode 100644
index 0000000..0c9654a
--- /dev/null
+++ b/src/eos/public/helm.h
@@ -0,0 +1,359 @@
+// 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
+
+#ifndef HELM_H
+#define HELM_H
+
+#define IMAX 541
+#define JMAX 201
+
+#include <array>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <string>
+
+/**
+ * @brief 2D array template alias.
+ * @tparam T Type of the array elements.
+ * @tparam J Number of columns.
+ * @tparam I Number of rows.
+ */
+template <typename T, std::size_t J, std::size_t I>
+using array2D = std::array<std::array<T, J>, 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<double>::max(), tMax = std::numeric_limits<double>::lowest();
+            double dMin = std::numeric_limits<double>::max(), dMax = std::numeric_limits<double>::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