mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-16 08:38:23 +01:00
Optimizations for pitzer.cpp
Change in PTEMP did not work for unknown reasons. git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/trunk@9793 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst
parent
0e78b17230
commit
4e07ab28e7
@ -524,6 +524,7 @@ public:
|
||||
struct theta_param *theta_param_search(LDBLE zj, LDBLE zk);
|
||||
struct theta_param *theta_param_alloc(void);
|
||||
int theta_param_init(struct theta_param *theta_param_ptr);
|
||||
void pitzer_make_lists(void);
|
||||
int gammas_pz(void);
|
||||
int model_pz(void);
|
||||
int pitzer(void);
|
||||
@ -543,12 +544,13 @@ public:
|
||||
LDBLE GP(LDBLE Y);
|
||||
int ETHETAS(LDBLE ZJ, LDBLE ZK, LDBLE I, LDBLE * etheta,
|
||||
LDBLE * ethetap);
|
||||
int BDK(LDBLE X);
|
||||
void ETHETA_PARAMS(LDBLE X, LDBLE& JAY, LDBLE& JPRIME );
|
||||
//int BDK(LDBLE X);
|
||||
int pitzer_initial_guesses(void);
|
||||
int pitzer_revise_guesses(void);
|
||||
int PTEMP(LDBLE TK);
|
||||
LDBLE JAY(LDBLE X);
|
||||
LDBLE JPRIME(LDBLE Y);
|
||||
//LDBLE JAY(LDBLE X);
|
||||
//LDBLE JPRIME(LDBLE Y);
|
||||
int jacobian_pz(void);
|
||||
|
||||
// pitzer_structures.cpp -------------------------------
|
||||
|
||||
575
pitzer.cpp
575
pitzer.cpp
@ -770,7 +770,9 @@ C
|
||||
DW0 = rho_0 = calc_rho_0(TK - 273.15, patm_x);
|
||||
VP = patm_x;
|
||||
for (i = 0; i < count_pitz_param; i++)
|
||||
//for (size_t j = 0; j < param_list.size(); j++)
|
||||
{
|
||||
//int i = param_list[j];
|
||||
calc_pitz_param(pitz_params[i], TK, TR);
|
||||
}
|
||||
calc_dielectrics(TK - 273.15, patm_x);
|
||||
@ -845,7 +847,7 @@ calc_pitz_param(struct pitz_param *pz_ptr, LDBLE TK, LDBLE TR)
|
||||
}
|
||||
return OK;
|
||||
}
|
||||
|
||||
#ifdef SKIP
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
pitzer(void)
|
||||
@ -1180,7 +1182,385 @@ pitzer(void)
|
||||
*/
|
||||
return (OK);
|
||||
}
|
||||
#endif
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
pitzer(void)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, i0, i1, i2;
|
||||
LDBLE param, l_alpha, z0, z1;
|
||||
LDBLE etheta, ethetap;
|
||||
/*
|
||||
LDBLE CONV, XI, XX, OSUM, BIGZ, DI, F, XXX, GAMCLM,
|
||||
CSUM, PHIMAC, OSMOT, BMXP, ETHEAP, CMX, BMX, PHI,
|
||||
BMXPHI, PHIPHI, AW, A, B;
|
||||
*/
|
||||
LDBLE CONV, XX, OSUM, BIGZ, DI, F, F1, F2, F_var, XXX, GAMCLM, CSUM, PHIMAC, OSMOT,
|
||||
B, B1, B2;
|
||||
LDBLE I, TK;
|
||||
/*
|
||||
C
|
||||
C INITIALIZE
|
||||
C
|
||||
*/
|
||||
CONV = 1.0 / log(10.0);
|
||||
XX = 0.0;
|
||||
OSUM = 0.0;
|
||||
/*n
|
||||
I = *I_X;
|
||||
TK = *TK_X;
|
||||
*/
|
||||
I = mu_x;
|
||||
TK = tk_x;
|
||||
/* DH_AB(TK, &A, &B); */
|
||||
/*
|
||||
C
|
||||
C TRANSFER DATA FROM TO M
|
||||
C
|
||||
*/
|
||||
double log_min = log10(MIN_TOTAL);
|
||||
for (size_t j = 0; j < s_list.size(); j++)
|
||||
{
|
||||
i = s_list[j];
|
||||
if (spec[i]->lm > log_min)
|
||||
{
|
||||
M[i] = under(spec[i]->lm);
|
||||
}
|
||||
else
|
||||
{
|
||||
M[i] = 0.0;
|
||||
}
|
||||
}
|
||||
#ifdef SKIP
|
||||
for (i = 0; i < 3 * count_s; i++)
|
||||
{
|
||||
IPRSNT[i] = FALSE;
|
||||
M[i] = 0.0;
|
||||
if (spec[i] != NULL && spec[i]->in == TRUE)
|
||||
{
|
||||
if (spec[i]->type == EX ||
|
||||
spec[i]->type == SURF || spec[i]->type == SURF_PSI)
|
||||
continue;
|
||||
M[i] = under(spec[i]->lm);
|
||||
if (M[i] > MIN_TOTAL)
|
||||
IPRSNT[i] = TRUE;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
if (ICON == TRUE)
|
||||
{
|
||||
IPRSNT[IC] = TRUE;
|
||||
}
|
||||
/*
|
||||
ICON = 0;
|
||||
M[1] = 1.40070736;
|
||||
M[4] = 2.52131086E-05;
|
||||
M[140] = 4.59985435E-09;
|
||||
*/
|
||||
|
||||
/*
|
||||
C
|
||||
C COMPUTE PITZER COEFFICIENTS' TEMPERATURE DEPENDENCE
|
||||
C
|
||||
*/
|
||||
PTEMP(TK);
|
||||
for (size_t j = 0; j < s_list.size(); j++)
|
||||
{
|
||||
int i = s_list[j];
|
||||
LGAMMA[i] = 0.0;
|
||||
XX = XX + M[i] * fabs(spec[i]->z);
|
||||
OSUM = OSUM + M[i];
|
||||
}
|
||||
//for (i = 0; i < 2 * count_s + count_anions; i++)
|
||||
//{
|
||||
// LGAMMA[i] = 0.0;
|
||||
// if (IPRSNT[i] == TRUE)
|
||||
// {
|
||||
// XX = XX + M[i] * fabs(spec[i]->z);
|
||||
// OSUM = OSUM + M[i];
|
||||
// }
|
||||
//}
|
||||
/*
|
||||
C
|
||||
C EQUATION (8)
|
||||
C
|
||||
*/
|
||||
BIGZ = XX;
|
||||
DI = sqrt(I);
|
||||
/*
|
||||
C
|
||||
C CALCULATE F & GAMCLM
|
||||
C
|
||||
*/
|
||||
B = 1.2;
|
||||
F = F1 = F2 = -A0 * (DI / (1.0 + B * DI) + 2.0 * log(1.0 + B * DI) / B);
|
||||
if (patm_x > 1.0)
|
||||
{
|
||||
LDBLE pap;
|
||||
pap = (7e-5 + 1.93e-9 * pow(TK - 250.0, 2.0)) * patm_x;
|
||||
B1 = B - (pap > 0.2 ? 0.2 : pap);
|
||||
pap = (9.65e-10 * pow(TK - 263.0, 2.773)) * pow(patm_x, 0.623);
|
||||
//pap = (-5.22e-4 + 7.19e-8 * pow(TK - 263.0, 2.0)) * pow(patm_x, 0.623);
|
||||
B2 = B - (pap > 0.2 ? 0.2 : pap);
|
||||
if (B1 != 0)
|
||||
F1 = -A0 * (DI / (1.0 + B1 * DI) + 2.0 * log(1.0 + B1 * DI) / B1);
|
||||
if (B2 != 0)
|
||||
F2 = -A0 * (DI / (1.0 + B2 * DI) + 2.0 * log(1.0 + B2 * DI) / B2);
|
||||
}
|
||||
XXX = 2.0 * DI;
|
||||
XXX = (1.0 - (1.0 + XXX - XXX * XXX * 0.5) * exp(-XXX)) / (XXX * XXX);
|
||||
/*GAMCLM=F+I*2.0e0*(BCX(1,IK,IC)+BCX(2,IK,IC)*XXX)+1.5e0*BCX(4,IK,IC)*I*I; */
|
||||
/*GAMCLM=F+I*2.0e0*(mcb0->U.b0 + mcb1->U.b1*XXX) + 1.5e0*mcc0->U.c0*I*I; */
|
||||
/*GAMCLM = F + I * 2.0e0 * (mcb0->p + mcb1->p * XXX) + 1.5e0 * mcc0->p * I * I; */
|
||||
GAMCLM = F1;
|
||||
if (mcb0 != NULL)
|
||||
GAMCLM += I * 2.0 * mcb0->p;
|
||||
if (mcb1 != NULL)
|
||||
GAMCLM += I * 2.0 * mcb1->p * XXX;
|
||||
if (mcc0 != NULL)
|
||||
GAMCLM += 1.5 * mcc0->p * I * I;
|
||||
CSUM = 0.0;
|
||||
OSMOT = -(A0) * pow(I, (LDBLE) 1.5) / (1.0 + B * DI);
|
||||
/*
|
||||
* Calculate ethetas
|
||||
*/
|
||||
if (use_etheta == TRUE)
|
||||
{
|
||||
for (i = 0; i < count_theta_param; i++)
|
||||
{
|
||||
z0 = theta_params[i]->zj;
|
||||
z1 = theta_params[i]->zk;
|
||||
ETHETAS(z0, z1, I, ðeta, ðetap);
|
||||
theta_params[i]->etheta = etheta;
|
||||
theta_params[i]->ethetap = ethetap;
|
||||
}
|
||||
}
|
||||
/*
|
||||
* Sums for F, LGAMMA, and OSMOT
|
||||
*/
|
||||
//for (i = 0; i < count_pitz_param; i++)
|
||||
for (size_t j = 0; j < param_list.size(); j++)
|
||||
{
|
||||
int i = param_list[j];
|
||||
i0 = pitz_params[i]->ispec[0];
|
||||
i1 = pitz_params[i]->ispec[1];
|
||||
//if (IPRSNT[i0] == FALSE || IPRSNT[i1] == FALSE)
|
||||
// continue;
|
||||
z0 = spec[i0]->z;
|
||||
z1 = spec[i1]->z;
|
||||
param = pitz_params[i]->p;
|
||||
l_alpha = pitz_params[i]->alpha;
|
||||
F_var = 0;
|
||||
switch (pitz_params[i]->type)
|
||||
{
|
||||
case TYPE_B0:
|
||||
LGAMMA[i0] += M[i1] * 2.0 * param;
|
||||
LGAMMA[i1] += M[i0] * 2.0 * param;
|
||||
OSMOT += M[i0] * M[i1] * param;
|
||||
break;
|
||||
case TYPE_B1:
|
||||
if (param != 0.0)
|
||||
{
|
||||
F_var = M[i0] * M[i1] * param * GP(l_alpha * DI) / I;
|
||||
LGAMMA[i0] += M[i1] * 2.0 * param * G(l_alpha * DI);
|
||||
LGAMMA[i1] += M[i0] * 2.0 * param * G(l_alpha * DI);
|
||||
OSMOT += M[i0] * M[i1] * param * exp(-l_alpha * DI);
|
||||
}
|
||||
break;
|
||||
case TYPE_B2:
|
||||
if (param != 0.0)
|
||||
{
|
||||
F_var = M[i0] * M[i1] * param * GP(l_alpha * DI) / I;
|
||||
LGAMMA[i0] += M[i1] * 2.0 * param * G(l_alpha * DI);
|
||||
LGAMMA[i1] += M[i0] * 2.0 * param * G(l_alpha * DI);
|
||||
OSMOT += M[i0] * M[i1] * param * exp(-l_alpha * DI);
|
||||
}
|
||||
break;
|
||||
case TYPE_C0:
|
||||
CSUM +=
|
||||
M[i0] * M[i1] * pitz_params[i]->p / (2.0 *
|
||||
sqrt(fabs(z0 * z1)));
|
||||
LGAMMA[i0] += M[i1] * BIGZ * param / (2.0 * sqrt(fabs(z0 * z1)));
|
||||
LGAMMA[i1] += M[i0] * BIGZ * param / (2.0 * sqrt(fabs(z0 * z1)));
|
||||
OSMOT +=
|
||||
M[i0] * M[i1] * BIGZ * param / (2.0 * sqrt(fabs(z0 * z1)));
|
||||
break;
|
||||
case TYPE_THETA:
|
||||
LGAMMA[i0] += 2.0 * M[i1] * (param /*+ ETHETA(z0, z1, I) */ );
|
||||
LGAMMA[i1] += 2.0 * M[i0] * (param /*+ ETHETA(z0, z1, I) */ );
|
||||
OSMOT += M[i0] * M[i1] * param;
|
||||
break;
|
||||
case TYPE_ETHETA:
|
||||
/*
|
||||
ETHETAS(z0, z1, I, ðeta, ðetap);
|
||||
*/
|
||||
if (use_etheta == TRUE)
|
||||
{
|
||||
etheta = pitz_params[i]->thetas->etheta;
|
||||
ethetap = pitz_params[i]->thetas->ethetap;
|
||||
F_var = M[i0] * M[i1] * ethetap;
|
||||
LGAMMA[i0] += 2.0 * M[i1] * etheta;
|
||||
LGAMMA[i1] += 2.0 * M[i0] * etheta;
|
||||
OSMOT += M[i0] * M[i1] * (etheta + I * ethetap);
|
||||
/*
|
||||
F += M[i0]*M[i1]*ETHETAP(z0, z1, I);
|
||||
LGAMMA[i0] += 2.0*M[i1]*(ETHETA(z0, z1, I) );
|
||||
LGAMMA[i1] += 2.0*M[i0]*(ETHETA(z0, z1, I) );
|
||||
OSMOT += M[i0]*M[i1]*(ETHETA(z0, z1, I) + I*ETHETAP(z0, z1, I) );
|
||||
*/
|
||||
}
|
||||
break;
|
||||
case TYPE_PSI:
|
||||
i2 = pitz_params[i]->ispec[2];
|
||||
if (IPRSNT[i2] == FALSE)
|
||||
continue;
|
||||
LGAMMA[i0] += M[i1] * M[i2] * param;
|
||||
LGAMMA[i1] += M[i0] * M[i2] * param;
|
||||
LGAMMA[i2] += M[i0] * M[i1] * param;
|
||||
OSMOT += M[i0] * M[i1] * M[i2] * param;
|
||||
break;
|
||||
case TYPE_LAMDA:
|
||||
LGAMMA[i0] += M[i1] * param * pitz_params[i]->ln_coef[0];
|
||||
LGAMMA[i1] += M[i0] * param * pitz_params[i]->ln_coef[1];
|
||||
OSMOT += M[i0] * M[i1] * param * pitz_params[i]->os_coef;
|
||||
break;
|
||||
case TYPE_ZETA:
|
||||
i2 = pitz_params[i]->ispec[2];
|
||||
if (IPRSNT[i2] == FALSE)
|
||||
continue;
|
||||
LGAMMA[i0] += M[i1] * M[i2] * param;
|
||||
LGAMMA[i1] += M[i0] * M[i2] * param;
|
||||
LGAMMA[i2] += M[i0] * M[i1] * param;
|
||||
OSMOT += M[i0] * M[i1] * M[i2] * param;
|
||||
break;
|
||||
case TYPE_MU:
|
||||
i2 = pitz_params[i]->ispec[2];
|
||||
if (IPRSNT[i2] == FALSE)
|
||||
continue;
|
||||
|
||||
LGAMMA[i0] += M[i1] * M[i2] * param * pitz_params[i]->ln_coef[0];
|
||||
LGAMMA[i1] += M[i0] * M[i2] * param * pitz_params[i]->ln_coef[1];
|
||||
LGAMMA[i2] += M[i0] * M[i1] * param * pitz_params[i]->ln_coef[2];
|
||||
OSMOT += M[i0] * M[i1] * M[i2] * param * pitz_params[i]->os_coef;
|
||||
break;
|
||||
case TYPE_ETA:
|
||||
i2 = pitz_params[i]->ispec[2];
|
||||
if (IPRSNT[i2] == FALSE)
|
||||
continue;
|
||||
LGAMMA[i0] += M[i1] * M[i2] * param;
|
||||
LGAMMA[i1] += M[i0] * M[i2] * param;
|
||||
LGAMMA[i2] += M[i0] * M[i1] * param;
|
||||
OSMOT += M[i0] * M[i1] * M[i2] * param;
|
||||
break;
|
||||
case TYPE_ALPHAS:
|
||||
break;
|
||||
case TYPE_Other:
|
||||
default:
|
||||
error_msg("TYPE_Other in pitz_param list.", STOP);
|
||||
break;
|
||||
}
|
||||
F += F_var;
|
||||
F1 += F_var;
|
||||
F2 += F_var;
|
||||
}
|
||||
|
||||
/*
|
||||
* Add F and CSUM terms to LGAMMA
|
||||
*/
|
||||
|
||||
for (size_t j = 0; j < ion_list.size(); j++)
|
||||
{
|
||||
int i = ion_list[j];
|
||||
z0 = fabs(spec[i]->z);
|
||||
F_var = (z0 == 1 ? F1 : (z0 == 2.0 ? F2 : F));
|
||||
LGAMMA[i] += z0 * z0 * F_var + z0 * CSUM;
|
||||
}
|
||||
//for (i = 0; i < count_cations; i++)
|
||||
//{
|
||||
// if (!IPRSNT[i])
|
||||
// continue;
|
||||
// z0 = fabs(spec[i]->z);
|
||||
// F_var = (z0 == 1 ? F1 : (z0 == 2.0 ? F2 : F));
|
||||
// LGAMMA[i] += z0 * z0 * F_var + z0 * CSUM;
|
||||
//}
|
||||
//for (i = 2 * count_s; i < 2 * count_s + count_anions; i++)
|
||||
//{
|
||||
// if (!IPRSNT[i])
|
||||
// continue;
|
||||
// z0 = fabs(spec[i]->z);
|
||||
// F_var = (z0 == 1 ? F1 : (z0 == 2.0 ? F2 : F));
|
||||
// LGAMMA[i] += z0 * z0 * F_var + z0 * CSUM;
|
||||
//}
|
||||
/*
|
||||
C
|
||||
C CONVERT TO MACINNES CONVENTION
|
||||
C
|
||||
*/
|
||||
if (ICON == TRUE)
|
||||
{
|
||||
PHIMAC = LGAMMA[IC] - GAMCLM;
|
||||
/*
|
||||
C
|
||||
C CORRECTED ERROR IN PHIMAC, NOVEMBER, 1989
|
||||
C
|
||||
*/
|
||||
for (size_t j = 0; j < s_list.size(); j++)
|
||||
{
|
||||
int i = s_list[j];
|
||||
LGAMMA[i] = LGAMMA[i] + spec[i]->z * PHIMAC;
|
||||
}
|
||||
//for (i = 0; i < 2 * count_s + count_anions; i++)
|
||||
//{
|
||||
// if (IPRSNT[i] == TRUE)
|
||||
// {
|
||||
// LGAMMA[i] = LGAMMA[i] + spec[i]->z * PHIMAC;
|
||||
// }
|
||||
//}
|
||||
}
|
||||
|
||||
COSMOT = 1.0 + 2.0 * OSMOT / OSUM;
|
||||
/*
|
||||
C
|
||||
C CALCULATE THE ACTIVITY OF WATER
|
||||
C
|
||||
*/
|
||||
AW = exp(-OSUM * COSMOT / 55.50837);
|
||||
/*
|
||||
if (AW > 1.0)
|
||||
AW = 1.0;
|
||||
*/
|
||||
/*s_h2o->la=log10(AW); */
|
||||
mu_x = I;
|
||||
for (size_t j = 0; j < s_list.size(); j++)
|
||||
{
|
||||
int i = s_list[j];
|
||||
spec[i]->lg_pitzer = LGAMMA[i] * CONV;
|
||||
}
|
||||
//for (i = 0; i < 2 * count_s + count_anions; i++)
|
||||
//{
|
||||
// if (IPRSNT[i] == FALSE)
|
||||
// continue;
|
||||
// /*spec[i]->lg=LGAMMA[i]*CONV; */
|
||||
// spec[i]->lg_pitzer = LGAMMA[i] * CONV;
|
||||
// /*
|
||||
// output_msg(sformatf( "%d %s:\t%e\t%e\t%e\t%e \n", i, spec[i]->name, M[i], spec[i]->la, spec[i]->lg_pitzer, spec[i]->lg));
|
||||
// */
|
||||
//}
|
||||
/*
|
||||
output_msg(sformatf( "OSUM: %e\n", OSUM));
|
||||
output_msg(sformatf( "OSMOT: %e\n", OSMOT));
|
||||
output_msg(sformatf( "COSMOT: %e\n", COSMOT));
|
||||
output_msg(sformatf( "F: %e\n", F));
|
||||
output_msg(sformatf( "AW: %e\n", AW));
|
||||
*/
|
||||
/*
|
||||
*I_X = I;
|
||||
*COSMOT_X = COSMOT;
|
||||
*/
|
||||
return (OK);
|
||||
}
|
||||
#ifdef SKIP
|
||||
/* ---------------------------------------------------------------------- */
|
||||
LDBLE Phreeqc::
|
||||
JAY(LDBLE X)
|
||||
@ -1217,7 +1597,6 @@ JPRIME(LDBLE L_Y)
|
||||
return (L_Y * (.25e0 + L_DZ * (DK[0] - DK[2]) / 2.0e0));
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
BDK(LDBLE X)
|
||||
@ -1277,6 +1656,7 @@ C
|
||||
}
|
||||
return OK;
|
||||
}
|
||||
#endif
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
LDBLE Phreeqc::
|
||||
@ -1305,7 +1685,7 @@ GP(LDBLE L_Y)
|
||||
}
|
||||
return d;
|
||||
}
|
||||
|
||||
#ifdef SKIP
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
ETHETAS(LDBLE ZJ, LDBLE ZK, LDBLE I, LDBLE * etheta, LDBLE * ethetap)
|
||||
@ -1340,6 +1720,124 @@ C
|
||||
JPRIME(XKK) / 2.0e0) / (8.0e0 * I * I) - *etheta / I;
|
||||
return (OK);
|
||||
}
|
||||
#endif
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
ETHETAS(LDBLE ZJ, LDBLE ZK, LDBLE I, LDBLE * etheta, LDBLE * ethetap)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
*etheta = 0.0;
|
||||
*ethetap = 0.0;
|
||||
|
||||
if (ZJ == ZK)
|
||||
return (OK);
|
||||
|
||||
const LDBLE XCON = 6.0e0 * A0 * sqrt(I);
|
||||
const LDBLE ZZ = ZJ * ZK;
|
||||
/*
|
||||
C
|
||||
C NEXT 3 ARE EQUATION (A1)
|
||||
C
|
||||
*/
|
||||
const LDBLE XJK = XCON * ZZ;
|
||||
const LDBLE XJJ = XCON * ZJ * ZJ;
|
||||
const LDBLE XKK = XCON * ZK * ZK;
|
||||
|
||||
/*
|
||||
C
|
||||
C EQUATION (A3)
|
||||
C
|
||||
*/
|
||||
LDBLE JAY_XJK;
|
||||
LDBLE JPRIME_XJK;
|
||||
ETHETA_PARAMS( XJK, JAY_XJK, JPRIME_XJK );
|
||||
|
||||
LDBLE JAY_XJJ;
|
||||
LDBLE JPRIME_XJJ;
|
||||
ETHETA_PARAMS( XJJ, JAY_XJJ, JPRIME_XJJ );
|
||||
|
||||
LDBLE JAY_XKK;
|
||||
LDBLE JPRIME_XKK;
|
||||
ETHETA_PARAMS( XKK, JAY_XKK, JPRIME_XKK );
|
||||
|
||||
*etheta =
|
||||
ZZ * (JAY_XJK - JAY_XJJ / 2.0e0 - JAY_XKK / 2.0e0) / (4.0e0 * I);
|
||||
*ethetap =
|
||||
ZZ * (JPRIME_XJK - JPRIME_XJJ / 2.0e0 -
|
||||
JPRIME_XKK / 2.0e0) / (8.0e0 * I * I) - *etheta / I;
|
||||
|
||||
return (OK);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void Phreeqc::
|
||||
ETHETA_PARAMS(LDBLE X, LDBLE& JAY, LDBLE& JPRIME )
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/*
|
||||
C
|
||||
C NUMERICAL APPROXIMATION TO THE INTEGRALS IN THE EXPRESSIONS FOR J0
|
||||
C AND J1. CHEBYSHEV APPROXIMATION IS USED. THE CONSTANTS 'AK' ARE
|
||||
C DEFINED IN BLOCK COMMON.
|
||||
C
|
||||
*/
|
||||
/*
|
||||
C
|
||||
C AK IS USED TO CALCULATE HIGHER ORDER ELECTROSTATIC TERMS IN
|
||||
C SUBROUTINE PITZER
|
||||
C
|
||||
*/
|
||||
{
|
||||
static const LDBLE AKX[42] = {
|
||||
1.925154014814667e0, -.060076477753119e0, -.029779077456514e0,
|
||||
-.007299499690937e0, 0.000388260636404e0, 0.000636874599598e0,
|
||||
0.000036583601823e0, -.000045036975204e0, -.000004537895710e0,
|
||||
0.000002937706971e0, 0.000000396566462e0, -.000000202099617e0,
|
||||
-.000000025267769e0, 0.000000013522610e0, 0.000000001229405e0,
|
||||
-.000000000821969e0, -.000000000050847e0, 0.000000000046333e0,
|
||||
0.000000000001943e0, -.000000000002563e0, -.000000000010991e0,
|
||||
0.628023320520852e0, 0.462762985338493e0, 0.150044637187895e0,
|
||||
-.028796057604906e0, -.036552745910311e0, -.001668087945272e0,
|
||||
0.006519840398744e0, 0.001130378079086e0, -.000887171310131e0,
|
||||
-.000242107641309e0, 0.000087294451594e0, 0.000034682122751e0,
|
||||
-.000004583768938e0, -.000003548684306e0, -.000000250453880e0,
|
||||
0.000000216991779e0, 0.000000080779570e0, 0.000000004558555e0,
|
||||
-.000000006944757e0, -.000000002849257e0, 0.000000000237816e0
|
||||
};
|
||||
/*
|
||||
LDBLE PRECISION AK, BK, DK
|
||||
COMMON / MX8 / AK(0:20,2),BK(0:22),DK(0:22)
|
||||
*/
|
||||
const LDBLE *AK;
|
||||
LDBLE L_Z = 0.0;
|
||||
LDBLE L_DZ = 0.0;
|
||||
|
||||
if ( X <= 1.0e0 )
|
||||
{
|
||||
const LDBLE powX0_2 = pow( X, 0.2 );
|
||||
L_Z = 4.0e0 * powX0_2 - 2.0e0;
|
||||
L_DZ = 0.8e0 * powX0_2 / 2.0e0;
|
||||
AK = &AKX[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
const LDBLE powXmin0_1 = pow( X, -0.1 );
|
||||
L_Z = ( 40.0e0 * powXmin0_1 - 22.0e0 ) / 9.0e0;
|
||||
L_DZ = -4.0e0 * powXmin0_1 / 18.0e0;
|
||||
AK = &AKX[21];
|
||||
}
|
||||
|
||||
BK[20] = AK[20];
|
||||
BK[19] = L_Z * AK[20] + AK[19];
|
||||
DK[19] = AK[20];
|
||||
for ( int i = 18; i >= 0; i-- )
|
||||
{
|
||||
BK[i] = L_Z * BK[i + 1] - BK[i + 2] + AK[i];
|
||||
DK[i] = BK[i + 1] + L_Z * DK[i + 1] - DK[i + 2];
|
||||
}
|
||||
|
||||
JAY = X / 4.0e0 - 1.0e0 + 0.5e0 * (BK[0] - BK[2]);
|
||||
JPRIME = X * .25e0 + L_DZ * (DK[0] - DK[2]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
@ -2222,3 +2720,74 @@ gammas_pz()
|
||||
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void Phreeqc::
|
||||
pitzer_make_lists(void)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
double log_min = log10(MIN_TOTAL);
|
||||
s_list.clear();
|
||||
cation_list.clear();
|
||||
neutral_list.clear();
|
||||
anion_list.clear();
|
||||
ion_list.clear();
|
||||
param_list.clear();
|
||||
for (int i = 0; i < 3 * count_s; i++)
|
||||
{
|
||||
IPRSNT[i] = FALSE;
|
||||
M[i] = 0.0;
|
||||
if (spec[i] != NULL && spec[i]->in == TRUE)
|
||||
{
|
||||
if (spec[i]->type == EX ||
|
||||
spec[i]->type == SURF || spec[i]->type == SURF_PSI)
|
||||
continue;
|
||||
IPRSNT[i] = TRUE;
|
||||
s_list.push_back(i);
|
||||
if (i < count_s)
|
||||
{
|
||||
cation_list.push_back(i);
|
||||
}
|
||||
if (i >= count_s && i < 2*count_s)
|
||||
{
|
||||
neutral_list.push_back(i);
|
||||
}
|
||||
if (i >= 2*count_s)
|
||||
{
|
||||
anion_list.push_back(i);
|
||||
}
|
||||
if (i < count_s || i >= 2*count_s)
|
||||
{
|
||||
ion_list.push_back(i);
|
||||
}
|
||||
if (spec[i]->lm > log_min)
|
||||
{
|
||||
M[i] = under(spec[i]->lm);
|
||||
}
|
||||
}
|
||||
}
|
||||
if (ICON == TRUE)
|
||||
{
|
||||
IPRSNT[IC] = TRUE;
|
||||
}
|
||||
for (int i = 0; i < count_pitz_param; i++)
|
||||
{
|
||||
/*
|
||||
TYPE_B0, TYPE_B1, TYPE_B2, TYPE_C0, TYPE_THETA, TYPE_LAMDA, TYPE_ZETA,
|
||||
TYPE_PSI, TYPE_ETHETA, TYPE_ALPHAS, TYPE_MU, TYPE_ETA, TYPE_Other,
|
||||
TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU
|
||||
*/
|
||||
int i0 = pitz_params[i]->ispec[0];
|
||||
int i1 = pitz_params[i]->ispec[1];
|
||||
if (IPRSNT[i0] == FALSE || IPRSNT[i1] == FALSE) continue;
|
||||
int i2 = pitz_params[i]->ispec[2];
|
||||
if (pitz_params[i]->type == TYPE_PSI ||
|
||||
pitz_params[i]->type == TYPE_ZETA ||
|
||||
pitz_params[i]->type == TYPE_MU ||
|
||||
pitz_params[i]->type == TYPE_ETA)
|
||||
{
|
||||
if (IPRSNT[i2] == FALSE)
|
||||
continue;
|
||||
}
|
||||
param_list.push_back(i);
|
||||
}
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user