diff --git a/src/eos/private/helmholtz.cpp b/src/eos/private/helmholtz.cpp
new file mode 100644
index 0000000..54dbf65
--- /dev/null
+++ b/src/eos/private/helmholtz.cpp
@@ -0,0 +1,865 @@
+// 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;
+}