github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-22 20:38:01 -05:00
Code Issues Packages Projects Releases Wiki Activity

Implement all the derivatives needed for critical point determinant

Browse Source
This commit is contained in:
Ian Bell
2015-05-30 23:49:04 -06:00
parent 454759b366
commit 3fd77d57ba
3 changed files with 113 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
src/Backends/Helmholtz/MixtureDerivatives.cpp
View File

@@ -328,6 +328,31 @@ CoolPropDbl MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixt
return line1 + line2 + line3 + line4;
}
CoolPropDbl MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return Kronecker_delta(i, j)/HEOS.mole_fractions[i]+nd2nalphardnidnj__constT_V(HEOS, i, j, xN_flag);
}
CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return d_ndalphardni_dTau(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dTau__constdelta_x(HEOS, i, j, xN_flag);
}
CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return d_ndalphardni_dDelta(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dDelta__consttau_x(HEOS, i, j, xN_flag);
}
CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){
return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag) + d_nd_ndalphardni_dnj_dxk__consttau_delta(HEOS, i, j, k, xN_flag);
}
CoolPropDbl MixtureDerivatives::nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){
double sum = d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag)
+ d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag)
+ d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, k, xN_flag);
std::size_t mmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ mmax--; }
for (unsigned int m = 0; m < mmax; ++m)
{
sum -= HEOS.mole_fractions[m]*d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, m, xN_flag);
}
return sum;
}
CoolPropDbl MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
// Gernert 3.115
@@ -525,10 +550,10 @@ CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dDelta__consttau_x(Helmholt
CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag)
{
double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, i, k, xN_flag)
double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
+ d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag);
double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, i, k, xN_flag)
double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
+ d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag);
double line3 = d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, j, k, xN_flag) - d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
@@ -808,7 +833,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
SECTION(ss00.str(),"")
{
x_N_dependency_flag xN_flag;
for (std::size_t xNxN = 0; xNxN <=1; xNxN++){
for (std::size_t xNxN = 1; xNxN > 0; xNxN--){
std::ostringstream ss000;
std::string xN_string;
if (xNxN == 0){
@@ -1404,7 +1429,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
// These derivatives depend on i,j, and k indices
for (std::size_t k = 0; k < Ncomp; ++k){
if (xN_flag == XN_DEPENDENT && k == Ncomp-1){ continue; }
HelmholtzEOSMixtureBackend & rHEOS_pluszk = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_xNindep[Ncomp][k].get()) : *(HEOS_plusz_xNdep[Ncomp][k].get());
HelmholtzEOSMixtureBackend & rHEOS_minuszk = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_xNindep[Ncomp][k].get()) : *(HEOS_minusz_xNdep[Ncomp][k].get());
@@ -1414,13 +1438,13 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss1; ss1 << "d3Trdxidxjdxk, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss1.str(), "")
{
if (xN_flag == XN_INDEPENDENT){ continue; }
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d3Trdxidxjdxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d2Trdxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d2Trdxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){ err = 0; }
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-7);
@@ -1428,8 +1452,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss2; ss2 << "d3rhormolardxidxjdxk, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss2.str(), "")
{
if (xN_flag == XN_INDEPENDENT){ continue; }
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d3rhormolardxidxjdxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);
@@ -1442,8 +1465,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss3; ss3 << "d2_ndTrdni_dxj_dxk__constxi, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss3.str(), "")
{
if (xN_flag == XN_INDEPENDENT){ continue; }
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d2_ndTrdni_dxj_dxk__constxi(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);
@@ -1456,8 +1478,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss4; ss4 << "d2_ndrhorbardni_dxj_dxk__constxi, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss4.str(), "")
{
if (xN_flag == XN_INDEPENDENT){ continue; }
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d2_ndrhorbardni_dxj_dxk__constxi(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);
@@ -1470,7 +1491,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss5; ss5 << "d2_PSI_T_dxj_dxk, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss5.str(), "")
{
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d2_PSI_T_dxj_dxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);
@@ -1483,7 +1504,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
std::ostringstream ss6; ss6 << "d2_PSI_rho_dxj_dxk, i=" << i << ", j=" << j << ", k=" << k;
SECTION(ss6.str(), "")
{
if (j == Ncomp-1){ break; }
if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; }
double analytic = rHEOS.Reducing->d2_PSI_rho_dxj_dxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag);
double v1 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag);
double v2 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag);

36
src/Backends/Helmholtz/MixtureDerivatives.h
View File

@@ -149,6 +149,42 @@ class MixtureDerivatives{
static CoolPropDbl dln_fugacity_i_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static CoolPropDbl dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static CoolPropDbl dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl nAij(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
return 1/RT*ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag);
}
static CoolPropDbl n2Aijk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
return 1/RT*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS, i, j, k, xN_flag);
}
static CoolPropDbl L1_star(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){
Eigen::Matrix2d L1;
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
L1(0, 0) = nAij(HEOS, 0, 0, xN_flag);
L1(1, 0) = nAij(HEOS, 1, 0, xN_flag);
L1(0, 1) = L1(1,0);
L1(1, 1) = nAij(HEOS, 1, 1, xN_flag);
return L1.determinant();
}
static CoolPropDbl M1_star(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){
Eigen::Matrix2d M1;
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
M1(0, 0) = nAij(HEOS, 0, 0, xN_flag);
M1(0, 1) = nAij(HEOS, 0, 1, xN_flag);
double dL1_dn0 = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 0, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag);
double dL1_dn1 = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 1, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag);
M1(1, 0) = dL1_dn0;
M1(1, 1) = dL1_dn1;
return M1.determinant();
}
/** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
*

76
src/Backends/Helmholtz/ReducingFunctions.cpp
View File

@@ -337,19 +337,11 @@ CoolPropDbl GERG2008ReducingFunction::d2Yrdxidxj(const std::vector<CoolPropDbl>
CoolPropDbl d2Yr_dxidxj = 2*Yc[N-1];
d2Yr_dxidxj += c_Y_ij(i, j, beta, gamma, Y_c_ij)*d2fYijdxidxj(x,i,j,beta);
for (std::size_t k = 0; k < N-1; k++)
{
CoolPropDbl beta_Y_kN = beta[k][N-1];
d2Yr_dxidxj += 2*c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*x[k]*x[k]*(1-beta_Y_kN*beta_Y_kN)/pow(beta_Y_kN*beta_Y_kN*x[k]+xN,2)*(xN/(beta_Y_kN*beta_Y_kN*x[k]+xN)-1);
}
{
CoolPropDbl beta_Y_iN = beta[i][N-1];
d2Yr_dxidxj += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*((1-beta_Y_iN*beta_Y_iN)*(2*x[i]*xN*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 3)-x[i]*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 2)) - (x[i]+xN)/(beta_Y_iN*beta_Y_iN*x[i]+xN));
}
{
CoolPropDbl beta_Y_jN = beta[j][N-1];
d2Yr_dxidxj += -c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*((1-beta_Y_jN*beta_Y_jN)*(2*x[j]*x[j]*xN*beta_Y_jN*beta_Y_jN/pow(beta_Y_jN*beta_Y_jN*x[j]+xN, 3)-x[j]*xN/pow(beta_Y_jN*beta_Y_jN*x[j]+xN, 2)) + (x[j]+xN)/(beta_Y_jN*beta_Y_jN*x[j]+xN));
for (std::size_t k = 0; k < N-1; k++){
d2Yr_dxidxj += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*d2fYkidxi2__constxk(x, k, N-1, beta);
}
d2Yr_dxidxj -= c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, i, N-1, beta);
d2Yr_dxidxj -= c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, j, N-1, beta);
return d2Yr_dxidxj;
}
else{
@@ -387,28 +379,44 @@ CoolPropDbl GERG2008ReducingFunction::d3Yrdxidxjdxk(const std::vector<CoolPropDb
}
}
else if (xN_flag == XN_DEPENDENT){
throw ValueError("d3Yrdxidxjdxk not implemented yet");
// Table S1 from Gernert, 2014, supplemental information
if (j == N - 1 || i == N - 1){ return 0.0; }
if (i == j){ return d2Yrdxi2__constxj(x, i, beta, gamma, Y_c_ij, Yc, xN_flag); }
CoolPropDbl xN = x[N - 1];
CoolPropDbl d2Yr_dxidxj = 2 * Yc[N - 1];
d2Yr_dxidxj += c_Y_ij(i, j, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, i, j, beta);
for (std::size_t k = 0; k < N - 1; k++)
{
CoolPropDbl beta_Y_kN = beta[k][N - 1];
d2Yr_dxidxj += 2 * c_Y_ij(k, N - 1, beta, gamma, Y_c_ij)*x[k] * x[k] * (1 - beta_Y_kN*beta_Y_kN) / pow(beta_Y_kN*beta_Y_kN*x[k] + xN, 2)*(xN / (beta_Y_kN*beta_Y_kN*x[k] + xN) - 1);
}
{
CoolPropDbl beta_Y_iN = beta[i][N - 1];
d2Yr_dxidxj += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij)*((1 - beta_Y_iN*beta_Y_iN)*(2 * x[i] * xN*xN / pow(beta_Y_iN*beta_Y_iN*x[i] + xN, 3) - x[i] * xN / pow(beta_Y_iN*beta_Y_iN*x[i] + xN, 2)) - (x[i] + xN) / (beta_Y_iN*beta_Y_iN*x[i] + xN));
}
{
CoolPropDbl beta_Y_jN = beta[j][N - 1];
d2Yr_dxidxj += -c_Y_ij(j, N - 1, beta, gamma, Y_c_ij)*((1 - beta_Y_jN*beta_Y_jN)*(2 * x[j] * x[j] * xN*beta_Y_jN*beta_Y_jN / pow(beta_Y_jN*beta_Y_jN*x[j] + xN, 3) - x[j] * xN / pow(beta_Y_jN*beta_Y_jN*x[j] + xN, 2)) + (x[j] + xN) / (beta_Y_jN*beta_Y_jN*x[j] + xN));
}
return d2Yr_dxidxj;
CoolPropDbl xN = x[N - 1];
CoolPropDbl summer = 0;
// Needed for all third partials
for (std::size_t m = 0; m < N-1; m++)
{
summer -= c_Y_ij(m, N-1, beta, gamma, Y_c_ij)*d3fYkidxi3__constxk(x, m, N-1, beta);
}
double summer0 = summer;
if (i != j && j != k && k != i){
summer += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, i, N-1, beta);
summer += c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, j, N-1, beta);
summer += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, k, N-1, beta);
}
else if (k == i && i != j){ // two i, one j
summer += c_Y_ij(i, j, beta, gamma, Y_c_ij)*d3fYijdxi2dxj(x, i, j, beta);
summer += c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, j, N-1, beta);
summer += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*(2*d3fYijdxidxj2(x, i, N-1, beta) - d3fYijdxi2dxj(x, i, N-1, beta));
}
else if (k == j && i != j){ // two j, one i
summer += c_Y_ij(i, j, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, i, j, beta);
summer += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, i, N-1, beta);
summer += c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*(2*d3fYijdxidxj2(x, j, N-1, beta) - d3fYijdxi2dxj(x, j, N-1, beta));
}
else if (i == j && i != k){ // two i, one k
summer += c_Y_ij(i, k, beta, gamma, Y_c_ij)*d3fYijdxi2dxj(x, i, k, beta);
summer += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*d3fYijdxidxj2(x, k, N-1, beta);
summer += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*(2*d3fYijdxidxj2(x, i, N-1, beta) - d3fYijdxi2dxj(x, i, N-1, beta));
}
else{
for (std::size_t m = 0; m < i; m++){
summer += c_Y_ij(m, i, beta, gamma, Y_c_ij)*d3fYkidxi3__constxk(x, m, i, beta);
}
for (std::size_t m = i + 1; m < N-1; m++){
summer += c_Y_ij(i, m, beta, gamma, Y_c_ij)*d3fYikdxi3__constxk(x, i, m, beta);
}
summer += c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*(3*d3fYijdxidxj2(x, i, N-1, beta)-3*d3fYijdxi2dxj(x, i, N-1, beta)+d3fYikdxi3__constxk(x, i, N-1, beta));
}
return summer;
}
else{
throw ValueError(format("xN dependency flag invalid"));