GitHub
434964
assets/eos/helm_table.dat
Normal file
434964
assets/eos/helm_table.dat
Normal file
File diff suppressed because it is too large
Load Diff
865
src/eos/private/helmholtz.cpp
Normal file
865
src/eos/private/helmholtz.cpp
Normal file
@@ -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;
|
||||||
|
}
|
||||||
Reference in New Issue
Block a user