LDBLE T = tc + 273.15; LDBLE u1 = 3.4279e2, u2 = -5.0866e-3, u3 = 9.469e-7, u4 = -2.0525, u5 = 3.1159e3, u6 = -1.8289e2, u7 = -8.0325e3, u8 = 4.2142e6, u9 = 2.1417; LDBLE d1000 = u1 * exp(T * (u2 + T * u3)); // relative dielectric constant at 1000 bar LDBLE c = u4 + u5 / (u6 + T); LDBLE b = u7 + u8 / T + u9 * T; LDBLE pb = pa * 1.01325; // pa in bar eps_r = d1000 + c * log((b + pb) / (b + 1e3)); // relative dielectric constant if (eps_r <= 0) { eps_r = 10.; warning_msg("Relative dielectric constant is negative.\nTemperature is out of range of parameterization."); } /* qe^2 / (eps_r * kB * T) = 4.803204e-10**2 / 1.38065e-16 / (eps_r * T) = 1.671008e-3 (esu^2 / (erg/K)) / (eps_r * T) */ LDBLE e2_DkT = 1.671008e-3 / (eps_r * T); DH_B = sqrt(8 * pi * AVOGADRO * e2_DkT * rho_0 / 1e3); // Debye length parameter, 1/cm(mol/kg)^-0.5 DH_A = DH_B * e2_DkT / (2. * LOG_10); //(mol/kg)^-0.5 /* A0 in pitzer */ if (pitzer_model || sit_model) { A0 = DH_B * e2_DkT / 6.0; if (pitzer_model && aphi != NULL) { calc_pitz_param(aphi, T, 298.15); A0 = aphi->p; } } /* Debye-Hueckel limiting slope = DH_B * e2_DkT * RT * (d(ln(eps_r)) / d(P) - compressibility) */ DH_Av = DH_B * e2_DkT * R_LITER_ATM * 1e3 * T * (c / (b + pb) * 1.01325 / eps_r - kappa_0 / 3.); // (cm3/mol)(mol/kg)^-0.5 DH_B /= 1e8; // kappa, 1/Angstrom(mol/kg)^-0.5 /* the Born functions, * 41.84 to give molal volumes in cm3/mol... */ ZBrn = (- 1 / eps_r + 1.0) * 41.84004; QBrn = c / (b + pb) / eps_r / eps_r * 41.84004;