/* ***********************************************************************
//
//   Copyright (C) 2025 -- The 4D-STAR Collaboration
//   File Author: Emily Boudreaux
//   Last Modified: March 17, 2025
//
//   4DSSE is free software; you can use it and/or modify
//   it under the terms and restrictions the GNU General Library Public
//   License version 3 (GPLv3) as published by the Free Software Foundation.
//
//   4DSSE is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//   See the GNU Library General Public License for more details.
//
//   You should have received a copy of the GNU Library General Public License
//   along with this software; if not, write to the Free Software
//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
// 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

#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <string>
#include <stdexcept>
#include <array>
#include <numbers>
#include <cerrno>



#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] = pow(10, tlo + tstp*j); }

    for (i=0; i<table.imax; i++) { table.d[i] = pow(10, dlo + dstp*i); }
    
    ifstream helm_table(filename);
    if (!helm_table) {
      int errorCode = errno;
      LOG_ERROR(logger, "read_helm_table : Error ({}) opening file {}", errorCode, filename);
      throw std::runtime_error("Error (" + std::to_string(errorCode) + ") 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. 
  }
}