#include "Phreeqc.h"
#include "phqalloc.h"
#include "Utils.h"
#include "Solution.h"

#define MAX_QUAD 20
#define K_POLY 5

#if defined(PHREEQCI_GUI)
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#endif

/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_all_g(void)
/* ---------------------------------------------------------------------- */
{
	int converge, converge1;
	LDBLE new_g, xd1;
	LDBLE epsilon;

	if (use.Get_surface_ptr() == NULL)
		return (OK);
	/*
	 *   calculate g for each surface
	 */
	epsilon = convergence_tolerance;
	if (convergence_tolerance >= 1e-8)
	{
		G_TOL = 1e-9;
	}
	else
	{
		G_TOL = 1e-10;
	}

	converge = TRUE;

	for (int j = 0; j < count_unknowns; j++)
	{
		if (x[j]->type != SURFACE_CB)
			continue;
		if (debug_diffuse_layer == TRUE)
			output_msg(sformatf("Calc_all_g, X[%d]\n", j));
		cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
		std::map<LDBLE, cxxSurfDL> temp_g_map;
		cxxSurfDL temp_g;
		charge_ptr->Get_g_map()[0] = temp_g;
		temp_g_map[0] = temp_g;
		xd_global = exp(-2 * x[j]->master[0]->s->la * LOG_10);
		/* alpha = 0.02935 @ 25;        (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
		/* 1000 J/kJ and 1000 L/m**3 */
		//alpha_global =	sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
		//		 tk_x * 0.5);
		alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
			tk_x * 0.5);
		/*
		 *   calculate g for given surface for each species
		 */
		for (int i = 0; i < (int)this->s_x.size(); i++)
		{
			if (s_x[i]->type > HPLUS)
				continue;
			if (temp_g_map.find(s_x[i]->z) != temp_g_map.end())
				continue;
			z_global = s_x[i]->z;
			if (charge_ptr->Get_grams() > 0.0)
			{

				if ((use.Get_surface_ptr()->Get_only_counter_ions() == false) ||
					(((x[j]->master[0]->s->la > 0) && (z_global < 0))
						|| ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
				{
					if (xd_global > 0.1)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, xd_global);
					}
					else if (xd_global > 0.01)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, xd_global);
					}
					else if (xd_global > 0.001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, xd_global);
					}
					else if (xd_global > 0.0001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, xd_global);
					}
					else if (xd_global > 0.00001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
						new_g += qromb_midpnt(charge_ptr, 0.0001, xd_global);
					}
					else if (xd_global > 0.000001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
						new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
						new_g += qromb_midpnt(charge_ptr, 0.00001, xd_global);
					}
					else if (xd_global > 0.0000001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
						new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
						new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
						new_g += qromb_midpnt(charge_ptr, 0.000001, xd_global);
					}
					else if (xd_global > 0.00000001)
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
						new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
						new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
						new_g += qromb_midpnt(charge_ptr, 0.000001, .0000001);
						new_g += qromb_midpnt(charge_ptr, 0.0000001, xd_global);
					}
					else
					{
						new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
						new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
						new_g += qromb_midpnt(charge_ptr, 0.01, .001);
						new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
						new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
						new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
						new_g += qromb_midpnt(charge_ptr, 0.000001, .0000001);
						new_g += qromb_midpnt(charge_ptr, 0.0000001, .00000001);
						new_g += qromb_midpnt(charge_ptr, 0.00000001, xd_global);
					}
				}
				else
				{
					new_g = 0;
				}
			}
			else
			{
				new_g = 0.0;
			}
			if ((use.Get_surface_ptr()->Get_only_counter_ions()) && new_g < 0)
				new_g = 0;
			converge1 = TRUE;
			if (fabs(new_g) >= 1.)
			{
				if (fabs((new_g - charge_ptr->Get_g_map()[z_global].Get_g()) / new_g) > epsilon)
				{
					converge1 = FALSE;
				}
			}
			else
			{
				if (fabs(new_g - charge_ptr->Get_g_map()[z_global].Get_g()) > epsilon)
				{
					converge1 = FALSE;
				}
			}
			if (converge1 == FALSE)
			{
				converge = FALSE;
				if (debug_diffuse_layer == TRUE)
				{
					output_msg(sformatf(
						"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
						(double)z_global,
						(double)charge_ptr->Get_g_map()[z_global].Get_g(),
						(double)new_g,
						(double)(new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
				}
			}
			charge_ptr->Get_g_map()[z_global].Set_g(new_g);
			if (new_g == 0)
			{
				charge_ptr->Get_g_map()[z_global].Set_dg(0.);
			}
			else
			{
				if (charge_ptr->Get_grams() > 0.0)
				{
					LDBLE dg = charge_ptr->Get_grams() *
						charge_ptr->Get_specific_area() * alpha_global *
						g_function(xd_global) / F_C_MOL;
					dg *=
						-2. / (exp(x[j]->master[0]->s->la * LOG_10) *
							exp(x[j]->master[0]->s->la * LOG_10));
					if ((xd_global - 1) < 0.0)
					{
						dg *= -1.0;
					}

					if (fabs(dg) < 1e-8)
					{
						xd1 = exp(-2 * 1e-3 * LOG_10);
						new_g = qromb_midpnt(charge_ptr, 1.0, xd1);
						dg = new_g / .001;
					}
					charge_ptr->Get_g_map()[z_global].Set_dg(dg);
				}
				else
				{
					charge_ptr->Get_g_map()[z_global].Set_dg(0.0);
				}
			}
			temp_g_map[z_global] = charge_ptr->Get_g_map()[z_global];
		}
		if (debug_diffuse_layer == TRUE)
		{
			output_msg(sformatf("\nSurface component %d: charge,\tg,\tdg/dlny,\txd\n",
				(int)charge_ptr->Get_g_map().size()));
			std::map<LDBLE, cxxSurfDL>::iterator it;
			for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
			{
				output_msg(sformatf(
					"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
					(double)it->first,
					(double)it->second.Get_g(),
					(double)it->second.Get_dg(),
					(double)xd_global));
			}
		}
	}
	return (converge);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
g_function(LDBLE x_value)
/* ---------------------------------------------------------------------- */
{
	LDBLE sum, return_value, sum1;
	int i;
	LDBLE ln_x_value;

	if (equal(x_value, 1.0, G_TOL * 100) == TRUE)
		return (0.0);
	sum = 0.0;
	ln_x_value = log(x_value);

	cxxSurfaceCharge *charge_ptr = &(use.Get_surface_ptr()->Get_surface_charges()[0]);
	std::map<LDBLE, cxxSurfDL>::iterator it = charge_ptr->Get_g_map().begin();
	for (; it != charge_ptr->Get_g_map().end(); it++)
	{
		it->second.Set_psi_to_z(exp(ln_x_value * it->first) - 1.0);
	}
	for (i = 0; i < (int)this->s_x.size(); i++)
	{
		if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
		{
			sum += s_x[i]->moles * charge_ptr->Get_g_map()[s_x[i]->z].Get_psi_to_z();
		}
	}
	if (sum < 0.0)
	{
		sum = 0.0;
		sum1 = 0.0;
		output_msg(sformatf(
			"Species\tmoles\tX**z-1\tsum\tsum charge\n"));
		for (i = 0; i < (int)this->s_x.size(); i++)
		{
			if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
			{
				sum += s_x[i]->moles * (pow(x_value, s_x[i]->z) - 1.0);
				sum1 += s_x[i]->moles * s_x[i]->z;
				output_msg(sformatf("%s\t%e\t%e\t%e\t%e\n",
					s_x[i]->name, (double)s_x[i]->moles,
					(double)(pow((LDBLE)x_value, (LDBLE)s_x[i]->z) -
						1.0), (double)sum, (double)sum1));
			}
		}
		error_string = sformatf("Negative sum in g_function, %e\t%e.",
			(double)sum, (double)x_value);
		error_msg(error_string, CONTINUE);
		error_string = sformatf(
			"Solutions must be charge balanced, charge imbalance is %e\n",
			(double)sum1);
		error_msg(error_string, STOP);
	}

	return_value =
		(exp(ln_x_value * z_global) -
			1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
	return (return_value);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
/* ---------------------------------------------------------------------- */
{
	int i, m, ns;
	LDBLE den, dif, dift, ho, hp, w;

	ns = 1;
	dif = fabs(xv - xa[1]);
	/*
	 *   Malloc work space
	 */
	std::vector<double> c, d;
	c.resize((size_t)n + 1);
	d.resize((size_t)n + 1);

	for (i = 1; i <= n; i++)
	{
		dift = fabs(xv - xa[i]);
		if (dift < dif)
		{
			ns = i;
			dif = dift;
		}
		c[i] = ya[i];
		d[i] = ya[i];
	}

	*yv = ya[ns--];
	for (m = 1; m < n; m++)
	{
		for (size_t i = 1; i <= n - m; i++)
		{
			ho = xa[i] - xv;
			hp = xa[i + m] - xv;
			w = c[i + 1] - d[i];
			if ((den = ho - hp) == 0.0)
			{
				error_msg("In subroutine polint.", STOP);
			}
			den = w / den;
			d[i] = hp * den;
			c[i] = ho * den;
		}
		if (2 * ns < (n - m))
		{
			*dy = c[(size_t)ns + 1];
		}
		else
		{
			*dy = d[ns--];
		}
		*yv += *dy;

		/*		*yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
	}
	return;
}

/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
midpnt(LDBLE x1, LDBLE x2, int n)
/* ---------------------------------------------------------------------- */
{
	LDBLE xv, tnm, sum, del, ddel;
	int it, j;

	if (n == 1)
	{
		midpoint_sv = (x2 - x1) * g_function(0.5 * (x1 + x2));
		return (midpoint_sv);
	}
	else
	{
		for (it = 1, j = 1; j < n - 1; j++)
			it *= 3;
		tnm = (LDBLE)it;
		del = (x2 - x1) / (3 * tnm);
		ddel = del + del;
		xv = x1 + 0.5 * del;
		sum = 0.0;
		for (j = 1; j <= it; j++)
		{
#if defined(PHREEQCI_GUI)
			PhreeqcIWait(this);
#endif
			sum += g_function(xv);
			xv += ddel;
			sum += g_function(xv);
			xv += del;
		}
		midpoint_sv = (midpoint_sv + (x2 - x1) * sum / tnm) / 3.0;
		return midpoint_sv;
	}
}

/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
/* ---------------------------------------------------------------------- */
{
	LDBLE ss, dss;
	LDBLE sv[MAX_QUAD + 2], h[MAX_QUAD + 2];
	int j;

	h[0] = 1.0;
	sv[0] = midpnt(x1, x2, 1);
	for (j = 1; j < MAX_QUAD; j++)
	{
		sv[j] = midpnt(x1, x2, j + 1);
		h[j] = h[j - 1] / 9.0;

		if (fabs(sv[j] - sv[j - 1]) <= G_TOL * fabs(sv[j]))
		{
			sv[j] *= charge_ptr->Get_grams() * charge_ptr->Get_specific_area() * alpha_global / F_C_MOL;	/* (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
			if ((x2 - 1) < 0.0)
				sv[j] *= -1.0;
			if (debug_diffuse_layer == TRUE)
			{
				output_msg(sformatf(
					"Iterations in qromb_midpnt: %d\n", j));
			}
			return (sv[j]);
		}

		if (j >= K_POLY - 1)
		{
			polint(&h[j - K_POLY], &sv[j - K_POLY], K_POLY, 0.0, &ss, &dss);
			if (fabs(dss) <= G_TOL * fabs(ss) || fabs(dss) < G_TOL)
			{
				ss *= charge_ptr->Get_grams() * charge_ptr->Get_specific_area() * alpha_global / F_C_MOL;	/* (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
				if ((x2 - 1) < 0.0)
					ss *= -1.0;
				if (debug_diffuse_layer == TRUE)
				{
					output_msg(sformatf(
						"Iterations in qromb_midpnt: %d\n", j));
				}
				return (ss);
			}
		}

	}
	error_string = sformatf(
		"\nToo many iterations integrating diffuse layer.\n");
	error_msg(error_string, STOP);
	return (-999.9);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_init_g(void)
/* ---------------------------------------------------------------------- */
{
	if (use.Get_surface_ptr() == NULL)
		return (OK);
	/*
	 *   calculate g for each surface
	 */
	for (int j = 0; j < count_unknowns; j++)
	{
		if (x[j]->type != SURFACE_CB)
			continue;
		cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
		xd_global = exp(-2 * x[j]->master[0]->s->la * LOG_10);
		/* alpha = 0.02935 @ 25;      (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
		/*  second 1000 is liters/m**3 */
		//alpha_global =	sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
		//	1000.0 * tk_x * 0.5);
		alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
			1000.0 * tk_x * 0.5);

		if (charge_ptr->Get_g_map().size() == 0)
		{
			cxxSurfDL temp_g;
			charge_ptr->Get_g_map()[0.0] = temp_g;
		}
		/*
		 *   calculate g for given surface for each species
		 */
		for (int i = 0; i < (int)this->s_x.size(); i++)
		{
			if (s_x[i]->type > HPLUS)
				continue;

			if (charge_ptr->Get_g_map().find(s_x[i]->z) == charge_ptr->Get_g_map().end())
			{
				cxxSurfDL temp_g;
				/* save g for charge */
				if (charge_ptr->Get_grams() > 0.0)
				{
					temp_g.Set_g(2 * alpha_global * sqrt(mu_x) * (pow(xd_global, s_x[i]->z / 2.0) - 1) *
						charge_ptr->Get_grams() *
						charge_ptr->Get_specific_area() / F_C_MOL);
					temp_g.Set_dg(-s_x[i]->z);
					if (use.Get_surface_ptr()->Get_only_counter_ions() &&
						temp_g.Get_g() < 0)
					{
						temp_g.Set_g(0);
						temp_g.Set_dg(0);
					}
				}
				else
				{
					temp_g.Set_g(0);
					temp_g.Set_dg(-s_x[i]->z);
				}
				charge_ptr->Get_g_map()[s_x[i]->z] = temp_g;
			}

			{
				int is = s_x[i]->number;
				assert(is < (int)s_diff_layer.size());
				// species found in diff_layer
				s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0);
				s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0);
			}
		}
		if (debug_diffuse_layer == TRUE)
		{
			output_msg(sformatf(
				"\nSurface component %d: charge,\tg,\tdg\n",
				(int)charge_ptr->Get_g_map().size()));
			std::map<LDBLE, cxxSurfDL>::iterator it;
			for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
			{
				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\n",
					(double)it->first,
					(double)it->second.Get_g(),
					(double)it->second.Get_dg()));
			}
		}
	}
	return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
initial_surface_water(void)
/* ---------------------------------------------------------------------- */
{
	/*
	 *   In initial surface calculation, need to calculate
	 *   mass of water in diffuse layer.
	 *   diffuse layer water + aqueous solution water = bulk water.
	 *   Ionic strength is fixed, so diffuse-layer water will not change
	 */
	LDBLE debye_length, b, r, rd, ddl_limit, rd_limit, fraction, sum_surfs, l_s;
	LDBLE damp_aq;
	/*
	 *   Debye  length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
	 *
	 *   1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
	 */
	 //debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
	 //	/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
	debye_length = (eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
		/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
	debye_length = sqrt(debye_length);

	/*  ddl is at most the fraction ddl_limit of bulk water */
	ddl_limit = use.Get_surface_ptr()->Get_DDL_limit();

	/*
	 *   Loop through all surface components, calculate each H2O surface (diffuse layer),
	 *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
	 */

	if (use.Get_surface_ptr()->Get_debye_lengths() > 0)
	{
		sum_surfs = 0.0;
		for (int i = 0; i < count_unknowns; i++)
		{
			if (x[i]->type != SURFACE_CB)
				continue;
			cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
			sum_surfs +=
				charge_ptr->Get_specific_area() *
				charge_ptr->Get_grams();
		}
		rd = debye_length * use.Get_surface_ptr()->Get_debye_lengths();
		use.Get_surface_ptr()->Set_thickness(rd);

		if (!sum_surfs)
		{
			for (int i = 0; i < count_unknowns; i++)
			{
				if (x[i]->type != SURFACE_CB)
					continue;
				cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
				charge_ptr->Set_mass_water(0);
			}
		}
		else if (state == INITIAL_SURFACE)
		{
			/* distribute water over DDL (rd) and free pore (r - rd) */
			/* find r: free pore (m3) = pi * (r - rd)^2 * L, where L = A / (2*pi*r),
			   A = sum_surfs = sum of the surface areas */
			b = -2 * (rd + use.Get_solution_ptr()->Get_mass_water() / (1000 * sum_surfs));
			r = 0.5 * (-b + sqrt(b * b - 4 * rd * rd));
			/* DDL (m3) = pi * (r^2 - (r - rd)^2) * L */
			rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
			/* rd should be smaller than r and the ddl limit */
			if (rd > rd_limit)
			{
				mass_water_surfaces_x =
					use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 - ddl_limit);
				r = 0.002 * (mass_water_surfaces_x +
					use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
				rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
				rd = rd_limit;
				use.Get_surface_ptr()->Set_thickness(rd);
			}
			else
				mass_water_surfaces_x =
				(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
			for (int i = 0; i < count_unknowns; i++)
			{
				if (x[i]->type != SURFACE_CB)
					continue;
				cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
				l_s = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
				charge_ptr->Set_mass_water(mass_water_surfaces_x * l_s / sum_surfs);
			}
		}
		else
		{
			r = 0.002 * mass_water_bulk_x / sum_surfs;
			rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
			if (rd > rd_limit)
			{
				rd = rd_limit;
				use.Get_surface_ptr()->Set_thickness(rd);
				fraction = ddl_limit;
			}
			else
				fraction = 1 - pow(r - rd, 2) / (r * r);
			damp_aq = 1.0;
			if (g_iterations > 10)
				damp_aq = 0.2;
			else if (g_iterations > 5)
				damp_aq = 0.5;
			mass_water_surfaces_x = damp_aq * fraction * mass_water_bulk_x +
				(1 - damp_aq) * mass_water_surfaces_x;
			for (int i = 0; i < count_unknowns; i++)
			{
				if (x[i]->type != SURFACE_CB)
					continue;
				cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
				l_s = charge_ptr->Get_specific_area() *
					charge_ptr->Get_grams();
				charge_ptr->Set_mass_water(mass_water_surfaces_x * l_s / sum_surfs);
			}
		}
	}
	else
	{
		/* take constant thickness of, default 1e-8 m (100 Angstroms) */
		mass_water_surfaces_x = 0.0;
		for (int i = 0; i < count_unknowns; i++)
		{
			if (x[i]->type != SURFACE_CB)
				continue;
			cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
			charge_ptr->Set_mass_water(charge_ptr->Get_specific_area() *
				charge_ptr->Get_grams() * use.Get_surface_ptr()->Get_thickness() *
				1000);
			mass_water_surfaces_x += charge_ptr->Get_mass_water();
		}
	}

	if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC)
		mass_water_bulk_x = mass_water_aq_x + mass_water_surfaces_x;
	else
	{
		/*  for variable distribution of water over DDL and bulk... */
		if (state > INITIAL_SURFACE)
			mass_water_aq_x = mass_water_bulk_x - mass_water_surfaces_x;
		else
			mass_water_bulk_x = mass_water_aq_x + mass_water_surfaces_x;
	}

	return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
/* ---------------------------------------------------------------------- */
{
	LDBLE mass_water_surface;
	LDBLE molality, moles_excess, moles_surface;

	if (use.Get_surface_ptr() == NULL)
		return (OK);
	/*
	 *   Find position of component in list of components
	 */

	 /*
	  *   Loop through all surface components, calculate each H2O surface (diffuse layer),
	  *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
	  */
	count_elts = 0;
	paren_count = 0;
	mass_water_surface = charge_ptr->Get_mass_water();
	for (int j = 0; j < (int)this->s_x.size(); j++)
	{
		if (s_x[j]->type > HPLUS)
			continue;
		molality = under(s_x[j]->lm);
		LDBLE g = charge_ptr->Get_g_map()[s_x[j]->z].Get_g();
		if (s_x[j]->erm_ddl != 1)
		{
			LDBLE ratio_aq = mass_water_surface / mass_water_aq_x;
			LDBLE g2 = g / ratio_aq + 1;
			g = ratio_aq * (g2 * s_x[j]->erm_ddl - 1);
		}
		moles_excess = mass_water_aq_x * molality * g;
		moles_surface = mass_water_surface * molality + moles_excess;
		/*
		 *   Accumulate elements in diffuse layer
		 */
		add_elt_list(s_x[j]->next_elt, moles_surface);
	}
	add_elt_list(s_h2o->next_elt, mass_water_surface / gfw_water);
	elt_list_combine();
	return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_all_donnan(void)
/* ---------------------------------------------------------------------- */
{
	bool converge;
	int cd_m;
	LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_surf_aq, co_ion;
	LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1, viscos;
	cxxSurface *surf_ptr = use.Get_surface_ptr();

	if (surf_ptr == NULL)
		return (OK);
	f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
		tk_x * mu_x);
	bool only_count = surf_ptr->Get_only_counter_ions();
	bool correct_D = surf_ptr->Get_correct_D();
	/*   calculate g for each surface...
	 */
	if (!calculating_deriv || surf_ptr->Get_debye_lengths() ||
		correct_D) // DL_pitz && correct_D
		initial_surface_water();
	// z1, z2, fr_cat2 are the counter-ions, z_1, z_2, fr_ani2 are for co-ions.
	LDBLE nDbl = 1, db_lim = 2, f_free, fr_cat2, fr_ani2;
	LDBLE z1, z2, z_1, z_2;
	z1 = z2 = z_1 = z_2 = f_free = fr_cat2 = fr_ani2 = 0;
	/*
	 *  sum eq of each charge number in solution...
	 */
	std::map<LDBLE, LDBLE>::iterator it;
	for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
	{
		it->second = 0.0;
	}
	charge_group_map.clear();
	for (int i = 0; i < (int)this->s_x.size(); i++)
	{
		if (s_x[i]->type > HPLUS)
			continue;
		charge_group_map[s_x[i]->z] += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
	}
	for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
	{
		if      (it->first < -1.5) { z_2 += it->second; continue; }
		else if (it->first < 0)    { z_1 += it->second; continue; }
		else if (it->first < 1.5)  { z1  += it->second; continue; }
		else					   { z2  += it->second; continue; }
	}
	if (correct_D)
	{
		if ((nDbl = surf_ptr->Get_debye_lengths()) == 0)
		{
			LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3);
			nDbl = surf_ptr->Get_thickness() / debye_length;
		}
		fr_ani2 = z_2 / (z_1 + z_2);
		fr_cat2 = z2  / (z1 + z2);
		db_lim = 2 - 0.5 * (fr_cat2 + fr_ani2);
		if (nDbl > db_lim)
		{
			f_free = 1 - db_lim / nDbl;
			if (f_free < 0) f_free = 0;
		}
	}

	converge = TRUE;
	viscos = 0;
	for (int j = 0; j < count_unknowns; j++)
	{
		if (x[j]->type != SURFACE_CB)
			continue;
		cxxSurfaceCharge *charge_ptr = surf_ptr->Find_charge(x[j]->surface_charge);

		if (debug_diffuse_layer == TRUE)
			output_msg(sformatf("Calc_all_g, X[%d]\n", j));

		/* find surface charge from potential... */
		A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
		if (surf_ptr->Get_type() == cxxSurface::CD_MUSIC)
		{
			f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10;	/* -FPsi/RT */
			f_psi = f_psi / 2;
			cd_m = 1;
		}
		else
		{
			f_psi = x[j]->master[0]->s->la * LOG_10;
			cd_m = -1;
		}
		surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
		LDBLE lim_seq = 5e3;
		if (correct_D) lim_seq = 5e3;
		if (fabs(surf_chrg_eq) > lim_seq)
		{
			surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq);
			var1 = surf_chrg_eq / (A_surf * f_sinh / F_C_MOL);
			var1 = (var1 + sqrt(var1 * var1 + 1));
			f_psi = (var1 > 1e-8 ? log(var1) : -18.4);
			surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
			x[j]->master[0]->s->la = f_psi / LOG_10;
		}
		/* also for the derivative... */
		dif = 1e-5;
		f_psi2 = f_psi + dif;
		surf_chrg_eq2 = A_surf * f_sinh * sinh(f_psi2) / F_C_MOL;

		/* find psi_avg that matches surface charge... */
		std::vector<LDBLE> zcorr(charge_group_map.size());
		std::vector<LDBLE> zcorr2(charge_group_map.size());
		//LDBLE fD = 0;
		psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, f_free, zcorr);
		psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, f_free, zcorr2);

		/*output_msg(sformatf( "psi's  %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */

		/* fill in g's */
		ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
		ratio_surf_aq = charge_ptr->Get_mass_water() / mass_water_surfaces_x;
		//ratio_surf_aq = charge_ptr->Get_mass_water() / mass_water_bulk_x;
		if (correct_D)
			ratio_aq *= (1 - f_free);
		int z_iter = 0;
		for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
		{
			LDBLE z = it->first, z1 = z;
			co_ion = surf_chrg_eq * z;
			if (correct_D)
				z1 = zcorr[z_iter];

			if (!ratio_aq)
			{
				charge_ptr->Get_g_map()[z].Set_g(0);
				charge_ptr->Get_g_map()[z].Set_dg(0);
				charge_ptr->Get_z_gMCD_map()[z] = 0;
				converge = true;
				continue;
			}
			new_g = ratio_aq * (exp(cd_m * z1 * psi_avg) - 1);
			if (only_count && co_ion > 0)
				new_g = -ratio_aq;
			if (new_g <= -ratio_aq)
				new_g = -ratio_aq + G_TOL * 1e-3;
			new_g2 = ratio_aq * (exp(cd_m * z1 * psi_avg2) - 1);
			if (only_count && co_ion > 0)
				new_g2 = -ratio_aq;
			if (new_g2 <= -ratio_aq)
				new_g2 = -ratio_aq + G_TOL * 1e-3;
			if (fabs(new_g) >= 1)
			{
				if (fabs((new_g - charge_ptr->Get_g_map()[z].Get_g()) / new_g) > convergence_tolerance)
				{
					converge = FALSE;
				}
			}
			else
			{
				if (fabs(new_g - charge_ptr->Get_g_map()[z].Get_g()) > convergence_tolerance)
				{
					converge = FALSE;
				}
			}
			charge_ptr->Get_g_map()[z].Set_g(new_g);
			if (new_g != 0)
			{
				charge_ptr->Get_g_map()[z].Set_dg((new_g2 - new_g) / dif);
			}
			else
			{
				charge_ptr->Get_g_map()[z].Set_dg(-z);
			}
			/* save Boltzmann factor * water fraction for MCD calc's in transport */
			if (converge)
			{
				if (only_count && co_ion > 0) // co-ions are not in the DL
						charge_ptr->Get_z_gMCD_map()[z] = 0;
				else
				{
					charge_ptr->Get_z_gMCD_map()[z] = (exp(cd_m * z1 * psi_avg) * (1 - f_free) + f_free) *
						ratio_surf_aq;// * s_x[]->moles == mol_DL in charge_ptr
				}
			}
			z_iter++;
		}

		charge_ptr->Set_f_free(f_free);
		if (debug_diffuse_layer == TRUE)
		{
			std::string name = x[j]->master[0]->elt->name;
			Utilities::replace("_psi", "", name);
			output_msg(sformatf(
				"\nDonnan all on %s (%d): charge, \tg, \tdg, \tzcorr, \tPsi_surface = %8f V, \tDebye lengths = %8f. \n",
				name.c_str(), (int)charge_ptr->Get_g_map().size(),
				x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ,
				nDbl));
			int i = 0;
			for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
				i_it != charge_ptr->Get_g_map().end(); i_it++)
			{
				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
					(double)i_it->first,
					(double)i_it->second.Get_g(),
					(double)i_it->second.Get_dg(),
					(double)zcorr[i]));
				i++;
			}
		}
	}
	return (converge);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_init_donnan(void)
/* ---------------------------------------------------------------------- */
{
	LDBLE f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq;
	cxxSurface *surf_ptr = use.Get_surface_ptr();

	if (surf_ptr == NULL)
		return (OK);
	f_sinh =
		//sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
		//	 tk_x * mu_x);
		sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
			tk_x * mu_x);
	if (convergence_tolerance >= 1e-8)
	{
		G_TOL = 1e-9;
	}
	else
	{
		G_TOL = 1e-13;
	}
	/*
	 *  sum eq of each charge number in solution...
	 */
	charge_group_map.clear();
	charge_group_map[0.0] = 0.0;

	for (int i = 0; i < (int)this->s_x.size(); i++)
	{
		if (s_x[i]->type > HPLUS)
			continue;
		if (charge_group_map.find(s_x[i]->z) != charge_group_map.end())
		{
			charge_group_map.find(s_x[i]->z)->second += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
		}
		else
		{
			charge_group_map[s_x[i]->z] = s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
		}
	}
	std::vector<LDBLE> zcorr(charge_group_map.size());
	/*
	 *   calculate g for each surface...
	 */
	for (int j = 0; j < count_unknowns; j++)
	{
		if (x[j]->type != SURFACE_CB)
			continue;
		cxxSurfaceCharge *charge_ptr = surf_ptr->Find_charge(x[j]->surface_charge);
		charge_ptr->Get_g_map().clear();

		/* find surface charge from potential... */
		A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
		if (surf_ptr->Get_type() == cxxSurface::CD_MUSIC)
		{
			f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10;	/* -FPsi/RT */
			f_psi = f_psi / 2;
		}
		else
			f_psi = x[j]->master[0]->s->la * LOG_10;
		surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;

		/* find psi_avg that matches surface charge... */
		psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, 0, zcorr);

		/* fill in g's */
		ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;

		std::map<LDBLE, LDBLE>::iterator kit;
		for (kit = charge_group_map.begin(); kit != charge_group_map.end(); kit++)
		{
			LDBLE z = kit->first;
			LDBLE eq = kit->second;

			charge_ptr->Get_g_map()[z].Set_g(ratio_aq * (exp(-z * psi_avg) - 1));

			if (surf_ptr->Get_only_counter_ions()
				&& ((surf_chrg_eq < 0 && z < 0)
					|| (surf_chrg_eq > 0 && z > 0)))
				charge_ptr->Get_g_map()[z].Set_g(-ratio_aq);

			if (charge_ptr->Get_g_map()[z].Get_g() != 0)
			{
				charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) /
					(eq * F_C_MOL));
			}
			else
			{
				charge_ptr->Get_g_map()[z].Set_dg(-z);
			}
			/* save g for species */
			for (int i = 0; i < (int)this->s_x.size(); i++)
			{
				int is = s_x[i]->number;
				assert(is < (int)s_diff_layer.size());

				s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0.0);
				s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0.0);
			}
		}
		if (debug_diffuse_layer == TRUE)
		{
			std::string name = x[j]->master[0]->elt->name;
			Utilities::replace("_psi", "", name);
			output_msg(sformatf(
				"\nDonnan init on %s : charge, \tg, \tdg, \tzcorr, Psi_surface = %8f V. \n",
				name.c_str(),
				x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
				tk_x / F_KJ_V_EQ));
			for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
				i_it != charge_ptr->Get_g_map().end(); i_it++)
			{
				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
					(double)i_it->first,
					(double)i_it->second.Get_g(),
					(double)i_it->second.Get_dg()));
			}
		}
	}
	return (OK);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, LDBLE f_free, std::vector<LDBLE> &zcorr)
/* ---------------------------------------------------------------------- */
{
	/*
	 * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
	 */
	LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, Z1, Z1_c, eq, co_ion, sum_counter, sum_co;
	LDBLE z1, z2, z_1, z_2;
	ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
	p = 0;
	if (surf_chrg_eq == 0 || ratio_aq == 0)
		return (0.0);
	else if (surf_chrg_eq < 0)
		p = -0.5 * log(-surf_chrg_eq * ratio_aq * (1 - f_free) / mu_x + 1);
	else if (surf_chrg_eq > 0)
		p = 0.5 * log(surf_chrg_eq * ratio_aq * (1 - f_free) / mu_x + 1);
	/*
	 * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
	 *  g(p) = exp(-p * z_i) * ratio_aq * (1 - f_free)
	 * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
	 *		      g(p) = (exp(-p *z_i) - 1) * ratio_aq
	 * with correct_D true and f_free > 0:
	     c_edl = c_free * (f_free + (1 - f_free) * exp(-p * z_i))
	 * with correct_D true and f_free == 0:
	 * correct ions to better match the integrated PB concentrations:
		 Gamma = abs(surf_chrg_eq / A_surf / 1e-6)
		 a = cgc[1] * nDbl**cgc[2]
		 b = Gamma**cgc[3] / abs(log10(I))
		 counter_ions...
		 z == 1? z1 = cgc[0] * I**(a * b)
		 z == 2? z2 = 2 * cgc[4] * I**(cgc[5] * a * b)
		 co_ions...
		 c = cgc[7] * nDbl**cgc[8] * Gamma**cgc[9]
		 z == -1? z_1 = -cgc[6] * I**(c)
		 z == -2? z_2 = -2 * cgc[10] * I**(c * cgc[11])
		 c_edl = c_free * exp(-p * z_i)
	 */

	cxxSurface *surf_ptr = use.Get_surface_ptr();
	bool correct_D = surf_ptr->Get_correct_D(),  local_correct_D = correct_D;
	bool only_count = surf_ptr->Get_only_counter_ions();
	LDBLE Gamma, cgc[12] = { 0.3805, -0.0106, 1.96, 0.812,
							0.395, 2.13,
							0.380, 0.0408, 0.799, 0.594,
							0.373, 1.181 };
	if (correct_D)
	{
		if (f_free)
		{
			z1 = z2 = z_1 = z_2 = 1;
		}
		else
		{
			if (!surf_ptr->Donnan_factors.empty())
			{
				std::copy(surf_ptr->Donnan_factors.begin(), surf_ptr->Donnan_factors.end(), cgc);
				z1 = cgc[0];
				z2 = cgc[1];
				z_1 = cgc[2];
				z_2 = cgc[3];
			}
			else
			{
				Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6;
				LDBLE a = cgc[1] * pow(nDbl, cgc[2]),
					b = pow(Gamma, cgc[3]) / abs(log10(mu_x));
				// counter_ions...
				z1 = cgc[0] * pow(mu_x, (a * b));
				z2 = cgc[4] * pow(mu_x, (cgc[5] * a * b));
				if (z1 > 1) z1 = 1;
				if (z2 > 1) z2 = 1;
				// co_ions...
				LDBLE c = cgc[7] * pow(nDbl, cgc[8]) * pow(Gamma, cgc[9]);
				z_1 = cgc[6] * pow(mu_x, c);
				z_2 = cgc[10] * pow(mu_x, (c * cgc[11]));
			}
		}
	}

	int l_iter = 0, z_iter;
	sum_co = sum_counter = 0;
	do
	{
		fd = surf_chrg_eq;
		fd1 = 0.0;
		z_iter = 0;
		if (l_iter == 1 && local_correct_D && fabs(sum_counter) < fabs(sum_co))
		{
			local_correct_D = false;
			l_iter = 0;
		}
		std::map<LDBLE, LDBLE>::iterator it;
		for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
		{
			z = it->first;
			Z1 = z;
			if (l_iter == 0) zcorr[z_iter] = z;
			co_ion = surf_chrg_eq * z;
			if (!z || (only_count && co_ion > 0))
			{
				z_iter++;
				continue;
			}
			if (nDbl && local_correct_D)
			{
				/*psi_DL = fabs(p * p_psi);*/
				if (co_ion < 0)
				{//counter-ion
					if (fabs(z) > 1.5) temp = z2;
					else               temp = z1;
					sum_counter += z * temp;
				}
				else
				{// co-ion
					if (fabs(z) > 1.5) temp = z_2;
					else               temp = z_1;
					sum_co += z * temp;
				}
				zcorr[z_iter] = z * temp;
			}
			Z1 = zcorr[z_iter];
			eq = it->second;
			temp = exp(-Z1 * p) * ratio_aq * (1 - f_free);

			fd += eq * temp;
			fd1 -= Z1 * eq * temp;
			if (z == 1) Z1_c = Z1;
			z_iter++;
		}
		fd /= -fd1;
		p += (fd > 1) ? 1 : ((fd < -1) ? -1 : fd);
		if (fabs(p) < G_TOL)
			p = 0.0;
		l_iter++;
		if (l_iter > 50)
		{
			pr.all = TRUE;
			pr.exchange = TRUE;
			pr.headings = TRUE;
			pr.pp_assemblage = TRUE;
			pr.species = TRUE;
			pr.surface = TRUE;
			pr.totals = TRUE;
			print_all();
			error_string = sformatf(
				"\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n",
				(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
			error_msg(error_string, STOP);
		}
	} while (fabs(fd) > 1e-12 && p != 0.0);
	if (debug_diffuse_layer == TRUE)
		output_msg(sformatf(
			"iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n",
			l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * Z1_c * p_psi)));

	return (p);
}