/* ***********************************************************************
//
//   Copyright (C) 2025 -- The 4D-STAR Collaboration
//   File Authors: Aaron Dotter, Emily Boudreaux
//   Last Modified: March 20, 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 <memory>
#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"

namespace serif::constant {
  class Constants;
}

using namespace std;

// interpolating polynomila function definitions
namespace serif::eos::helmholtz {
  double** heap_allocate_contiguous_2D_memory(const int rows, const int cols) {
    if (rows <= 0 || cols <= 0) {
      throw std::invalid_argument("Rows and columns must be positive integers.");
    }
    const auto 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;
  }
  double psi0(const double z) { return z*z*z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; }

  double dpsi0(const double z) { return z*z * ( z * (-30.0 * z + 60.0) - 30.0); }

  double ddpsi0(const double z) { return z * ( z * (-120.0 * z + 180.0) - 60.0); }

  double psi1(const double z) { return  z * ( z*z * ( z * (-3.0*z + 8.0) - 6.0) + 1.0); }

  double dpsi1(const double z) { return z*z * ( z * (-15.0*z + 32.0) - 18.0) +1.0; }

  double ddpsi1(const double z) { return z * (z * (-60.0*z + 96.0) - 36.0); }

  double psi2(const double z) { return 0.5 * z*z * ( z* ( z * (-z + 3.0) - 3.0) + 1.0); }

  double dpsi2(const double z) { return 0.5 * z * ( z * (z * (-5.0*z + 12.0) - 9.0) + 2.0); }

  double ddpsi2(const double z) { return 0.5*(z*( z * (-20.*z + 36.) - 18.) + 2.); }

  double xpsi0(const double z) { return z * z * (2.0 * z - 3.0) + 1.0; }

  double xdpsi0(const double z) { return z * (6.0*z - 6.0); }

  double xpsi1(const double z) { return z * ( z * (z - 2.0) + 1.0); }

  double xdpsi1(const double z) { return z * (3.0 * z - 4.0) + 1.0; }

  double h3(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w0mt, const double w1mt,
      const double w0d, const double w1d, const double w0md, const 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(const std::array<double, 36> &fi, const double w0t, const double w1t, const double w2t, const double w0mt,
      const double w1mt, const double w2mt, const double w0d, const double w1d, const double w2d,
      const double w0md, const double w1md, const 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
  std::unique_ptr<HELMTable> read_helm_table(const std::string &filename) {
    serif::config::Config& config = serif::config::Config::getInstance();
    auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
    serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
    quill::Logger* logger = logManager.getLogger(logFile);
    LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);

    // Make a unique pointer to the HELMTable
    auto table = std::make_unique<serif::eos::helmholtz::HELMTable>();
    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 
  ***/
  serif::eos::helmholtz::EOS get_helm_EOS(serif::eos::helmholtz::EOSInput &q, const serif::eos::helmholtz::HELMTable &table) {
    serif::config::Config& config = serif::config::Config::getInstance();
    auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
    serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
    quill::Logger* logger = logManager.getLogger(logFile);

    serif::constant::Constants& constants = serif::constant::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;
    constexpr 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);
    constexpr double qe = 4.8032042712e-10;
    constexpr double esqu = qe * qe;
    constexpr 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 = (log10(din)-dlo)*dstpi;
    iat = max(0, min(iat, table.imax-1)); // bounds

    //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

    constexpr double a1 = -0.898004;
    constexpr double b1 =  0.96786;
    constexpr double c1 =  0.220703;
    constexpr double d1 = -0.86097;
    constexpr double e1 =  2.5269;
    constexpr double a2 =  0.29561;
    constexpr double b2 =  1.9885;
    constexpr 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:
    serif::eos::helmholtz::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. 
  }
}
