34 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
36 m_IionicMolality(0.0),
37 m_maxIionicStrength(100.0),
38 m_TempPitzerRef(298.15),
39 m_IionicMolalityStoich(0.0),
40 m_form_A_Debye(A_DEBYE_WATER),
45 m_molalitiesAreCropped(false),
48 IMS_gamma_o_min_(1.0E-5),
49 IMS_gamma_k_min_(10.0),
69 CROP_ln_gamma_o_min(-6.0),
70 CROP_ln_gamma_o_max(3.0),
71 CROP_ln_gamma_k_min(-5.0),
72 CROP_ln_gamma_k_max(15.0),
75 for (
size_t i = 0; i < 17; i++) {
84 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
86 m_IionicMolality(0.0),
87 m_maxIionicStrength(100.0),
88 m_TempPitzerRef(298.15),
89 m_IionicMolalityStoich(0.0),
90 m_form_A_Debye(A_DEBYE_WATER),
95 m_molalitiesAreCropped(false),
98 IMS_gamma_o_min_(1.0E-5),
99 IMS_gamma_k_min_(10.0),
119 CROP_ln_gamma_o_min(-6.0),
120 CROP_ln_gamma_o_max(3.0),
121 CROP_ln_gamma_k_min(-5.0),
122 CROP_ln_gamma_k_max(15.0),
125 for (
int i = 0; i < 17; i++) {
135 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
137 m_IionicMolality(0.0),
138 m_maxIionicStrength(100.0),
139 m_TempPitzerRef(298.15),
140 m_IionicMolalityStoich(0.0),
141 m_form_A_Debye(A_DEBYE_WATER),
144 m_densWaterSS(1000.),
146 m_molalitiesAreCropped(false),
148 IMS_X_o_cutoff_(0.2),
149 IMS_gamma_o_min_(1.0E-5),
150 IMS_gamma_k_min_(10.0),
170 CROP_ln_gamma_o_min(-6.0),
171 CROP_ln_gamma_o_max(3.0),
172 CROP_ln_gamma_k_min(-5.0),
173 CROP_ln_gamma_k_max(15.0),
176 for (
int i = 0; i < 17; i++) {
186 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
188 m_IionicMolality(0.0),
189 m_maxIionicStrength(100.0),
190 m_TempPitzerRef(298.15),
191 m_IionicMolalityStoich(0.0),
192 m_form_A_Debye(A_DEBYE_WATER),
195 m_densWaterSS(1000.),
197 m_molalitiesAreCropped(false),
199 IMS_X_o_cutoff_(0.2),
200 IMS_gamma_o_min_(1.0E-5),
201 IMS_gamma_k_min_(10.0),
221 CROP_ln_gamma_o_min(-6.0),
222 CROP_ln_gamma_o_max(3.0),
223 CROP_ln_gamma_k_min(-5.0),
224 CROP_ln_gamma_k_max(15.0),
253 throw CanteraError(
"HMWSoln::operator=()",
"Dynamic cast to PDSS_Water failed");
367 IMS_dfCut_ = b.IMS_dfCut_;
368 IMS_efCut_ = b.IMS_efCut_;
369 IMS_afCut_ = b.IMS_afCut_;
370 IMS_bfCut_ = b.IMS_bfCut_;
372 IMS_dgCut_ = b.IMS_dgCut_;
373 IMS_egCut_ = b.IMS_egCut_;
374 IMS_agCut_ = b.IMS_agCut_;
375 IMS_bgCut_ = b.IMS_bgCut_;
379 MC_dpCut_ = b.MC_dpCut_;
380 MC_epCut_ = b.MC_epCut_;
381 MC_apCut_ = b.MC_apCut_;
382 MC_bpCut_ = b.MC_bpCut_;
383 MC_cpCut_ = b.MC_cpCut_;
384 CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
385 CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
386 CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
387 CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
398 m_formPitzerTemp(PITZER_TEMP_CONSTANT),
400 m_IionicMolality(0.0),
401 m_maxIionicStrength(100.0),
402 m_TempPitzerRef(298.15),
403 m_IionicMolalityStoich(0.0),
404 m_form_A_Debye(A_DEBYE_WATER),
407 m_densWaterSS(1000.),
409 m_molalitiesAreCropped(false),
411 IMS_X_o_cutoff_(0.2),
412 IMS_gamma_o_min_(1.0E-5),
413 IMS_gamma_k_min_(10.0),
433 CROP_ln_gamma_o_min(-6.0),
434 CROP_ln_gamma_o_max(3.0),
435 CROP_ln_gamma_k_min(-5.0),
436 CROP_ln_gamma_k_max(15.0),
440 printf(
"unknown test problem\n");
448 size_t n = i *
m_kk + j;
490 double param = -0.004;
584 for (
size_t k = 0; k <
m_kk; k++) {
596 size_t kcation =
npos;
597 double xcation = 0.0;
598 size_t kanion =
npos;
599 for (
size_t k = 0; k <
m_kk; k++) {
605 }
else if (
charge(k) < 0.0) {
606 if (
m_tmpV[k] > xcation) {
612 if (kcation ==
npos || kanion ==
npos) {
615 double xuse = xcation;
617 if (xanion < xcation) {
619 if (
charge(kcation) != 1.0) {
623 if (
charge(kanion) != 1.0) {
627 xuse = xuse / factor;
636 return hh - pres * molarV;
664 return cp - beta * beta * tt * molarV / kappa_t;
683 double* vbar = &
m_pp[0];
687 doublereal vtotal = 0.0;
688 for (
size_t i = 0; i <
m_kk; i++) {
689 vtotal += vbar[i] * x[i];
697 throw CanteraError(
"HMWSoln::isothermalCompressibility",
719 if (rho != dens_old) {
721 "Density is not an independent variable");
728 "Density is not an independent variable");
773 for (
size_t k = 1; k <
m_kk; k++) {
786 return 1.0 / mvSolvent;
792 return log(c_solvent);
797 for (
int i = 0; i < sizeUA; i++) {
802 uA[1] = -int(
nDim());
826 s_update_lnMolalityActCoeff();
830 for (
size_t k = 0; k <
m_kk; k++) {
849 s_update_lnMolalityActCoeff();
851 for (
size_t k = 0; k <
m_kk; k++) {
852 acMolality[k] = exp(acMolality[k]);
874 s_update_lnMolalityActCoeff();
877 for (
size_t k = 0; k <
m_kk; k++) {
899 for (
size_t k = 0; k <
m_kk; k++) {
906 s_update_lnMolalityActCoeff();
909 for (
size_t k = 0; k <
m_kk; k++) {
926 for (
size_t k = 0; k <
m_kk; k++) {
933 s_update_lnMolalityActCoeff();
939 for (
size_t k = 0; k <
m_kk; k++) {
955 for (
size_t k = 0; k <
m_kk; k++) {
969 s_update_lnMolalityActCoeff();
973 for (
size_t k = 0; k <
m_kk; k++) {
986 for (
size_t k = 0; k <
m_kk; k++) {
993 s_update_lnMolalityActCoeff();
999 for (
size_t k = 0; k <
m_kk; k++) {
1038 if (tempArg != -1.0) {
1042 if (presArg != -1.0) {
1055 printf(
"shouldn't be here\n");
1064 if (tempArg != -1.0) {
1068 if (presArg != -1.0) {
1081 printf(
"shouldn't be here\n");
1090 if (tempArg != -1.0) {
1094 if (presArg != -1.0) {
1106 printf(
"shouldn't be here\n");
1115 double dAphidT = dAdT /3.0;
1117 if (tempArg != -1.0) {
1126 double dAphidP = dAdP /3.0;
1128 if (tempArg != -1.0) {
1137 if (tempArg != -1.0) {
1142 double d2Aphi = d2 / 3.0;
1143 return 2.0 * A_L / T + 4.0 *
GasConstant * T * T *d2Aphi;
1149 if (tempArg != -1.0) {
1153 if (presArg != -1.0) {
1165 printf(
"shouldn't be here\n");
1186 "Unfinished func called: " + msg);
1207 size_t maxCounterIJlen = 1 + (
m_kk-1) * (
m_kk-2) / 2;
1213 int TCoeffLength = 1;
1287 m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1299 m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1308 m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1321 void HMWSoln::s_update_lnMolalityActCoeff()
const
1339 for (
size_t k = 0; k <
m_kk; k++) {
1345 double zs_k2 = z_k - zs_k1;
1369 double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1370 double lnxs = log(xx);
1372 for (
size_t k = 1; k <
m_kk; k++) {
1409 doublereal Imax = 0.0, Itmp;
1413 for (
size_t k = 0; k <
m_kk; k++) {
1424 if (cropMethod == 0) {
1435 for (
size_t i = 1; i < (m_kk - 1); i++) {
1436 double charge_i =
charge(i);
1437 double abs_charge_i = fabs(charge_i);
1438 if (charge_i == 0.0) {
1441 for (
size_t j = (i+1); j <
m_kk; j++) {
1442 double charge_j =
charge(j);
1443 double abs_charge_j = fabs(charge_j);
1452 if (charge_i * charge_j < 0) {
1457 if (Imax > Iac_max) {
1461 if (Imax > Iac_max) {
1466 if (Imax > Iac_max) {
1470 if (Imax > Iac_max) {
1482 for (
int times = 0; times< 10; times++) {
1483 double anion_charge = 0.0;
1484 double cation_charge = 0.0;
1485 size_t anion_contrib_max_i =
npos;
1486 double anion_contrib_max = -1.0;
1487 size_t cation_contrib_max_i =
npos;
1488 double cation_contrib_max = -1.0;
1489 for (
size_t i = 0; i <
m_kk; i++) {
1490 double charge_i =
charge(i);
1491 if (charge_i < 0.0) {
1493 anion_charge += anion_contrib ;
1494 if (anion_contrib > anion_contrib_max) {
1495 anion_contrib_max = anion_contrib;
1496 anion_contrib_max_i = i;
1498 }
else if (charge_i > 0.0) {
1500 cation_charge += cation_contrib ;
1501 if (cation_contrib > cation_contrib_max) {
1502 cation_contrib_max = cation_contrib;
1503 cation_contrib_max_i = i;
1507 double total_charge = cation_charge - anion_charge;
1508 if (total_charge > 1.0E-8) {
1509 double desiredCrop = total_charge/
charge(cation_contrib_max_i);
1511 if (desiredCrop < maxCrop) {
1517 }
else if (total_charge < -1.0E-8) {
1518 double desiredCrop = total_charge/
charge(anion_contrib_max_i);
1520 if (desiredCrop < maxCrop) {
1532 if (cropMethod == 1) {
1542 double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1543 double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1546 for (
size_t k = 0; k <
m_kk; k++) {
1554 for (
size_t k = 0; k <
m_kk; k++) {
1559 for (
size_t k = 0; k <
m_kk; k++) {
1574 for (i = 0; i <
m_kk; i++) {
1580 for (i = 1; i < (m_kk - 1); i++) {
1583 for (j = (i+1); j <
m_kk; j++) {
1596 size_t i, j, n, counterIJ;
1597 const double* beta0MX_coeff;
1598 const double* beta1MX_coeff;
1599 const double* beta2MX_coeff;
1600 const double* CphiMX_coeff;
1601 const double* Theta_coeff;
1604 double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1609 tquad = T * T - Tr * Tr;
1611 tinv = 1.0/T - 1.0/Tr;
1614 for (i = 1; i < (m_kk - 1); i++) {
1615 for (j = (i+1); j <
m_kk; j++) {
1630 case PITZER_TEMP_CONSTANT:
1632 case PITZER_TEMP_LINEAR:
1635 + beta0MX_coeff[1]*tlin;
1640 + beta1MX_coeff[1]*tlin;
1645 + beta2MX_coeff[1]*tlin;
1650 + CphiMX_coeff[1]*tlin;
1654 m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1660 case PITZER_TEMP_COMPLEX1:
1662 + beta0MX_coeff[1]*tlin
1663 + beta0MX_coeff[2]*tquad
1664 + beta0MX_coeff[3]*tinv
1665 + beta0MX_coeff[4]*tln;
1668 + beta1MX_coeff[1]*tlin
1669 + beta1MX_coeff[2]*tquad
1670 + beta1MX_coeff[3]*tinv
1671 + beta1MX_coeff[4]*tln;
1674 + beta2MX_coeff[1]*tlin
1675 + beta2MX_coeff[2]*tquad
1676 + beta2MX_coeff[3]*tinv
1677 + beta2MX_coeff[4]*tln;
1680 + CphiMX_coeff[1]*tlin
1681 + CphiMX_coeff[2]*tquad
1682 + CphiMX_coeff[3]*tinv
1683 + CphiMX_coeff[4]*tln;
1686 + Theta_coeff[1]*tlin
1687 + Theta_coeff[2]*tquad
1688 + Theta_coeff[3]*tinv
1689 + Theta_coeff[4]*tln;
1692 + beta0MX_coeff[2]*2.0*T
1693 - beta0MX_coeff[3]/(T*T)
1694 + beta0MX_coeff[4]/T;
1697 + beta1MX_coeff[2]*2.0*T
1698 - beta1MX_coeff[3]/(T*T)
1699 + beta1MX_coeff[4]/T;
1702 + beta2MX_coeff[2]*2.0*T
1703 - beta2MX_coeff[3]/(T*T)
1704 + beta2MX_coeff[4]/T;
1707 + CphiMX_coeff[2]*2.0*T
1708 - CphiMX_coeff[3]/(T*T)
1709 + CphiMX_coeff[4]/T;
1712 + Theta_coeff[2]*2.0*T
1713 - Theta_coeff[3]/(T*T)
1719 + beta0MX_coeff[2]*2.0
1720 + 2.0*beta0MX_coeff[3]/(T*T*T)
1721 - beta0MX_coeff[4]/(T*T);
1724 + beta1MX_coeff[2]*2.0
1725 + 2.0*beta1MX_coeff[3]/(T*T*T)
1726 - beta1MX_coeff[4]/(T*T);
1729 + beta2MX_coeff[2]*2.0
1730 + 2.0*beta2MX_coeff[3]/(T*T*T)
1731 - beta2MX_coeff[4]/(T*T);
1734 + CphiMX_coeff[2]*2.0
1735 + 2.0*CphiMX_coeff[3]/(T*T*T)
1736 - CphiMX_coeff[4]/(T*T);
1739 + Theta_coeff[2]*2.0
1740 + 2.0*Theta_coeff[3]/(T*T*T)
1741 - Theta_coeff[4]/(T*T);
1763 for (i = 1; i <
m_kk; i++) {
1765 for (j = 1; j <
m_kk; j++) {
1769 case PITZER_TEMP_CONSTANT:
1772 case PITZER_TEMP_LINEAR:
1773 m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1777 case PITZER_TEMP_COMPLEX1:
1779 + Lambda_coeff[1]*tlin
1780 + Lambda_coeff[2]*tquad
1781 + Lambda_coeff[3]*tinv
1782 + Lambda_coeff[4]*tln;
1785 + Lambda_coeff[2]*2.0*T
1786 - Lambda_coeff[3]/(T*T)
1787 + Lambda_coeff[4]/T;
1791 + 2.0*Lambda_coeff[3]/(T*T*T)
1792 - Lambda_coeff[4]/(T*T);
1798 case PITZER_TEMP_CONSTANT:
1801 case PITZER_TEMP_LINEAR:
1802 m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1806 case PITZER_TEMP_COMPLEX1:
1820 + 2.0*Mu_coeff[3]/(T*T*T)
1821 - Mu_coeff[4]/(T*T);
1829 for (i = 1; i <
m_kk; i++) {
1830 for (j = 1; j <
m_kk; j++) {
1831 for (
size_t k = 1; k <
m_kk; k++) {
1832 n = i * m_kk *m_kk + j * m_kk + k ;
1835 case PITZER_TEMP_CONSTANT:
1838 case PITZER_TEMP_LINEAR:
1839 m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1843 case PITZER_TEMP_COMPLEX1:
1846 + Psi_coeff[2]*tquad
1851 + Psi_coeff[2]*2.0*T
1852 - Psi_coeff[3]/(T*T)
1857 + 2.0*Psi_coeff[3]/(T*T*T)
1858 - Psi_coeff[4]/(T*T);
1874 printf(
"Wrong index solvent value!\n");
1884 std::string sni, snj, snk;
1911 double etheta[5][5], etheta_prime[5][5], sqrtIs;
1920 double molarcharge = 0.0;
1925 double molalitysumUncropped = 0.0;
1941 double Aphi, F, zsqF;
1942 double sum1, sum2, sum3, sum4, sum5, term1;
1943 double sum_m_phi_minus_1, osmotic_coef, lnwateract;
1946 size_t n, i, j, m, counterIJ, counterIJ2;
1950 printf(
"\n Debugging information from hmw_act \n");
1961 for (n = 1; n <
m_kk; n++) {
1965 molarcharge += fabs(
charge(n)) * molality[n];
1977 printf(
" Step 1: \n");
1978 printf(
" ionic strenth = %14.7le \n total molar "
1979 "charge = %14.7le \n", Is, molarcharge);
1997 printf(
" Step 2: \n");
2000 for (z1 = 1; z1 <=4; z1++) {
2001 for (z2 =1; z2 <=4; z2++) {
2002 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
2005 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2006 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2014 printf(
" Step 3: \n");
2015 printf(
" Species Species g(x) "
2025 for (i = 1; i < (m_kk - 1); i++) {
2026 for (j = (i+1); j <
m_kk; j++) {
2040 x1 = sqrtIs * alpha1MX[counterIJ];
2041 if (x1 > 1.0E-100) {
2042 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2043 hfunc[counterIJ] = -2.0 *
2044 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
2046 gfunc[counterIJ] = 0.0;
2047 hfunc[counterIJ] = 0.0;
2050 if (beta2MX[counterIJ] != 0.0) {
2051 x2 = sqrtIs * alpha2MX[counterIJ];
2052 if (x2 > 1.0E-100) {
2053 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2054 h2func[counterIJ] = -2.0 *
2055 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2057 g2func[counterIJ] = 0.0;
2058 h2func[counterIJ] = 0.0;
2062 gfunc[counterIJ] = 0.0;
2063 hfunc[counterIJ] = 0.0;
2069 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2070 gfunc[counterIJ], hfunc[counterIJ]);
2082 printf(
" Step 4: \n");
2083 printf(
" Species Species BMX "
2084 "BprimeMX BphiMX \n");
2088 for (i = 1; i < m_kk - 1; i++) {
2089 for (j = i+1; j <
m_kk; j++) {
2098 if (counterIJ == 2) {
2101 printf(
"beta0MX[%d] = %g\n", (
int) counterIJ, beta0MX[counterIJ]);
2102 printf(
"beta1MX[%d] = %g\n", (
int) counterIJ, beta1MX[counterIJ]);
2103 printf(
"beta2MX[%d] = %g\n", (
int) counterIJ, beta2MX[counterIJ]);
2112 BMX[counterIJ] = beta0MX[counterIJ]
2113 + beta1MX[counterIJ] * gfunc[counterIJ]
2114 + beta2MX[counterIJ] * g2func[counterIJ];
2117 printf(
"%d %g: %g %g %g %g\n",
2118 (
int) counterIJ, BMX[counterIJ], beta0MX[counterIJ],
2119 beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]);
2122 if (Is > 1.0E-150) {
2123 BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is +
2124 beta2MX[counterIJ] * h2func[counterIJ]/Is);
2126 BprimeMX[counterIJ] = 0.0;
2128 BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ];
2130 BMX[counterIJ] = 0.0;
2131 BprimeMX[counterIJ] = 0.0;
2132 BphiMX[counterIJ] = 0.0;
2138 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
2139 sni.c_str(), snj.c_str(),
2140 BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ]);
2152 printf(
" Step 5: \n");
2153 printf(
" Species Species CMX \n");
2156 for (i = 1; i < m_kk-1; i++) {
2157 for (j = i+1; j <
m_kk; j++) {
2168 CMX[counterIJ] = CphiMX[counterIJ]/
2171 CMX[counterIJ] = 0.0;
2175 if (counterIJ == 2) {
2178 printf(
"CphiMX[%d] = %g\n", (
int) counterIJ, CphiMX[counterIJ]);
2186 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
2199 printf(
" Step 6: \n");
2200 printf(
" Species Species Phi_ij "
2201 " Phiprime_ij Phi^phi_ij \n");
2204 for (i = 1; i < m_kk-1; i++) {
2205 for (j = i+1; j <
m_kk; j++) {
2216 z1 = (int) fabs(
charge(i));
2217 z2 = (int) fabs(
charge(j));
2218 Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
2219 Phiprime[counterIJ] = etheta_prime[z1][z2];
2220 Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
2222 Phi[counterIJ] = 0.0;
2223 Phiprime[counterIJ] = 0.0;
2224 Phiphi[counterIJ] = 0.0;
2230 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
2231 sni.c_str(), snj.c_str(),
2232 Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ]);
2244 printf(
" Step 7: \n");
2252 F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2253 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2256 printf(
"Aphi = %20.13g\n", Aphi);
2261 printf(
" initial value of F = %10.6f \n", F);
2264 for (i = 1; i < m_kk-1; i++) {
2265 for (j = i+1; j <
m_kk; j++) {
2276 F = F + molality[i]*molality[j] * BprimeMX[counterIJ];
2283 F = F + molality[i]*molality[j] * Phiprime[counterIJ];
2287 printf(
" F = %10.6f \n", F);
2294 printf(
" Step 8: Summing in All Contributions to Activity Coefficients \n");
2298 for (i = 1; i <
m_kk; i++) {
2310 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2317 printf(
" Unary term: z*z*F = %10.5f\n", zsqF);
2325 for (j = 1; j <
m_kk; j++) {
2334 sum1 = sum1 + molality[j]*
2335 (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]);
2339 printf(
" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2340 molality[j]*2.0*BMX[counterIJ]);
2341 printf(
" m_j Z CMX = %10.5f\n",
2342 molality[j]* molarcharge*CMX[counterIJ]);
2351 for (
size_t k = j+1; k <
m_kk; k++) {
2354 n = k + j * m_kk + i * m_kk *
m_kk;
2355 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2358 if (psi_ijk[n] != 0.0) {
2360 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2361 molality[j]*molality[k]*psi_ijk[n]);
2374 sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2377 if ((molality[j] * Phi[counterIJ])!= 0.0) {
2379 printf(
" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(),
2380 molality[j]*(2.0*Phi[counterIJ]));
2385 for (
size_t k = 1; k <
m_kk; k++) {
2389 n = k + j * m_kk + i * m_kk *
m_kk;
2390 sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2393 if (psi_ijk[n] != 0.0) {
2395 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2396 molality[j]*molality[k]*psi_ijk[n]);
2405 sum4 = sum4 + (fabs(
charge(i))*
2406 molality[j]*molality[k]*CMX[counterIJ2]);
2409 if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2411 printf(
" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2412 fabs(
charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2429 printf(
" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2437 for (
size_t k = 1; k <
m_kk; k++) {
2441 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2442 double zeta = psi_ijk[n];
2444 sum5 = sum5 + molality[j]*molality[k]*zeta;
2448 printf(
" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(),
2449 molality[j]*molality[k]*psi_ijk[n]);
2466 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
2482 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2490 printf(
" Unary term: z*z*F = %10.5f\n", zsqF);
2498 for (j = 1; j <
m_kk; j++) {
2509 sum1 = sum1 + molality[j]*
2510 (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]);
2514 printf(
" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2515 molality[j]*2.0*BMX[counterIJ]);
2516 printf(
" m_j Z CMX = %10.5f\n",
2517 molality[j]* molarcharge*CMX[counterIJ]);
2521 for (
size_t k = j+1; k <
m_kk; k++) {
2524 n = k + j * m_kk + i * m_kk *
m_kk;
2525 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2528 if (psi_ijk[n] != 0.0) {
2530 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2531 molality[j]*molality[k]*psi_ijk[n]);
2546 sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2549 if ((molality[j] * Phi[counterIJ])!= 0.0) {
2551 printf(
" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(),
2552 molality[j]*(2.0*Phi[counterIJ]));
2557 for (
size_t k = 1; k <
m_kk; k++) {
2560 n = k + j * m_kk + i * m_kk *
m_kk;
2561 sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2564 if (psi_ijk[n] != 0.0) {
2566 printf(
" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2567 molality[j]*molality[k]*psi_ijk[n]);
2578 molality[j]*molality[k]*CMX[counterIJ2]);
2581 if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2583 printf(
" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2584 fabs(
charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2601 printf(
" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2609 for (
size_t k = 1; k <
m_kk; k++) {
2614 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2615 double zeta = psi_ijk[n];
2617 sum5 = sum5 + molality[j]*molality[k]*zeta;
2621 printf(
" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(),
2622 molality[j]*molality[k]*psi_ijk[n]);
2635 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2649 printf(
" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2654 for (j = 1; j <
m_kk; j++) {
2660 printf(
" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(),
2670 for (
size_t k = 1; k <
m_kk; k++) {
2672 n = k + j * m_kk + i * m_kk *
m_kk;
2673 sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2676 if (psi_ijk[n] != 0.0) {
2678 printf(
" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2679 molality[j]*molality[k]*psi_ijk[n]);
2687 sum2 = 3.0 * molality[i]* molality[i] *
m_Mu_nnn[i];
2690 if (m_Mu_nnn[i] != 0.0) {
2691 printf(
" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
2692 3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
2702 printf(
" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2711 printf(
" Step 9: \n");
2733 term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2735 for (j = 1; j <
m_kk; j++) {
2740 for (
size_t k = 1; k <
m_kk; k++) {
2748 sum1 = sum1 + molality[j]*molality[k]*
2749 (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]);
2753 for (
size_t k = j+1; k <
m_kk; k++) {
2754 if (j == (m_kk-1)) {
2756 printf(
"logic error 1 in Step 9 of hmw_act");
2766 sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ];
2767 for (m = 1; m <
m_kk; m++) {
2770 n = m + k * m_kk + j * m_kk *
m_kk;
2772 molality[j]*molality[k]*molality[m]*psi_ijk[n];
2783 for (
size_t k = j+1; k <
m_kk; k++) {
2786 printf(
"logic error 2 in Step 9 of hmw_act");
2797 sum3 = sum3 + molality[j]*molality[k]*Phiphi[counterIJ];
2798 for (m = 1; m <
m_kk; m++) {
2800 n = m + k * m_kk + j * m_kk *
m_kk;
2802 molality[j]*molality[k]*molality[m]*psi_ijk[n];
2813 for (
size_t k = 1; k <
m_kk; k++) {
2815 sum4 = sum4 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
2818 sum5 = sum5 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
2822 sum6 = sum6 + molality[j]*molality[k]*
m_Lambda_nj(j,k);
2823 }
else if (k == j) {
2824 sum6 = sum6 + 0.5 * molality[j]*molality[k]*
m_Lambda_nj(j,k);
2829 for (m = 1; m <
m_kk; m++) {
2832 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
2833 double zeta = psi_ijk[n];
2835 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2841 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn[j];
2844 sum_m_phi_minus_1 = 2.0 *
2845 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2850 if (molalitysumUncropped > 1.0E-150) {
2851 osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2857 printf(
"OsmCoef - 1 = %20.13g\n", osmotic_coef - 1.0);
2862 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
2863 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
2864 term1, sum1, sum2, sum3, sum4, sum5);
2865 printf(
" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n",
2866 sum_m_phi_minus_1, osmotic_coef);
2870 printf(
" Step 10: \n");
2873 lnwateract = -(
m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2888 double wateract = exp(lnwateract);
2890 printf(
" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
2891 printf(
" ln_a_water=%10.6f a_water=%10.6f\n\n",
2892 lnwateract, wateract);
2912 for (
size_t k = 1; k <
m_kk; k++) {
2946 printf(
"Wrong index solvent value!\n");
2950 std::string sni, snj, snk;
2965 double etheta[5][5], etheta_prime[5][5], sqrtIs;
2974 double molarcharge = 0.0;
2979 double molalitysum = 0.0;
2994 double dFdT, zsqdFdT;
2995 double sum1, sum2, sum3, sum4, sum5, term1;
2996 double sum_m_phi_minus_1, d_osmotic_coef_dT, d_lnwateract_dT;
2999 size_t n, i, j, m, counterIJ, counterIJ2;
3003 printf(
"\n Debugging information from "
3004 "s_Pitzer_dlnMolalityActCoeff_dT()\n");
3015 for (n = 1; n <
m_kk; n++) {
3019 molarcharge += fabs(
charge(n)) * molality[n];
3020 molalitysum += molality[n];
3031 printf(
" Step 1: \n");
3032 printf(
" ionic strenth = %14.7le \n total molar "
3033 "charge = %14.7le \n", Is, molarcharge);
3051 printf(
" Step 2: \n");
3054 for (z1 = 1; z1 <=4; z1++) {
3055 for (z2 =1; z2 <=4; z2++) {
3056 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
3059 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3060 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3068 printf(
" Step 3: \n");
3069 printf(
" Species Species g(x) "
3079 for (i = 1; i < (m_kk - 1); i++) {
3080 for (j = (i+1); j <
m_kk; j++) {
3093 x1 = sqrtIs * alpha1MX[counterIJ];
3094 if (x1 > 1.0E-100) {
3095 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3096 hfunc[counterIJ] = -2.0 *
3097 (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
3099 gfunc[counterIJ] = 0.0;
3100 hfunc[counterIJ] = 0.0;
3103 if (beta2MX_L[counterIJ] != 0.0) {
3104 x2 = sqrtIs * alpha2MX[counterIJ];
3105 if (x2 > 1.0E-100) {
3106 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3107 h2func[counterIJ] = -2.0 *
3108 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3110 g2func[counterIJ] = 0.0;
3111 h2func[counterIJ] = 0.0;
3115 gfunc[counterIJ] = 0.0;
3116 hfunc[counterIJ] = 0.0;
3122 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3123 gfunc[counterIJ], hfunc[counterIJ]);
3136 printf(
" Step 4: \n");
3137 printf(
" Species Species BMX "
3138 "BprimeMX BphiMX \n");
3142 for (i = 1; i < m_kk - 1; i++) {
3143 for (j = i+1; j <
m_kk; j++) {
3154 BMX_L[counterIJ] = beta0MX_L[counterIJ]
3155 + beta1MX_L[counterIJ] * gfunc[counterIJ]
3156 + beta2MX_L[counterIJ] * gfunc[counterIJ];
3159 printf(
"%d %g: %g %g %g %g\n",
3160 (
int) counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ],
3161 beta1MX_L[counterIJ], beta2MX_L[counterIJ], gfunc[counterIJ]);
3164 if (Is > 1.0E-150) {
3165 BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is +
3166 beta2MX_L[counterIJ] * h2func[counterIJ]/Is);
3168 BprimeMX_L[counterIJ] = 0.0;
3170 BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ];
3172 BMX_L[counterIJ] = 0.0;
3173 BprimeMX_L[counterIJ] = 0.0;
3174 BphiMX_L[counterIJ] = 0.0;
3180 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
3181 sni.c_str(), snj.c_str(),
3182 BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]);
3194 printf(
" Step 5: \n");
3195 printf(
" Species Species CMX \n");
3198 for (i = 1; i < m_kk-1; i++) {
3199 for (j = i+1; j <
m_kk; j++) {
3210 CMX_L[counterIJ] = CphiMX_L[counterIJ]/
3213 CMX_L[counterIJ] = 0.0;
3219 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
3232 printf(
" Step 6: \n");
3233 printf(
" Species Species Phi_ij "
3234 " Phiprime_ij Phi^phi_ij \n");
3237 for (i = 1; i < m_kk-1; i++) {
3238 for (j = i+1; j <
m_kk; j++) {
3249 z1 = (int) fabs(
charge(i));
3250 z2 = (int) fabs(
charge(j));
3252 Phi_L[counterIJ] = thetaij_L[counterIJ];
3254 Phiprime[counterIJ] = 0.0;
3255 Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
3257 Phi_L[counterIJ] = 0.0;
3258 Phiprime[counterIJ] = 0.0;
3259 Phiphi_L[counterIJ] = 0.0;
3265 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
3266 sni.c_str(), snj.c_str(),
3267 Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ]);
3278 printf(
" Step 7: \n");
3287 double dAphidT = dA_DebyedT /3.0;
3294 dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3295 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3298 printf(
" initial value of dFdT = %10.6f \n", dFdT);
3301 for (i = 1; i < m_kk-1; i++) {
3302 for (j = i+1; j <
m_kk; j++) {
3313 dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ];
3320 dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ];
3324 printf(
" dFdT = %10.6f \n", dFdT);
3331 printf(
" Step 8: \n");
3335 for (i = 1; i <
m_kk; i++) {
3349 for (j = 1; j <
m_kk; j++) {
3358 sum1 = sum1 + molality[j]*
3359 (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3366 for (
size_t k = j+1; k <
m_kk; k++) {
3369 n = k + j * m_kk + i * m_kk *
m_kk;
3370 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3380 sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3382 for (
size_t k = 1; k <
m_kk; k++) {
3386 n = k + j * m_kk + i * m_kk *
m_kk;
3387 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3393 sum4 = sum4 + (fabs(
charge(i))*
3394 molality[j]*molality[k]*CMX_L[counterIJ2]);
3408 for (
size_t k = 1; k <
m_kk; k++) {
3412 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3413 double zeta_L = psi_ijk_L[n];
3414 if (zeta_L != 0.0) {
3415 sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3425 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3430 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3432 printf(
" %12g %12g %12g %12g %12g %12g\n",
3433 zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3450 for (j = 1; j <
m_kk; j++) {
3461 sum1 = sum1 + molality[j]*
3462 (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3464 for (
size_t k = j+1; k <
m_kk; k++) {
3467 n = k + j * m_kk + i * m_kk *
m_kk;
3468 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3480 sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3482 for (
size_t k = 1; k <
m_kk; k++) {
3485 n = k + j * m_kk + i * m_kk *
m_kk;
3486 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3494 molality[j]*molality[k]*CMX_L[counterIJ2]);
3504 for (
size_t k = 1; k <
m_kk; k++) {
3509 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3510 double zeta_L = psi_ijk_L[n];
3511 if (zeta_L != 0.0) {
3512 sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3519 zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3524 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
3526 printf(
" %12g %12g %12g %12g %12g %12g\n",
3527 zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3539 for (j = 1; j <
m_kk; j++) {
3545 for (
size_t k = 1; k <
m_kk; k++) {
3547 n = k + j * m_kk + i * m_kk *
m_kk;
3548 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3553 sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_L[i];
3559 printf(
" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3568 printf(
" Step 9: \n");
3590 term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3592 for (j = 1; j <
m_kk; j++) {
3597 for (
size_t k = 1; k <
m_kk; k++) {
3605 sum1 = sum1 + molality[j]*molality[k]*
3606 (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3610 for (
size_t k = j+1; k <
m_kk; k++) {
3611 if (j == (m_kk-1)) {
3613 printf(
"logic error 1 in Step 9 of hmw_act");
3623 sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3624 for (m = 1; m <
m_kk; m++) {
3627 n = m + k * m_kk + j * m_kk *
m_kk;
3629 molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3640 for (
size_t k = j+1; k <
m_kk; k++) {
3643 printf(
"logic error 2 in Step 9 of hmw_act");
3654 sum3 = sum3 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3655 for (m = 1; m <
m_kk; m++) {
3657 n = m + k * m_kk + j * m_kk *
m_kk;
3659 molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3670 for (
size_t k = 1; k <
m_kk; k++) {
3680 }
else if (k == j) {
3681 sum6 = sum6 + 0.5 * molality[j]*molality[k]*
m_Lambda_nj_L(j,k);
3686 for (m = 1; m <
m_kk; m++) {
3689 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
3690 double zeta_L = psi_ijk_L[n];
3691 if (zeta_L != 0.0) {
3692 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
3698 sum7 += molality[j]*molality[j]*molality[j]*
m_Mu_nnn_L[j];
3701 sum_m_phi_minus_1 = 2.0 *
3702 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3707 if (molalitysum > 1.0E-150) {
3708 d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3710 d_osmotic_coef_dT = 0.0;
3714 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3715 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3716 term1, sum1, sum2, sum3, sum4, sum5);
3717 printf(
" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n",
3718 sum_m_phi_minus_1, d_osmotic_coef_dT);
3722 printf(
" Step 10: \n");
3725 d_lnwateract_dT = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
3739 double d_wateract_dT = exp(d_lnwateract_dT);
3740 printf(
" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
3741 d_lnwateract_dT, d_wateract_dT);
3761 for (
size_t k = 1; k <
m_kk; k++) {
3786 printf(
"Wrong index solvent value!\n");
3790 std::string sni, snj, snk;
3805 double etheta[5][5], etheta_prime[5][5], sqrtIs;
3814 double molarcharge = 0.0;
3819 double molalitysum = 0.0;
3835 double d2FdT2, zsqd2FdT2;
3836 double sum1, sum2, sum3, sum4, sum5, term1;
3837 double sum_m_phi_minus_1, d2_osmotic_coef_dT2, d2_lnwateract_dT2;
3840 size_t n, i, j, m, counterIJ, counterIJ2;
3844 printf(
"\n Debugging information from "
3845 "s_Pitzer_d2lnMolalityActCoeff_dT2()\n");
3857 for (n = 1; n <
m_kk; n++) {
3861 molarcharge += fabs(
charge(n)) * molality[n];
3862 molalitysum += molality[n];
3873 printf(
" Step 1: \n");
3874 printf(
" ionic strenth = %14.7le \n total molar "
3875 "charge = %14.7le \n", Is, molarcharge);
3893 printf(
" Step 2: \n");
3896 for (z1 = 1; z1 <=4; z1++) {
3897 for (z2 =1; z2 <=4; z2++) {
3898 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
3901 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3902 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3910 printf(
" Step 3: \n");
3911 printf(
" Species Species g(x) "
3922 for (i = 1; i < (m_kk - 1); i++) {
3923 for (j = (i+1); j <
m_kk; j++) {
3936 x1 = sqrtIs * alpha1MX[counterIJ];
3937 if (x1 > 1.0E-100) {
3938 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
3939 hfunc[counterIJ] = -2.0*
3940 (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
3942 gfunc[counterIJ] = 0.0;
3943 hfunc[counterIJ] = 0.0;
3946 if (beta2MX_LL[counterIJ] != 0.0) {
3947 x2 = sqrtIs * alpha2MX[counterIJ];
3948 if (x2 > 1.0E-100) {
3949 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3950 h2func[counterIJ] = -2.0 *
3951 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3953 g2func[counterIJ] = 0.0;
3954 h2func[counterIJ] = 0.0;
3958 gfunc[counterIJ] = 0.0;
3959 hfunc[counterIJ] = 0.0;
3965 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3966 gfunc[counterIJ], hfunc[counterIJ]);
3978 printf(
" Step 4: \n");
3979 printf(
" Species Species BMX "
3980 "BprimeMX BphiMX \n");
3984 for (i = 1; i < m_kk - 1; i++) {
3985 for (j = i+1; j <
m_kk; j++) {
3996 BMX_LL[counterIJ] = beta0MX_LL[counterIJ]
3997 + beta1MX_LL[counterIJ] * gfunc[counterIJ]
3998 + beta2MX_LL[counterIJ] * g2func[counterIJ];
4001 printf(
"%d %g: %g %g %g %g\n",
4002 (
int) counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ],
4003 beta1MX_LL[counterIJ], beta2MX_LL[counterIJ], gfunc[counterIJ]);
4006 if (Is > 1.0E-150) {
4007 BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is +
4008 beta2MX_LL[counterIJ] * h2func[counterIJ]/Is);
4010 BprimeMX_LL[counterIJ] = 0.0;
4012 BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ];
4014 BMX_LL[counterIJ] = 0.0;
4015 BprimeMX_LL[counterIJ] = 0.0;
4016 BphiMX_LL[counterIJ] = 0.0;
4022 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
4023 sni.c_str(), snj.c_str(),
4024 BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]);
4036 printf(
" Step 5: \n");
4037 printf(
" Species Species CMX \n");
4040 for (i = 1; i < m_kk-1; i++) {
4041 for (j = i+1; j <
m_kk; j++) {
4052 CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/
4055 CMX_LL[counterIJ] = 0.0;
4061 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4074 printf(
" Step 6: \n");
4075 printf(
" Species Species Phi_ij "
4076 " Phiprime_ij Phi^phi_ij \n");
4079 for (i = 1; i < m_kk-1; i++) {
4080 for (j = i+1; j <
m_kk; j++) {
4091 z1 = (int) fabs(
charge(i));
4092 z2 = (int) fabs(
charge(j));
4095 Phi_LL[counterIJ] = thetaij_LL[counterIJ];
4097 Phiprime[counterIJ] = 0.0;
4100 Phiphi_LL[counterIJ] = Phi_LL[counterIJ];
4102 Phi_LL[counterIJ] = 0.0;
4103 Phiprime[counterIJ] = 0.0;
4104 Phiphi_LL[counterIJ] = 0.0;
4113 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
4114 sni.c_str(), snj.c_str(),
4115 Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]);
4126 printf(
" Step 7: \n");
4146 d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4147 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4150 printf(
" initial value of d2FdT2 = %10.6f \n", d2FdT2);
4153 for (i = 1; i < m_kk-1; i++) {
4154 for (j = i+1; j <
m_kk; j++) {
4165 d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ];
4172 d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ];
4176 printf(
" d2FdT2 = %10.6f \n", d2FdT2);
4183 printf(
" Step 8: \n");
4187 for (i = 1; i <
m_kk; i++) {
4201 for (j = 1; j <
m_kk; j++) {
4210 sum1 = sum1 + molality[j]*
4211 (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4218 for (
size_t k = j+1; k <
m_kk; k++) {
4221 n = k + j * m_kk + i * m_kk *
m_kk;
4222 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4232 sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4234 for (
size_t k = 1; k <
m_kk; k++) {
4238 n = k + j * m_kk + i * m_kk *
m_kk;
4239 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4245 sum4 = sum4 + (fabs(
charge(i))*
4246 molality[j]*molality[k]*CMX_LL[counterIJ2]);
4259 for (
size_t k = 1; k <
m_kk; k++) {
4263 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
4264 double zeta_LL = psi_ijk_LL[n];
4265 if (zeta_LL != 0.0) {
4266 sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4277 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4281 printf(
" %-16s d2lngammadT2[i]=%10.6f \n",
4283 printf(
" %12g %12g %12g %12g %12g %12g\n",
4284 zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4302 for (j = 1; j <
m_kk; j++) {
4313 sum1 = sum1 + molality[j]*
4314 (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4316 for (
size_t k = j+1; k <
m_kk; k++) {
4319 n = k + j * m_kk + i * m_kk *
m_kk;
4320 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4332 sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
4334 for (
size_t k = 1; k <
m_kk; k++) {
4337 n = k + j * m_kk + i * m_kk *
m_kk;
4338 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
4346 molality[j]*molality[k]*CMX_LL[counterIJ2]);
4359 for (
size_t k = 1; k <
m_kk; k++) {
4364 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
4365 double zeta_LL = psi_ijk_LL[n];
4366 if (zeta_LL != 0.0) {
4367 sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
4374 zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
4378 printf(
" %-16s d2lngammadT2[i]=%10.6f\n",
4380 printf(
" %12g %12g %12g %12g %12g %12g\n",
4381 zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
4393 for (j = 1; j <
m_kk; j++) {
4399 for (
size_t k = 1; k <
m_kk; k++) {
4401 n = k + j * m_kk + i * m_kk *
m_kk;
4402 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
4407 sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_LL[i];
4412 printf(
" %-16s d2lngammadT2[i]=%10.6f \n",
4421 printf(
" Step 9: \n");
4443 term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
4445 for (j = 1; j <
m_kk; j++) {
4450 for (
size_t k = 1; k <
m_kk; k++) {
4458 sum1 = sum1 + molality[j]*molality[k]*
4459 (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
4463 for (
size_t k = j+1; k <
m_kk; k++) {
4464 if (j == (m_kk-1)) {
4466 printf(
"logic error 1 in Step 9 of hmw_act");
4476 sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
4477 for (m = 1; m <
m_kk; m++) {
4480 n = m + k * m_kk + j * m_kk *
m_kk;
4482 molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
4493 for (
size_t k = j+1; k <
m_kk; k++) {
4496 printf(
"logic error 2 in Step 9 of hmw_act");
4507 sum3 = sum3 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
4508 for (m = 1; m <
m_kk; m++) {
4510 n = m + k * m_kk + j * m_kk *
m_kk;
4512 molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
4523 for (
size_t k = 1; k <
m_kk; k++) {
4533 }
else if (k == j) {
4539 for (m = 1; m <
m_kk; m++) {
4542 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
4543 double zeta_LL = psi_ijk_LL[n];
4544 if (zeta_LL != 0.0) {
4545 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
4552 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_LL[j];
4555 sum_m_phi_minus_1 = 2.0 *
4556 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4561 if (molalitysum > 1.0E-150) {
4562 d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4564 d2_osmotic_coef_dT2 = 0.0;
4568 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
4569 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
4570 term1, sum1, sum2, sum3, sum4, sum5);
4571 printf(
" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n",
4572 sum_m_phi_minus_1, d2_osmotic_coef_dT2);
4576 printf(
" Step 10: \n");
4579 d2_lnwateract_dT2 = -(
m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
4593 double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
4594 printf(
" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
4595 d2_lnwateract_dT2, d2_wateract_dT2);
4605 for (
size_t k = 1; k <
m_kk; k++) {
4627 printf(
"Wrong index solvent value!\n");
4631 std::string sni, snj, snk;
4646 double etheta[5][5], etheta_prime[5][5], sqrtIs;
4655 double molarcharge = 0.0;
4660 double molalitysum = 0.0;
4675 double dFdP, zsqdFdP;
4676 double sum1, sum2, sum3, sum4, sum5, term1;
4677 double sum_m_phi_minus_1, d_osmotic_coef_dP, d_lnwateract_dP;
4680 size_t n, i, j, m, counterIJ, counterIJ2;
4687 printf(
"\n Debugging information from "
4688 "s_Pitzer_dlnMolalityActCoeff_dP()\n");
4699 for (n = 1; n <
m_kk; n++) {
4703 molarcharge += fabs(
charge(n)) * molality[n];
4704 molalitysum += molality[n];
4715 printf(
" Step 1: \n");
4716 printf(
" ionic strenth = %14.7le \n total molar "
4717 "charge = %14.7le \n", Is, molarcharge);
4736 printf(
" Step 2: \n");
4739 for (z1 = 1; z1 <=4; z1++) {
4740 for (z2 =1; z2 <=4; z2++) {
4741 calc_thetas(z1, z2, ðeta[z1][z2], ðeta_prime[z1][z2]);
4744 printf(
" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
4745 z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
4753 printf(
" Step 3: \n");
4754 printf(
" Species Species g(x) "
4765 for (i = 1; i < (m_kk - 1); i++) {
4766 for (j = (i+1); j <
m_kk; j++) {
4779 x1 = sqrtIs * alpha1MX[counterIJ];
4780 if (x1 > 1.0E-100) {
4781 gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
4782 hfunc[counterIJ] = -2.0*
4783 (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
4785 gfunc[counterIJ] = 0.0;
4786 hfunc[counterIJ] = 0.0;
4789 if (beta2MX_P[counterIJ] != 0.0) {
4790 x2 = sqrtIs * alpha2MX[counterIJ];
4791 if (x2 > 1.0E-100) {
4792 g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
4793 h2func[counterIJ] = -2.0 *
4794 (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
4796 g2func[counterIJ] = 0.0;
4797 h2func[counterIJ] = 0.0;
4801 gfunc[counterIJ] = 0.0;
4802 hfunc[counterIJ] = 0.0;
4808 printf(
" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
4809 gfunc[counterIJ], hfunc[counterIJ]);
4822 printf(
" Step 4: \n");
4823 printf(
" Species Species BMX "
4824 "BprimeMX BphiMX \n");
4828 for (i = 1; i < m_kk - 1; i++) {
4829 for (j = i+1; j <
m_kk; j++) {
4840 BMX_P[counterIJ] = beta0MX_P[counterIJ]
4841 + beta1MX_P[counterIJ] * gfunc[counterIJ]
4842 + beta2MX_P[counterIJ] * g2func[counterIJ];
4845 printf(
"%d %g: %g %g %g %g\n",
4846 (
int) counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ],
4847 beta1MX_P[counterIJ], beta2MX_P[counterIJ], gfunc[counterIJ]);
4850 if (Is > 1.0E-150) {
4851 BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is +
4852 beta2MX_P[counterIJ] * h2func[counterIJ]/Is);
4854 BprimeMX_P[counterIJ] = 0.0;
4856 BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ];
4858 BMX_P[counterIJ] = 0.0;
4859 BprimeMX_P[counterIJ] = 0.0;
4860 BphiMX_P[counterIJ] = 0.0;
4866 printf(
" %-16s %-16s %11.7f %11.7f %11.7f \n",
4867 sni.c_str(), snj.c_str(),
4868 BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]);
4879 printf(
" Step 5: \n");
4880 printf(
" Species Species CMX \n");
4883 for (i = 1; i < m_kk-1; i++) {
4884 for (j = i+1; j <
m_kk; j++) {
4895 CMX_P[counterIJ] = CphiMX_P[counterIJ]/
4898 CMX_P[counterIJ] = 0.0;
4904 printf(
" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4917 printf(
" Step 6: \n");
4918 printf(
" Species Species Phi_ij "
4919 " Phiprime_ij Phi^phi_ij \n");
4922 for (i = 1; i < m_kk-1; i++) {
4923 for (j = i+1; j <
m_kk; j++) {
4934 z1 = (int) fabs(
charge(i));
4935 z2 = (int) fabs(
charge(j));
4937 Phi_P[counterIJ] = thetaij_P[counterIJ];
4939 Phiprime[counterIJ] = 0.0;
4940 Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ];
4942 Phi_P[counterIJ] = 0.0;
4943 Phiprime[counterIJ] = 0.0;
4944 Phiphi_P[counterIJ] = 0.0;
4950 printf(
" %-16s %-16s %10.6f %10.6f %10.6f \n",
4951 sni.c_str(), snj.c_str(),
4952 Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ]);
4963 printf(
" Step 7: \n");
4972 double dAphidP = dA_DebyedP /3.0;
4979 dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4980 + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4983 printf(
" initial value of dFdP = %10.6f \n", dFdP);
4986 for (i = 1; i < m_kk-1; i++) {
4987 for (j = i+1; j <
m_kk; j++) {
4998 dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ];
5005 dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ];
5009 printf(
" dFdP = %10.6f \n", dFdP);
5016 printf(
" Step 8: \n");
5021 for (i = 1; i <
m_kk; i++) {
5034 for (j = 1; j <
m_kk; j++) {
5043 sum1 = sum1 + molality[j]*
5044 (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5051 for (
size_t k = j+1; k <
m_kk; k++) {
5054 n = k + j * m_kk + i * m_kk *
m_kk;
5055 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5065 sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5067 for (
size_t k = 1; k <
m_kk; k++) {
5071 n = k + j * m_kk + i * m_kk *
m_kk;
5072 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5078 sum4 = sum4 + (fabs(
charge(i))*
5079 molality[j]*molality[k]*CMX_P[counterIJ2]);
5092 for (
size_t k = 1; k <
m_kk; k++) {
5096 n = izeta * m_kk * m_kk + jzeta * m_kk + k;
5097 double zeta_P = psi_ijk_P[n];
5098 if (zeta_P != 0.0) {
5099 sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5111 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5116 printf(
" %-16s lngamma[i]=%10.6f \n",
5118 printf(
" %12g %12g %12g %12g %12g %12g\n",
5119 zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5136 for (j = 1; j <
m_kk; j++) {
5147 sum1 = sum1 + molality[j]*
5148 (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5150 for (
size_t k = j+1; k <
m_kk; k++) {
5153 n = k + j * m_kk + i * m_kk *
m_kk;
5154 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5166 sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
5168 for (
size_t k = 1; k <
m_kk; k++) {
5171 n = k + j * m_kk + i * m_kk *
m_kk;
5172 sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
5180 molality[j]*molality[k]*CMX_P[counterIJ2]);
5193 for (
size_t k = 1; k <
m_kk; k++) {
5198 n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
5199 double zeta_P = psi_ijk_P[n];
5200 if (zeta_P != 0.0) {
5201 sum5 = sum5 + molality[j]*molality[k]*zeta_P;
5208 zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
5212 printf(
" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
5214 printf(
" %12g %12g %12g %12g %12g %12g\n",
5215 zsqdFdP, sum1, sum2, sum3, sum4, sum5);
5226 for (j = 1; j <
m_kk; j++) {
5232 for (
size_t k = 1; k <
m_kk; k++) {
5234 n = k + j * m_kk + i * m_kk *
m_kk;
5235 sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
5240 sum2 = 3.0 * molality[i] * molality[i] *
m_Mu_nnn_P[i];
5245 printf(
" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
5254 printf(
" Step 9: \n");
5276 term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
5278 for (j = 1; j <
m_kk; j++) {
5283 for (
size_t k = 1; k <
m_kk; k++) {
5291 sum1 = sum1 + molality[j]*molality[k]*
5292 (BphiMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
5296 for (
size_t k = j+1; k <
m_kk; k++) {
5297 if (j == (m_kk-1)) {
5299 printf(
"logic error 1 in Step 9 of hmw_act");
5309 sum2 = sum2 + molality[j]*molality[k]*Phiphi_P[counterIJ];
5310 for (m = 1; m <
m_kk; m++) {
5313 n = m + k * m_kk + j * m_kk *
m_kk;
5315 molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
5327 for (
size_t k = j+1; k <
m_kk; k++) {
5330 printf(
"logic error 2 in Step 9 of hmw_act");
5341 sum3 = sum3 + molality[j]*molality[k]*Phiphi_P[counterIJ];
5342 for (m = 1; m <
m_kk; m++) {
5344 n = m + k * m_kk + j * m_kk *
m_kk;
5346 molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
5357 for (
size_t k = 1; k <
m_kk; k++) {
5367 }
else if (k == j) {
5368 sum6 = sum6 + 0.5 * molality[j]*molality[k]*
m_Lambda_nj_P(j,k);
5373 for (m = 1; m <
m_kk; m++) {
5376 n = k + jzeta * m_kk + izeta * m_kk *
m_kk;
5377 double zeta_P = psi_ijk_P[n];
5378 if (zeta_P != 0.0) {
5379 sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
5386 sum7 += molality[j] * molality[j] * molality[j] *
m_Mu_nnn_P[j];
5389 sum_m_phi_minus_1 = 2.0 *
5390 (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
5396 if (molalitysum > 1.0E-150) {
5397 d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
5399 d_osmotic_coef_dP = 0.0;
5403 printf(
" term1=%10.6f sum1=%10.6f sum2=%10.6f "
5404 "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
5405 term1, sum1, sum2, sum3, sum4, sum5);
5406 printf(
" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n",
5407 sum_m_phi_minus_1, d_osmotic_coef_dP);
5411 printf(
" Step 10: \n");
5414 d_lnwateract_dP = -(
m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
5429 double d_wateract_dP = exp(d_lnwateract_dP);
5430 printf(
" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
5431 d_lnwateract_dP, d_wateract_dP);
5439 double aphi, dj, jfunc, jprime, t, x, zprod;
5447 double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
5452 printf(
" Is = %g\n", is);
5455 if (is < 1.0E-150) {
5456 for (i = 0; i < 17; i++) {
5467 for (i=1; i<=4; i++) {
5468 for (j=i; j<=4; j++) {
5477 x = 6.0* zprod * aphi * sqrt(is);
5479 jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4)));
5481 t = c3 * c4 * pow(x,c4);
5482 dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
5483 jprime = (jfunc/x)*(1.0 + jfunc*dj);
5485 elambda[ij] = zprod*jfunc / (4.0*is);
5486 elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
5490 printf(
" ij = %d, elambda = %g, elambda1 = %g\n",
5499 double* etheta,
double* etheta_prime)
const
5512 if (i > 4 || j > 4) {
5513 printf(
"we shouldn't be here\n");
5518 if ((i == 0) || (j == 0)) {
5519 printf(
"ERROR calc_thetas called with one species being neutral\n");
5529 *etheta_prime = 0.0;
5535 f1 = (double)i / (2.0 * j);
5536 f2 = (double)j / (2.0 * i);
5555 for (
size_t k = 1; k <
m_kk; k++) {
5562 for (
size_t k = 1; k <
m_kk; k++) {
5569 for (
size_t k = 1; k <
m_kk; k++) {
5579 double xminus2 = xminus * xminus;
5580 double xminus3 = xminus2 * xminus;
5584 double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
5585 double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
5587 double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
5588 double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
5594 double h1_f = h1 * alpha;
5595 double h1_f_prime = h1_prime * alpha;
5597 double f = h2 + h1_f;
5598 double f_prime = h2_prime + h1_f_prime;
5600 double g = h2 + h1_g;
5601 double g_prime = h2_prime + h1_g_prime;
5603 tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
5604 double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5605 double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5607 tmp = log(xmolSolvent) + lngammak;
5608 for (
size_t k = 1; k <
m_kk; k++) {
5617 for (
size_t k = 1; k <
m_kk; k++) {
5625 double eterm = std::exp(-xoverc);
5627 double fptmp = IMS_bfCut_ - IMS_afCut_ /
IMS_cCut_ - IMS_bfCut_*xoverc
5628 + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
5629 double f_prime = 1.0 + eterm*fptmp;
5630 double f = xmolSolvent + IMS_efCut_
5631 + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
5633 double gptmp = IMS_bgCut_ - IMS_agCut_ /
IMS_cCut_ - IMS_bgCut_*xoverc
5634 + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
5635 double g_prime = 1.0 + eterm*gptmp;
5636 double g = xmolSolvent + IMS_egCut_
5637 + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
5639 tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
5640 double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5641 double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5643 tmp = log(xx) + lngammak;
5644 for (
size_t k = 1; k <
m_kk; k++) {
5656 std::string sni, snj;
5667 printf(
"Index Name MoleF MolalityCropped Charge\n");
5668 for (k = 0; k <
m_kk; k++) {
5670 printf(
"%2s %-16s %14.7le %14.7le %5.1f \n",
5671 int2str(k).c_str(), sni.c_str(), moleF[k], molality[k],
charge(k));
5674 printf(
"\n Species Species beta0MX "
5675 "beta1MX beta2MX CphiMX alphaMX thetaij \n");
5676 for (i = 1; i < m_kk - 1; i++) {
5678 for (j = i+1; j <
m_kk; j++) {
5680 size_t n = i * m_kk + j;
5682 printf(
" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
5683 sni.c_str(), snj.c_str(),
5692 printf(
"\n Species Species Species "
5694 for (i = 1; i <
m_kk; i++) {
5696 for (j = 1; j <
m_kk; j++) {
5698 for (k = 1; k <
m_kk; k++) {
5700 size_t n = k + j * m_kk + i * m_kk *
m_kk;
5702 printf(
" %-16s %-16s %-16s %9.5f \n",
5703 sni.c_str(), snj.c_str(),
5719 doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5720 for (
size_t k = 0; k <
m_kk; k++) {
5721 acMolality[k] *= exp(
charge(k) * afac);
5734 doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5735 for (
size_t k = 0; k <
m_kk; k++) {
5749 doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
5750 for (
size_t k = 0; k <
m_kk; k++) {
5764 doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
5765 for (
size_t k = 0; k <
m_kk; k++) {
5779 doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
5780 for (
size_t k = 0; k <
m_kk; k++) {
5789 doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
5790 return lnGammaClMs2;
5797 return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
5804 return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
5811 return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P.
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T.
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P.
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT.
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
virtual doublereal pressure() const
Pressure.
virtual doublereal density() const
Density (kg/m^3).
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
vector_fp m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
virtual void getParameters(int &n, doublereal *const c) const
Get the equation of state parameters in a vector.
WaterProps * m_waterProps
Pointer to the water property calculator.
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P.
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T.
XML_Node * findXMLPhase(XML_Node *root, const std::string &idtarget)
Search an XML_Node tree for a named phase XML_Node.
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT.
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
virtual doublereal logStandardConc(size_t k=0) const
Returns the natural logarithm of the standard concentration of the kth species.
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P.
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
const size_t npos
index returned by functions to indicate "no position"
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T.
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
int IMS_typeCutoff_
IMS Cutoff type.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
virtual doublereal satPressure(doublereal T)
saturation pressure
void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
size_t m_indexCLM
Index of the phScale species.
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P.
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T.
double elambda[17]
This is elambda, MEC.
Class XML_Node is a tree-based representation of the contents of an XML file.
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
int debugPrinting()
Return int specifying the amount of debug printing.
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P.
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P.
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
vector_fp m_speciesSize
Vector of species sizes.
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Header for a class used to house several approximation routines for properties of water...
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
doublereal m_xmolSolventMIN
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
doublereal molarDensity() const
Molar density (kmol/m^3).
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
double elambda1[17]
This is elambda1, MEC.
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P.
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT.
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT.
ThermoPhase * duplMyselfAsThermoPhase() const
Duplicator from the ThermoPhase parent class.
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
virtual void setTemperature(const doublereal temp)
Set the temperature (K)
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
virtual doublereal intEnergy_mole() const
Molar internal energy. Units: J/kmol.
virtual doublereal thermalExpansionCoeff() const
The thermal expansion coefficient.
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Base class for a phase with thermodynamic properties.
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
HMWSoln()
Default Constructor.
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
doublereal MC_X_o_min_
gamma_o value for the cutoff process at the zero solvent point
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P.
PDSS * m_waterSS
Water standard state calculator.
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
vector_fp m_Alpha1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty thermophase object.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T.
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature...
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients...
vector_fp m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
vector_fp m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
doublereal molarVolume() const
Molar volume (m^3/kmol).
The WaterProps class is used to house several approximation routines for properties of water...
#define PITZERFORM_BASE
Major Parameters: The form of the Pitzer expression refers to the form of the Gibbs free energy expre...
vector_fp m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT.
vector_fp m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT.
virtual void getUnitsStandardConc(double *uA, int k=0, int sizeUA=6) const
Returns the units of the standard and generalized concentrations.
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
MolalityVPSSTP & operator=(const MolalityVPSSTP &b)
Assignment operator.
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
doublereal IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
doublereal IMS_slopefCut_
Parameter in the polyExp cutoff treatment.
vector_fp m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT.
virtual void setParameters(int n, doublereal *const c)
Set the equation of state parameters.
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients...
vector_int m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species...
doublereal MC_slopepCut_
Parameter in the Molality Exp cutoff treatment.
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature...
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Base class for exceptions thrown by Cantera classes.
vector_fp m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
vector_fp m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T.
void calc_lambdas(double is) const
Calculate the lambda interactions.
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer...
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
vector_fp m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T.
vector_fp m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers.
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
virtual doublereal density() const
Return the standard state density at standard state.
double m_densWaterSS
density of standard-state water
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
virtual ~HMWSoln()
Destructor.
vector_fp m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT.
int m_formPitzer
This is the form of the Pitzer parameterization used in this model.
doublereal IMS_gamma_k_min_
gamma_k minimum for the cutoff process at the zero solvent point
vector_fp m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the Debye Huckel constant as a function of temperature and pressure.
vector_int m_electrolyteSpeciesType
Vector containing the electrolyte species type.
vector_fp m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
vector_fp m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T.
size_t nSpecies() const
Returns the number of species in the phase.
const int cHMWSoln0
eosTypes returned for this ThermoPhase Object
doublereal IMS_gamma_o_min_
gamma_o value for the cutoff process at the zero solvent point
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature as a function of tem...
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
vector_fp m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT.
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
doublereal temperature() const
Temperature (K).
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
int m_formGC
Format for the generalized concentration:
vector_fp m_molalities
Current value of the molalities of the species in the phase.
double m_A_Debye
A_Debye -> this expression appears on the top of the ln actCoeff term in the general Debye-Huckel exp...
std::vector< int > CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
void initLengths()
Initialize all of the species-dependent lengths in the object.
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized activity concentrations.
const doublereal SmallNumber
smallest number to compare to zero.
doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale...
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
virtual void applyphScale(doublereal *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
vector_fp m_gamma_tmp
Intermediate storage of the activity coefficient itself.
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
virtual int eosType() const
Equation of state type flag.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
doublereal err(const std::string &msg) const
Local error routine.
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
doublereal m_weightSolvent
Molecular weight of the Solvent.
doublereal IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
void calcMolalitiesCropped() const
Calculate the cropped molalities.
vector_fp m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
vector_fp m_pp
Temporary array used in equilibrium calculations.
Class HMWSoln represents a dilute or concentrated liquid electrolyte phase which obeys the Pitzer for...
doublereal ADebye(doublereal T, doublereal P, int ifunc)
ADebye calculates the value of A_Debye as a function of temperature and pressure according to relatio...
vector_fp m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper This is the basic cation - anion interaction...
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
virtual void setState_TP(doublereal t, doublereal p)
Set the temperature (K) and pressure (Pa)
size_t m_kk
Number of species in the phase.
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
vector_fp m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
doublereal m_Pcurrent
Current value of the pressure - state variable.
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
vector_fp m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
vector_fp m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
virtual doublereal density() const
Returns the current value of the density.
size_t m_indexSolvent
Index of the solvent.
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model...
void printCoeffs() const
Print out all of the input Pitzer coefficients.
HMWSoln & operator=(const HMWSoln &right)
Assignment operator.
vector_fp m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T.
vector_fp m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT.
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual doublereal satPressure(doublereal T)
Get the saturation pressure for a given temperature.
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
vector_fp m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
std::string speciesName(size_t k) const
Name of the species with index k.
vector_fp m_Aionic
a_k = Size of the ionic species in the DH formulation units = meters
vector_fp m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P.
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
vector_fp m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase Note the density of a phase is an independent...
vector_fp m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P.
virtual doublereal isothermalCompressibility() const
The isothermal compressibility.
vector_fp m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT.
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...