/Software/PhreeqC/modelengine.cpp

Large files files are truncated, but you can click here to view the full file

#include <math.h>

#include <float.h>

#include "modelengine.h"

#include "conio.h"



//===========================

// Constructors & Destructors

//===========================

ModelEngine::ModelEngine(GlobalData *gd, ModelData *md):

SSAssemblageModel(gd, md)

{

	this->gd = gd;

	this->md = md;



	SetData(md);



	r_temp = NULL;

	r_temp = new Reaction;



	count_unknowns = 0;



	error = NULL;

	error_string[0] = '\0';

}



ModelEngine::~ModelEngine()

{

	delete r_temp;

}



//===========================

// Public functions

//===========================

bool ModelEngine::ConvertUnits(Solution *sol_p)

{

  int i;

  LDBLE sum_solutes;

  Master *m;

  Conc *conc;

	

  sum_solutes = exp(-sol_p->ph * md->LOG_10);



  for (i = 0; i < sol_p->totals->Count(); i++)

	{

    conc = (*sol_p->totals)[i];



    conc->moles = 0.0;

    if (conc->name == "H(1)" || conc->name == "E")

      continue;



		if (conc->input_conc <= 0)

      continue;



    if (conc->gfw <= 0.0)

    {

			if (!conc->as.IsEmpty())

      {

				if (!ComputeGFW(conc->as, conc->gfw))

					return false;



				if (conc->name == "Alkalinity" && conc->as == "CaCO3")

					conc->gfw /= 2.0;

      }

      else

      {

				if ((m = gd->master_list.Search(&conc->name(0, " "), true)) != NULL)

					conc->gfw = m->gfw;

				else

					return false;

      }

    }



    conc->moles = conc->input_conc;



    if (Check_l_(sol_p->units))

      conc->moles *= (LDBLE)1.0 / (sol_p->density);



		conc->moles *= ConversionFactorForMoles(conc->units);



    if (Check_g_kgs_(conc->units) || Check_g_l_(conc->units))

      sum_solutes += conc->moles;

    else if (Check_mol_kgs_(conc->units) || Check_mol_l_(conc->units) || Check_eq_l_(conc->units))

      sum_solutes += (conc->moles * conc->gfw);



		if (Check_g_(conc->units) && conc->gfw != 0.0)

      conc->moles /= conc->gfw;

  }



	if (Check_kgs_(sol_p->units) || Check_l_(sol_p->units))

  {

    md->mass_water_aq_x = (LDBLE)1.0 - ((LDBLE)1e-3 * sum_solutes);

    for (i = 0; i < sol_p->totals->Count(); i++)

      (*sol_p->totals)[i]->moles /= md->mass_water_aq_x;

  }



  md->mass_water_aq_x = sol_p->mass_water;

  for (i = 0; i < sol_p->totals->Count(); i++)

    (*sol_p->totals)[i]->moles *= md->mass_water_aq_x;



	return true;

}





bool ModelEngine::Clear()

{

	int i;

	Solution *sol_p = md->use.sol_p;



	// Clear species solution-dependent data

	for (i = 0; i < gd->species_list.Count(); i++)

	{

    gd->species_list[i]->in = false;

		//gd->species_list[i]->rewrite = false;

	}



	// Set pe structure

	sol_p->pe->CopyTo(md->pe_x);

  md->default_pe_x = sol_p->default_pe;



	// Clear master species solution-dependent data

	Master *m;

	PEData *ped;

	for (i = 0; i < gd->master_list.Count(); i++)

  {

		m = gd->master_list[i];

    

		m->in = false;

		m->rewrite = false;

		m->u = NULL;



		ped = (*md->pe_x)[sol_p->default_pe];

		m->pe_rxn = &(ped->rxn);



		// copy primary reaction to secondary reaction

		m->rxn_primary->CopyTo(m->rxn_secondary);

  }



  if (md->state == INITIAL_SOLUTION)

  {

    gd->s_h2o->secondary->in = true;

		gd->s_h2o->secondary->rewrite = false;

    gd->s_hplus->secondary->in = true;

		gd->s_hplus->secondary->rewrite = false;

  }

  else

  {

    gd->s_h2o->primary->in = true;

		gd->s_h2o->primary->rewrite = false;

    gd->s_hplus->primary->in = true;

		gd->s_hplus->primary->rewrite = false;

  }

  gd->s_eminus->primary->in = true;

	gd->s_eminus->primary->rewrite = false;



	// Set all unknown pointers to NULL

  md->mb_unknown = NULL;

  md->ah2o_unknown = NULL;

  md->mass_hydrogen_unknown = NULL;

  md->mass_oxygen_unknown = NULL;

  md->mu_unknown = NULL;

  md->alkalinity_unknown = NULL;

  md->carbon_unknown = NULL;

  md->ph_unknown = NULL;

  md->pe_unknown = NULL;

  md->charge_balance_unknown = NULL;

  md->solution_phase_boundary_unknown = NULL;

  md->pure_phase_unknown = NULL;

  md->exchange_unknown = NULL;

  md->surface_unknown = NULL;

  md->gas_unknown = NULL;

  md->s_s_unknown = NULL;

 

	// Free arrays used in model   

  FreeModelAllocs(); //free_model_allocs ();



  return true;

}





//===========================

// Private functions

//===========================

bool ModelEngine::GetMasterList(String &list, Master *m, ListOfPointers<Master> *copy)

{

	String *token;	

	Master *m_ptr, *m_1;

	int i;



	tl.SetString(list);

	copy->Clear();



	// Make list of master species pointers

	if (m == m->s->primary)

	{

		// First in list is primary species

    for (i = 0; i < gd->master_list.Count(); i++)

		{

			m_1 = gd->master_list[i];

			if (m_1 == m)

				break;

		}

    i++;		

		

		if (i >= gd->master_list.Count())

		{

			// Element has only one valence

			copy->AddNew(m);

		}

		else

		{

			m_1 = gd->master_list[i];



			if (m != m_1->e->primary)

			{

				// Element has only one valence

				copy->AddNew(m);

			}

			else

			{

				// Element has multiple valences

				if (m->s->secondary == NULL)

					return false;



				copy->AddNew(m->s->secondary);



				while (i < gd->master_list.Count() && (m_ptr = gd->master_list[i])->e->primary == m)

				{

					if (m_ptr->s->primary == NULL)

						copy->AddNew(m_ptr);



					i++;

				}

			}

		}

	}

	else

	{

		// First in list is secondary species, Include all valences from input

    copy->AddNew(m);



		while (!tl.EOTL())

    {

			token = tl.NextToken();

			m_ptr = gd->master_list.Search(token, true);



      if (m_ptr != NULL)

				copy->AddNew(m_ptr);

    }

  }



	return true;

}



bool ModelEngine::RewriteMasterToSecondary(Master *m1, Master *m2)

{

	//

	// Write equation for secondary master species in terms of another secondary master species

	// Store result in rxn_secondary of master_ptr.

	//



  LDBLE coef1, coef2;

  Master *m_p1, *m_p2;



	// Check that the two master species have the same primary master species

	m_p1 = m1->e->primary;

  m_p2 = m2->e->primary;



  if (m_p1 != m_p2 || m_p1 == NULL)

		return false;



	// Find coefficient of primary master in reaction

	if (!GetReactionCoef(m1->rxn_primary, m_p1->s->name, coef1) || !GetReactionCoef(m2->rxn_primary, m_p1->s->name, coef2))

		return false;



  if (Equal(coef1, (LDBLE)0.0, (LDBLE)TOLERANCE) || Equal(coef2, (LDBLE)0.0, (LDBLE)TOLERANCE))

		return false;



	// Rewrite equation to secondary master species

	r_temp->AddTRXN(m1->rxn_primary, 1.0, false);

  r_temp->AddTRXN(m2->rxn_primary, -coef1 / coef2, true);



  return false;

}



bool ModelEngine::GetReactionCoef(Reaction *r, String name, LDBLE &coef, int start)

{

	//

	// Finds coefficient of token in reaction.

	// input: r, pointer to a reaction structure

	//        name, string to find as reaction token

	//        start, search start position in r (default = 1)

	//

	// Return: 0.0, if token not found

	//         coefficient of token, if found.

	//



  coef = 0.0;



	ReactionToken *r_token;	

	for (int i = start; i < r->token_list->Count(); i++)

  {

		r_token = r->token_list->Element(i);



    if (r_token->s->name == name)

    {

      coef = r_token->coef;

      break;

    }

  }



	return true;

}

bool ModelEngine::CheckIn()

{

	//

	// Routine goes through trxn to determine if each master species is in this model.

  // Assumes equation is written in terms of primary and secondary species

  // Checks to see if in is TRUE or REWRITE for each species

  // Returns TRUE if in model FALSE if not

	//



  ReactionToken *t;



	for (int i = 1; i < r_temp->token_list->Count(); i++)

  {

		t = r_temp->token_list->Element(i);



		// Check primary master species in

    if (t->s->primary != NULL && t->s->primary->in)

      continue;



		// Check secondary master species

    if (t->s->secondary != NULL	&& (t->s->secondary->in || t->s->secondary->rewrite))

      continue;



		// Must be primary master species that is out

		return false;

  }



  return true;

}



bool ModelEngine::WriteMassActionEqnX()

{

	//

	// Reduce mass-action equation to the master species that are in the model

	//



  LDBLE coef_e;

  int count;

  int i, count_rxn_orig;





	// Rewrite any secondary master species flagged REWRITE. Replace pe if necessary

	count = 0;

  bool repeat = true;

  while (repeat)

  {

    count++;

    if (count > MAX_ADD_EQUATIONS)

			return false;



    repeat = false;



		ReactionToken *r;

		count_rxn_orig = r_temp->token_list->Count();

    for (i = 1; i < count_rxn_orig; i++)

    {

			r = r_temp->token_list->Element(i);

      if (r->s->secondary == NULL)

				continue;



      if (r->s->secondary->rewrite)

      {

				repeat = true;



				if (!GetReactionCoef(r->s->secondary->rxn_secondary, "e-", coef_e))

					return false;



				r_temp->AddTRXN(r->s->secondary->rxn_secondary,	r->coef, false);



				if (!Equal(coef_e, (LDBLE)0.0, (LDBLE)TOLERANCE))

					r_temp->AddTRXN(*r->s->secondary->pe_rxn, r->coef * coef_e, false);

      }

    }

    

		r_temp->TRXNCombine();

  }



  return true;

}



bool ModelEngine::AddPotentialFactor()

{

	//

	// Add the potential factor to surface mass-action equations.

	// Factor is essentially the activity coefficient, representing

	// the work required to bring charged ions to the surface

	//



  int i;

  LDBLE sum_z;

  Master *master_ptr;

  Unknown *unknown_ptr;



  if (md->use.sur_p->type != DDL)

    return true;



  sum_z = 0.0;

  master_ptr = NULL;



	// Find sum of charge of aqueous species and surface master species

	ReactionToken *r;

	for (i = 1; i < r_temp->token_list->Count(); i++)

  {

		r = r_temp->token_list->Element(i);



    if (r->s->type == AQ || r->s == gd->s_hplus || r->s == gd->s_eminus)

      sum_z += r->s->z * r->coef;



		if (r->s->type == SURF)

      master_ptr = r->s->primary;

  }



	if (master_ptr == NULL)

		return false;

		

	// Find potential unknown for surface species

  unknown_ptr = FindSurfaceChargeUnknown(master_ptr->name, SURF_PSI);



  if (unknown_ptr == NULL)

		return false;



	master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include psi in mass action equation

  if (master_ptr != NULL)

  {

		r = r_temp->token_list->AddNew();



    r->name = master_ptr->s->name;

    r->s = master_ptr->s;

    r->coef = (LDBLE)-2.0 * sum_z;

  }

  else

		return false;



	return true;

}



bool ModelEngine::AddCDMusicFactors(int n)

{

	//

	// Add the potential factors for cd_music to surface mass-action equations.

	// Factors are essentially the activity coefficient, representing

	// the work required to bring charged ions to the three charge layers 

	// of the cd_music model

	//



  int i;

	String name;

  Master *master_ptr;

  Unknown *unknown_ptr;

	ReactionToken *r;



  if (md->use.sur_p->type != CD_MUSIC)

    return true;



  master_ptr = NULL;



	// Find sum of charge of aqueous species and surface master species

	for (i = 1; i < r_temp->token_list->Count(); i++)

	{

		r = r_temp->token_list->Element(i);



    if (r->s->type == SURF)

      master_ptr = r->s->primary;

	}



	if (master_ptr == NULL)

		return false;



	name = master_ptr->e->name;



	// Find potential unknown for surface species

	unknown_ptr = FindSurfaceChargeUnknown (name, SURF_PSI);



  if (unknown_ptr == NULL)

		return false;



	master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include psi in mass action equation

	r = r_temp->token_list->AddNew();



  r->name = master_ptr->s->name;

  r->s = master_ptr->s;

	r->coef = gd->species_list[n]->dz[0];

  

	// Plane 1

  unknown_ptr = FindSurfaceChargeUnknown(name, SURF_PSI1);



  if (unknown_ptr == NULL)

		return false;



  master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include psi in mass action equation

	r = r_temp->token_list->AddNew();



  r->name = master_ptr->s->name;

  r->s = master_ptr->s;

	r->coef = gd->species_list[n]->dz[1];

  

	// Plane 2

  unknown_ptr = FindSurfaceChargeUnknown(name, SURF_PSI2);



  if (unknown_ptr == NULL)

		return false;



  master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include psi in mass action equation

	r = r_temp->token_list->AddNew();



  r->name = master_ptr->s->name;

  r->s = master_ptr->s;

	r->coef = gd->species_list[n]->dz[2];



	return true;

}



Unknown *ModelEngine::FindSurfaceChargeUnknown(String &name, int plane)

{

  int i;

  String token;



	token = name;

	token.Replace ("_", " ");

	token = token(0, " ");



  if (plane == SURF_PSI)

		token += "_CB";

  else if (plane == SURF_PSI1)

		token += "_CBb";

  else if (plane == SURF_PSI2)

		token += "_CBd";



	for (i = 0; i < count_unknowns; i++)

		if (token == (*unknown_list)[i]->name)

      return (*unknown_list)[i];



  return NULL;

}

bool ModelEngine::WriteMBEqnX()

{

	//

	// Rewrite any secondary master species flagged REWRITE

	// Don't add in any pe reactions

	//



  int count;

	bool repeat;

  int i, count_rxn_orig;

  int j, k;

	char *ptr, token[DEFAULT_STRING_LENGTH];

  Master *master_ptr;

	ReactionToken *rxn_token;



  count = 0;

  repeat = true;

  while (repeat)

  {

    count++;

    if (count > MAX_ADD_EQUATIONS)

			return false;



		repeat = false;



		count_rxn_orig = r_temp->token_list->Count();

    for (i = 1; i < count_rxn_orig; i++)

    {

			rxn_token = r_temp->token_list->Element(i);



      if (rxn_token->s->secondary == NULL)

				continue;



			if (rxn_token->s->secondary->rewrite)

      {

				repeat = true;

				r_temp->AddTRXN(rxn_token->s->secondary->rxn_secondary, rxn_token->coef, false);

      }

    }



		r_temp->TRXNCombine();

  }



	eos_list->Clear();

  parent_count = 0;



	ReactionToken *r;

	ElementOfSpecies *eos_p;

  for (i = 1; i < r_temp->token_list->Count(); i++)

  {

		r = r_temp->token_list->Element(i);

    j = eos_list->Count();



		

		r->s->name.Copy(token);

		ptr = &token[0];

		if (!GetElementsInSpecies(&ptr, r->coef))

			return false;



		for (k = j; k < eos_list->Count(); k++)

    {

			eos_p = (*eos_list)[k];

      if (r->s->secondary != NULL)

				master_ptr = r->s->secondary->e->primary;   

      else

      	master_ptr = r->s->primary;

      

      if (eos_p->e == master_ptr->e)

      {

				eos_p->coef = 0.0;

				break;

      }

    }



    if (r->s->secondary == NULL)

			r->s->primary->e->name.Copy(token);

    else

			r->s->secondary->e->name.Copy(token);



		ptr = &token[0];

		if (!GetSecondaryElementsInSpecies(&ptr, r->coef))

			return false;

  }



	CombineElements();



  return true;

}



bool ModelEngine::AddSurfaceChargeBalance()

{

	//

	// Include charge balance in list for mass-balance equations

	//



  int i;

  char *ptr;

	char token[DEFAULT_STRING_LENGTH];



  Master *master_ptr;

  Unknown *unknown_ptr;

	ElementOfSpecies *eos_p;



  if (md->use.sur_p->type != DDL)

    return true;



  master_ptr = NULL;



	// Find master species

	for (i = 0; i < eos_list->Count(); i++)

	{

		eos_p = (*eos_list)[i];



    if (eos_p->e->primary->s->type == SURF)

    {

      master_ptr = eos_p->e->primary;

      break;

    }

	}



	if (i >= eos_list->Count())

		throw exception();



	// Find potential unknown for surface species

	unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI);



  if (unknown_ptr == NULL)

		throw exception();



	master_ptr = (*unknown_ptr->master)[0];	// potential for surface component 



	master_ptr->e->name.Copy(token);

	ptr = &token[0];

  

	// Include charge balance in list for mass-balance equations

	LDBLE coef = 1.0;

	if (!GetSecondaryElementsInSpecies(&ptr, coef))

		return false;



  return true;

}



bool ModelEngine::AddCDMusicChargeBalances(int n)

{

	//

	// Add the potential factor to surface mass-action equations.

	// Factor is essentially the activity coefficient, representing

	// the work required to bring charged ions to the surface

	//



  int i;

  char *ptr;

	char token[DEFAULT_STRING_LENGTH];



  Master *master_ptr;

  Unknown *unknown_ptr;

	ElementOfSpecies *eos_p;



  if (md->use.sur_p->type != CD_MUSIC)

    return true;



  master_ptr = NULL;



	// Find master species

	for (i = 0; i < eos_list->Count(); i++)

	{

		eos_p = (*eos_list)[i];



		if (eos_p->e->primary->s->type == SURF)

    {

      master_ptr = eos_p->e->primary;

      break;

    }

	}

  

  if (i >= eos_list->Count())

		throw exception();



	// Find potential unknown for plane 0

  unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI);

  master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include charge balance in list for mass-balance equations

	master_ptr->e->name.Copy(token);

	ptr = &token[0];

	if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[0]))

		return false;



	// Find potential unknown for plane 1

  unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI1);

  master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include charge balance in list for mass-balance equations

	master_ptr->e->name.Copy(token);

	ptr = &token[0];

	if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[1]))

		return false;



	// Find potential unknown for plane 2

  unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI2);

  master_ptr = (*unknown_ptr->master)[0];	// potential for surface component



	// Include charge balance in list for mass-balance equations

	master_ptr->e->name.Copy(token);

	ptr = &token[0];

	if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[2]))

		return false;



  return true;

}



bool ModelEngine::MBForSpeciesAQ(int n)

{

	//

	// Make list of mass balance and charge balance equations in which

  //   to insert species n. 

  //

  // count_mb_unknowns - number of equations and summation relations

  // mb_unknowns.unknown - pointer to unknown which contains row number

  // mb_unknowns.source - pointer to the LDBLE number to be multiplied

  //                      by coef, usually moles.

  // mb_unknowns.coef - coefficient of s[n] in equation or relation



  int i, j;

	Species *s;

  Master *master_ptr;

  Unknown *unknown_ptr;



	md->mb_unknowns->Clear();



	s = gd->species_list[n];



	// e- does not appear in any mass balances

	if (s->type == EMINUS)

    return true;



	// Do not include diffuse layer in cb, alk, ah2o, mu

  if (md->charge_balance_unknown != NULL && s->type < H2O)

    StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);



	if (md->alkalinity_unknown != NULL && s->type < H2O)

    StoreMBUnknowns (md->alkalinity_unknown, &s->moles, s->alk, &s->dg);



	if (md->ah2o_unknown != NULL && s->type < H2O)

    StoreMBUnknowns (md->ah2o_unknown, &s->moles, 1.0, &s->dg);



	if (md->mu_unknown != NULL && s->type < H2O)

    StoreMBUnknowns (md->mu_unknown, &s->moles, s->z * s->z, &s->dg);



	// Include diffuse layer in hydrogen and oxygen mass balance

	if (md->mass_hydrogen_unknown != NULL)

  {

    if (md->dl_type_x != NO_DL && md->state >= REACTION)

      StoreMBUnknowns (md->mass_hydrogen_unknown, &s->tot_g_moles, s->h - 2 * s->o, &s->dg_total_g);

    else

      StoreMBUnknowns (md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);

  }



  if (md->mass_oxygen_unknown != NULL)

  {

    if (md->dl_type_x != NO_DL && md->state >= REACTION)

      StoreMBUnknowns (md->mass_oxygen_unknown, &s->tot_g_moles, s->o, &s->dg_total_g);

    else

      StoreMBUnknowns (md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);

  }



	// Sum diffuse layer charge into (surface + DL) charge balance

	SpeciesDiffLayer *sdl_p;

	if (md->use.sur_p != NULL && s->type < H2O && md->dl_type_x != NO_DL)

  {

    j = 0;

		for (i = 0; i < count_unknowns; i++)

    {

			unknown_ptr = (*unknown_list)[i];



      if (unknown_ptr->type == SURFACE_CB)

      {		

				if (md->use.sur_p->type == CD_MUSIC)

					unknown_ptr = (*unknown_list)[i + 2];

	

				sdl_p = (*s->diff_layer)[j];

				StoreMBUnknowns(unknown_ptr, &sdl_p->g_moles, s->z, &sdl_p->dg_g_moles);

				j++;

      }

    }

  }



	// Other mass balances

	ElementOfSpecies *eos;

	for (i = 0; i < eos_list->Count(); i++)

  {

		eos = (*eos_list)[i];



    if (eos->e->master->s->type > AQ && eos->e->master->s->type < SOLID)

      continue;



    master_ptr = eos->e->master;

    

		if (master_ptr->primary && master_ptr->s->secondary != NULL)

			master_ptr = master_ptr->s->secondary;



    if (master_ptr->u == md->ph_unknown)

      continue;

    else if (master_ptr->u == md->pe_unknown)

      continue;

    else if (master_ptr->u == md->charge_balance_unknown)

      continue;

    else if (master_ptr->u == md->alkalinity_unknown)

      continue;

    else if (master_ptr->u->type == SOLUTION_PHASE_BOUNDARY)

      continue;



		if (md->dl_type_x != NO_DL && md->state >= REACTION)

      StoreMBUnknowns(master_ptr->u, &s->tot_g_moles, eos->coef * master_ptr->coef, &s->dg_total_g);

    else

      StoreMBUnknowns(master_ptr->u, &s->moles, eos->coef * master_ptr->coef, &s->dg);

  }



  return true;

}



bool ModelEngine::StoreMBUnknowns(Unknown *unknown_ptr, LDBLE *LDBLE_ptr, LDBLE coef, LDBLE *gamma_ptr)

{

	//

	// Takes an unknown pointer and a coefficient and puts in

	// list of mb_unknowns

	//



  if (Equal(coef, (LDBLE)0.0, (LDBLE)TOLERANCE))

    return true;



	UnknownInfo *m = md->mb_unknowns->AddNew();



  m->u = unknown_ptr;

  m->source = LDBLE_ptr;

  m->gamma_source = gamma_ptr;

  m->coef = coef;



	return true;

}



bool ModelEngine::MBForSpeciesEX(int n)

{

	//

  // Make list of mass balance and charge balance equations in which

  // to insert exchange species n. 

  //

  //      count_mb_unknowns - number of equations and summation relations

  //      mb_unknowns.source - pointer to the LDBLE number to be multiplied

  //                           by coef, usually moles.

  //      mb_unknowns.unknown - pointer to unknown which contains row number

  //      mb_unknowns.coef - coefficient of s[n] in equation or relation

	//



  int i;

  Master *master_ptr;

	Species *s;



 	md->mb_unknowns->Clear();



	// Master species for exchange do not appear in any mass balances

	s = gd->species_list[n];

	if (s->type == EX && s->primary != NULL)

    return true;



	// Include diffuse layer in hydrogen and oxygen mass balance

  if (md->charge_balance_unknown != NULL)

		StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);



	if (md->mass_hydrogen_unknown != NULL)

		StoreMBUnknowns(md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);



	if (md->mass_oxygen_unknown != NULL)

    StoreMBUnknowns(md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);



	// Other mass balances

	ElementOfSpecies *eos_p;

  for (i = 0; i < eos_list->Count(); i++)

  {

		eos_p = (*eos_list)[i];

    if (eos_p->e->master->s->type > AQ && eos_p->e->master->s->type < SOLID)

      continue;



    master_ptr = eos_p->e->master;



		if (master_ptr->primary && master_ptr->s->secondary != NULL)

			master_ptr = master_ptr->s->secondary;

    

		// Special for ph_unknown, pe_unknown, and alkalinity_unknown

    if (master_ptr->u == md->ph_unknown)

      continue;

    else if (master_ptr->u == md->pe_unknown)

      continue;

    else if (master_ptr->u == md->alkalinity_unknown)

      continue;



		// EX, sum exchange species only into EXCH mass balance in initial calculation

		// into all mass balances in reaction calculation

    if (md->state >= REACTION || master_ptr->s->type == EX)

      StoreMBUnknowns(master_ptr->u, &s->moles, eos_p->coef * master_ptr->coef, &s->dg);

  }



  return true;

}



bool ModelEngine::MBForSpeciesSURF(int n)

{

	//

  // Make list of mass balance and charge balance equations in which

  // to insert species n. 

  //

  //      count_mb_unknowns - number of equations and summation relations

  //      mb_unknowns.source - pointer to the LDBLE number to be multiplied

  //                           by coef, usually moles.

  //      mb_unknowns.unknown - pointer to unknown which contains row number

  //      mb_unknowns.coef - coefficient of s[n] in equation or relation

	//



  int i;

  Master *master_ptr;



	md->mb_unknowns->Clear();



	// Include in charge balance, if diffuse_layer_x == false

	Species *s = gd->species_list[n];

  if (md->charge_balance_unknown != NULL && md->dl_type_x == NO_DL)

    StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);



	// Include diffuse layer in hydrogen and oxygen mass balance

  if (md->mass_hydrogen_unknown != NULL)

    StoreMBUnknowns (md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);



	if (md->mass_oxygen_unknown != NULL)

    StoreMBUnknowns (md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);



	// Other mass balances

	ElementOfSpecies *eos_p;

  for (i = 0; i < eos_list->Count(); i++)

  {

		eos_p = (*eos_list)[i];



		// Skip H+, e-, and H2O

    if (eos_p->e->master->s->type > AQ && eos_p->e->master->s->type < SOLID)

      continue;



    master_ptr = eos_p->e->master;

    if (master_ptr->primary && master_ptr->s->secondary != NULL)

			master_ptr = master_ptr->s->secondary;



		// SURF_PSI, sum surface species in (surface + DL) charge balance

    if (master_ptr->s->type == SURF_PSI && md->use.sur_p->type != CD_MUSIC)

    {

      StoreMBUnknowns (master_ptr->u, &s->moles, s->z, &s->dg);

      continue;

    }



    if (master_ptr->s->type == SURF_PSI && md->use.sur_p->type == CD_MUSIC)

    {

      StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[0], &s->dg);

      continue;

    }



    if (master_ptr->s->type == SURF_PSI1)

    {

      StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[1], &s->dg);

      continue;

    }



    if (master_ptr->s->type == SURF_PSI2)

    {

      StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[2], &s->dg);

      continue;

    }



		// Special for ph_unknown, pe_unknown, and alkalinity_unknown

    if (master_ptr->u == md->ph_unknown)

      continue;

    else if (master_ptr->u == md->pe_unknown)

      continue;

    else if (master_ptr->u == md->alkalinity_unknown)

      continue;



		// SURF, sum surface species only into SURFACE mass balance in initial calculation

		// into all mass balances in reaction calculation

    if (md->state >= REACTION || master_ptr->s->type == SURF)

    {

      StoreMBUnknowns (master_ptr->u, &s->moles, eos_p->coef * master_ptr->coef, &s->dg);

    }

  }



  return true;

}



bool ModelEngine::BuildMBSums()

{

	//

	// Function builds lists sum_mb1 and sum_mb2  that describe how to sum molalities

	// to calculate mass balance sums, including activity of water, ionic strength,

	// charge balance, and alkalinity.

	//



  int i;

  //LDBLE *target;

	UnknownInfo *ui;

	int count = md->mb_unknowns->Count();



	for (i = 0; i < count; i++)

  {

		ui = (*md->mb_unknowns)[i];

    StoreMB(ui->source, &(ui->u->f), ui->coef);



#ifdef DEBUG_MOHID

		if (d.debug_status)

			fprintf(d.buildjacobiansums_f, "BuildMBSums-1) i:%d %d %s %20.20e\n",

																			i,

																			ui->u->number,

																			ui->u->name.CharPtr(),

																			ui->coef);

#endif

  }



	return true;

}



bool ModelEngine::BuildJacobianSums(int k)

{

	//

	// Function builds lists sum_jacob1 and sum_jacob2 that describe how to sum molalities

	// to form jacobian.

	//



	int i, j, kk;

	int index;

  int count_g;

  LDBLE coef;

  //LDBLE *source, *target;

	UnknownInfo *mb_unknown;

	SpeciesDiffLayer *sdl_p;

	Unknown *u, *uprior, *u_kk;



	Species *s = gd->species_list[k];



	// Calculate jacobian coefficients for each mass balance equation

	for (i = 0; i < md->mb_unknowns->Count(); i++)

  {

		mb_unknown = (*md->mb_unknowns)[i];



		// Store d(moles) for a mass balance equation



		// initial solution only

    if (mb_unknown->u->type == SOLUTION_PHASE_BOUNDARY)

      continue;



		coef = mb_unknown->coef;



		StoreDN(k, mb_unknown->source, mb_unknown->u->number, coef, mb_unknown->gamma_source);



#ifdef DEBUG_MOHID

		if (d.debug_status)

			fprintf(d.buildjacobiansums_f, "1) k:%d i:%d u_number:%d: %20.20e\n",

																			k,

																			i, 

																			mb_unknown->u->number,

																			coef);

#endif



		// Add extra terms for change in dg/dx in diffuse layer model

		if (s->type >= H2O || md->dl_type_x == NO_DL)

      continue;

    else if ((mb_unknown->u->type == MB || mb_unknown->u->type == MH || mb_unknown->u->type == MH2O) && md->state >= REACTION)

    {

      if (md->mass_oxygen_unknown != NULL)

      {

				// term for water, sum of all surfaces

				index = mb_unknown->u->number * (count_unknowns + 1) + md->mass_oxygen_unknown->number;

				StoreJacob(&s->tot_dh2o_moles, &arr[index], coef);



#ifdef DEBUG_MOHID

				if (d.debug_status)

					fprintf(d.buildjacobiansums_f, "2) index:%d s_number:%d u_number:%d: %d %20.20e\n",

																					index,

																					s->number, 

																					md->mass_oxygen_unknown->number,

																					coef);

#endif

      }



			// terms for psi, one for each surface

      count_g = 0;

      for (j = 0; j < count_unknowns; j++)

      {

				u = (*unknown_list)[j];

				sdl_p = (*s->diff_layer)[count_g];



				if (u->type != SURFACE_CB)

					continue;



				index = mb_unknown->u->number * (count_unknowns + 1) + u->number;

				StoreJacob(&sdl_p->dx_moles, &arr[index], coef);

				count_g++;



#ifdef DEBUG_MOHID

				if (d.debug_status)

					fprintf(d.buildjacobiansums_f, "3) index:%d %d %d %d %20.20e\n",

																					index,

																					mb_unknown->u->number,

																					u->number,

																					count_g, 

																					coef);

#endif



				if (count_g >= md->use.sur_p->charge->Count())

					break;

      }



			// terms for related phases

      count_g = 0;

      for (j = 0; j < count_unknowns; j++)

      {

				u = (*unknown_list)[j];

				sdl_p = (*s->diff_layer)[count_g];



				if (u->type != SURFACE_CB)

					continue;



				uprior = (*unknown_list)[j - 1];



				// has related phase

				if (uprior->surface_comp->phase_name.IsEmpty())

					continue;



				// now find the related phase

				for (kk = count_unknowns - 1; kk >= 0; kk--)

				{

					u_kk = (*unknown_list)[kk];



					if (u_kk->type != PP)

						continue;



					if (u_kk->p->name == uprior->surface_comp->phase_name)

						break;

				}



				if (kk >= 0)

				{

					index = mb_unknown->u->number * (count_unknowns + 1) + (*unknown_list)[kk]->number;

					StoreJacob(&sdl_p->drelated_moles, &arr[index], coef);



#ifdef DEBUG_MOHID

					if (d.debug_status)

						fprintf(d.buildjacobiansums_f, "4) index:%d %d %d %20.20e\n",

																						index,

																						mb_unknown->u->number,

																						(*unknown_list)[kk]->number,

																						coef);

#endif

				}



				count_g++;

				if (count_g >= md->use.sur_p->charge->Count())

					break;

      }

    }

    else if (mb_unknown->u->type == SURFACE_CB)

    {

      count_g = 0;



			for (j = 0; j < count_unknowns; j++)

      {

				u = (*unknown_list)[j];

				sdl_p = (*s->diff_layer)[count_g];



				if (u->type != SURFACE_CB)

					continue;



				if (mb_unknown->u->number == u->number)

				{

					index = mb_unknown->u->number * (count_unknowns + 1) + u->number;

					StoreJacob(&sdl_p->dx_moles, &arr[index], coef);



					uprior = (*unknown_list)[j - 1];



					// term for related phase 

					// has related phase 

					if (!uprior->surface_comp->phase_name.IsEmpty())

					{

						// now find the related phase

						for (kk = count_unknowns - 1; kk >= 0; kk--)

						{

							u_kk = (*unknown_list)[kk];



							if (u_kk->type != PP)

								continue;



							if (u_kk->p->name == uprior->surface_comp->phase_name)

								break;

						}



						if (kk >= 0)

						{

							index = mb_unknown->u->number * (count_unknowns + 1) + u_kk->number;

							StoreJacob(&sdl_p->drelated_moles, &arr[index], coef);

						}

					}



					if (md->mass_oxygen_unknown != NULL)

					{

						// term for water, for same surfaces

						index = mb_unknown->u->number * (count_unknowns + 1) + md->mass_oxygen_unknown->number;

						StoreJacob(&sdl_p->dh2o_moles, &arr[index], coef);

					}

					break;

				}



				count_g++;



				if (count_g >= md->use.sur_p->charge->Count())

					break;

      }

    }

  }



  return true;

}



bool ModelEngine::WriteMBForSpeciesList(int n)

{

	//

	// Sets up data to add to species_list

	// Original secondary redox states are retained

	//



  int i;

	char *ptr, token[DEFAULT_STRING_LENGTH];

	ElementOfSpecies *new_eos;



	r_temp->Reset();

	r_temp->AddTRXN(gd->species_list[n]->rxn_s, 1.0, false);



	//Copy to elt_list

	eos_list->Clear();

  parent_count = 0;



	ReactionToken *r;

  for (i = 1; i < r_temp->token_list->Count(); i++)

  {

		r = r_temp->token_list->Element(i);



    if (r->s->secondary == NULL)

    {

			r->s->primary->e->name.Copy(token);

      ptr = &token[0];

			GetSecondaryElementsInSpecies(&ptr, r->coef);

    }

    else

    {

			r->s->secondary->e->name.Copy(token);

      ptr = &token[0];

			GetSecondaryElementsInSpecies(&ptr, r->coef);

    }

  }



	ElementOfSpecies *eos;

  for (i = 0; i < eos_list->Count(); i++)

  {

		eos = (*eos_list)[i];

    if (eos->e->name == "O(-2)")

    {

			new_eos = eos_list->AddNew();



			new_eos->e = gd->e_h_one;

      new_eos->coef = eos->coef * 2;

    }

  }



	CombineElements();



	eos_list->CopyTo(gd->species_list[n]->eos_sys_total);

	

  return true;

}



bool ModelEngine::BuildSpeciesList(int n)

{

	//

	// Builds a list that includes an entry for each master species in each

	// secondary reaction. Used for summing species of each element and 

	// printing results.

	//



  int j;

	SpeciesInfo *new_s;

	ElementOfSpecies *eos_p;

	Species *s = gd->species_list[n];



	// Treat species made only with H+, e-, and H2O specially 

	if (IsSpecial(s))

  {

		new_s = md->species_info_list->AddNew();



    new_s->master_s = gd->s_hplus;

    new_s->s = s;

    new_s->coef = 0.0;



		return true;

  }



	// Treat exchange species specially

  if (s->type == EX)

  {

    if (s->primary != NULL)

      return true;		// master species has md->zero molality



    for (j = 0; j < eos_list->Count(); j++)

    {

			eos_p = (*eos_list)[j];



      if (eos_p->e->master->s->type != EX)

				continue; 



			new_s = md->species_info_list->AddNew();



			new_s->master_s = eos_p->e->master->s;

			new_s->s = s;

			new_s->coef = eos_p->e->master->coef * eos_p->coef;

    }



    return true;

  }



	// Treat surface species specially

  if (s->type == SURF_PSI)

    return true;



  if (s->type == SURF)

  {

    for (j = 0; j < eos_list->Count(); j++)

    {

			eos_p = (*eos_list)[j];



      if (eos_p->e->master->s->type != SURF)

				continue;



      new_s = md->species_info_list->AddNew();



			new_s->master_s = eos_p->e->master->s;

			new_s->s = s;

			new_s->coef = eos_p->e->master->coef * eos_p->coef;

    }



    return true;

  }



	// Other aqueous species

	Master *master_ptr;

  for (j = 0; j < eos_list->Count(); j++)

  {

		eos_p = (*eos_list)[j];



    if (IsSpecial(eos_p->e->master->s))

      continue;



    if (eos_p->e->master->s->secondary != NULL)

      master_ptr = eos_p->e->master->s->secondary;

    else

      master_ptr = eos_p->e->master->s->primary;



    new_s = md->species_info_list->AddNew();



		new_s->master_s = master_ptr->s;

		new_s->s = s;

		new_s->coef = master_ptr->coef * eos_p->coef;

  }



  return true;

}



bool ModelEngine::StoreMB(LDBLE * source, LDBLE * target, LDBLE coef)

{

	//

	// Adds item to list sum_mb1 or sum_mb2

	// If coef is 1.0, adds to sum_mb1, which does not require a multiply

	// else, adds to sum_mb2, which will multiply by coef

	//



	STCoef *new_item;



  if (Equal(coef, (LDBLE)1.0, (LDBLE)TOLERANCE))

  {

		new_item = md->sum_mb1->AddNew();



		new_item->source = source;

		new_item->target = target;

  }

  else

  {

		new_item = md->sum_mb2->AddNew();



		new_item->source = source;

		new_item->target = target;

		new_item->coef = coef;

  }



  return true;

}



bool ModelEngine::WritePhaseSysTotal(int n)

{

	//

	// Sets up data to add to species_list

	// Original secondary redox states are retained

	//



  int i;

	ReactionToken *r;

	char *ptr, token[DEFAULT_STRING_LENGTH];

	Phase *p = gd->phase_list[n];



	// Start with secondary reaction

	r_temp->Reset();

	r_temp->AddTRXN(p->rxn_s, 1.0, false);



	// Copy to elt_list

	eos_list->Clear();

  parent_count = 0;



	for (i = 1; i < r_temp->token_list->Count(); i++)

  {

		r = r_temp->token_list->Element(i);



    if (r->s->secondary == NULL)

			r->s->primary->name.Copy(token); //master and element have the same name (change made in the load of database)

    else

			r->s->secondary->name.Copy(token);



    ptr = &token[0];

		GetSecondaryElementsInSpecies(&ptr, r->coef);

  }



	ElementOfSpecies *eos, *new_eos;

  for (i = 0; i < eos_list->Count(); i++)

  {

		eos = (*eos_list)[i];



    if (eos->e->name == "O(-2)")

    {

			new_eos = eos_list->AddNew();



      new_eos->e = gd->e_h_one;

      new_eos->coef = eos->coef * 2;

    }

  }



	CombineElements();



	eos_list->CopyTo(p->eos_sys_total);



	return true;

}



bool ModelEngine::BuildSolutionPhaseBoundaries()

{

	//

	// Build into sums the logic to calculate inverse saturation indices for

	// solution phase boundaries

	//



  int i, j;

  Master *master_ptr;

  ReactionToken *rxn_ptr;

	Unknown *x;



  if (md->solution_phase_boundary_unknown == NULL)

    return true;



	// Calculate inverse saturation index

  for (i = 0; i < count_unknowns; i++)

  {

		x = (*unknown_list)[i];



    if (x->type != SOLUTION_PHASE_BOUNDARY)

      continue;



    StoreMB(&x->p->lk, &x->f, 1.0);

    StoreMB(&x->si, &x->f, 1.0);



    if (!x->p->in)

			return false;



		for (j = 1; j < x->p->rxn_x->token_list->Count(); j++)

		{

			rxn_ptr = x->p->rxn_x->token_list->Element(j); 

      StoreMB(&rxn_ptr->s->la, &x->f, -rxn_ptr->coef);

		}

  }

	

	// Put coefficients into array

  for (i = 0; i < count_unknowns; i++)

  {

		x = (*unknown_list)[i];



    if (x->type != SOLUTION_PHASE_BOUNDARY)

      continue;



		for (j = 1; j < x->p->rxn_x->token_list->Count(); j++)

		{

			rxn_ptr = x->p->rxn_x->token_list->Element(j);



			if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)

				master_ptr = rxn_ptr->s->secondary;

      else

				master_ptr = rxn_ptr->s->primary;



			if (master_ptr->u == NULL)

				continue;



      StoreJacob0(x->number, master_ptr->u->number, rxn_ptr->coef);

    }

  }



  return true;

}







bool ModelEngine::BuildMinExch()

{

	//

	// Defines proportionality factor between mineral and exchanger to jacob0

	//



  int i, j, k, jj;

  int row;

  ExchComp *comp_ptr;

  Master *master_ptr, *m0;

  Unknown *unknown_ptr, *u_j, *u_k;

	ElementOfSpecies *eos_p;

  LDBLE coef;



	if (md->use.exc_p == NULL)

    return true;

		

  if (!md->use.exc_p->related_phases)

    return true;

		

	for (i = 0; i < md->use.exc_p->comps->Count(); i++)

  {

		comp_ptr = (*md->use.exc_p->comps)[i];



		if (comp_ptr->phase_name.IsEmpty())

      continue;

			

    // find unknown number

    for (j = count_unknowns - 1; j >= 0; j--)

    {

			u_j = (*unknown_list)[j];

			m0 = (*u_j->master)[0];



			if (u_j->type != EXCH)

				continue;

				

      if (m0 == comp_ptr->master)

				break;

    }

		

    for (k = count_unknowns - 1; k >= 0; k--)

    {

			u_k = (*unknown_list)[k];



			if (u_k->type != PP)

				continue;

				

      if (u_k->p->name == comp_ptr->phase_name)

				break;

    }

		

    if (j == -1)

			return false;

		

    if (k == -1)

      continue;

			

		// Build jacobian



    // charge balance

    StoreJacob0(md->charge_balance_unknown->number, u_k->number, comp_ptr->formula_z * comp_ptr->phase_proportion);

		StoreSumDeltas(&delta[k], &md->charge_balance_unknown->delta, -comp_ptr->formula_z * comp_ptr->phase_proportion);

		

    // mole balance balance

    eos_list->Clear();

    parent_count = 0;

    

		CopyToTempEOSList(comp_ptr->formula_totals, 1.0);

    ChangeHydrogenInEOSList(0);

		

    for (jj = 0; jj < eos_list->Count(); jj++)

    {

			eos_p = (*eos_list)[jj];

      master_ptr = (*eos_list)[jj]->e->primary;

			

			if (!master_ptr->in)

				master_ptr = master_ptr->s->secondary;

			

      if (master_ptr == NULL)

				return false;

			

      if (master_ptr->s->type == EX)

      {

				if (!Equal(u_j->moles, u_k->moles * eos_p->coef * comp_ptr->phase_proportion, (LDBLE)5.0 * md->convergence_tolerance))

				  u_j->moles = u_k->moles * eos_p->coef * comp_ptr->phase_proportion;

      }

			

      coef = eos_p->coef;

			

      if (master_ptr->s == gd->s_hplus)

      {

				row = md->mass_hydrogen_unknown->number;

				unknown_ptr = md->mass_hydrogen_unknown;

      }

      else if (master_ptr->s == gd->s_h2o)

      {

				row = md->mass_oxygen_unknown->number;

				unknown_ptr = md->mass_oxygen_unknown;

      }

      else

      {

				row = master_ptr->u->number;

				unknown_ptr = master_ptr->u;

      }

			

      StoreJacob0 (row, u_k->number, coef * comp_ptr->phase_proportion);

      StoreSumDeltas(&delta[k], &unknown_ptr->delta, -coef * comp_ptr->phase_proportion);

    }

  }

	

  return true;

}



bool ModelEngine::BuildMinSurface()

{

	//

	//   Defines proportionality factor between mineral and surface to jacob0

	//



  int i, j, k, jj, row;

  List<ElementOfSpecies> *next_elt;

	ElementOfSpecies *eos_p;

  SurfaceComp *comp_ptr, *sc_p;

  Unknown *unknown_ptr, *u_j, *u_k, *u;

  Master *master_ptr, *m;

  LDBLE coef;



  if (md->use.sur_p == NULL)

    return true;

		

	if (!md->use.sur_p->related_phases)

		return true;

		

	for (i = 0; i < md->use.sur_p->comps->Count(); i++)

  {

		sc_p = (*md->use.sur_p->comps)[i];



		if (sc_p->phase_name.IsEmpty())

      continue;

			

    // find unknown number

    for (j = count_unknowns - 1; j >= 0; j--)

    {

			u_j = (*unknown_list)[j];

			m = (*u_j->master)[0];



      if (u_j->type != SURFACE)

				continue;

				

      if (m == sc_p->master)

				break;

    }

		

    for (k = count_unknowns - 1; k >= 0; k--)

    {

			u_k = (*unknown_list)[k];



      if (u_k->type != PP)

				continue;

				

      if (u_k->p->name == sc_p->phase_name)

				break;

    }

		

    if (j == -1)

			return false;

		

    if (k == -1)

      continue;

			

    comp_ptr = u_j->surface_comp;

		

    if (md->use.sur_p->type == CD_MUSIC)

    {

      // Add formula for CD_MUSIC

      next_elt = comp_ptr->formula_totals;

    }

    else

    {

      // Add master species for non CD_MUSIC

			next_elt = (*u_j->master)[0]->s->eos_list;

    }

		

    // update grams == moles in this case

		u = (*unknown_list)[j + 1];

    if (j < count_unknowns - 1 && u->type == SURFACE_CB)

      StoreSumDeltas (&delta[k], &u->related_moles, -1.0);

		

    // charge balance

    StoreJacob0 (md->charge_balance_unknown->number, u_k->number, comp_ptr->formula_z * comp_ptr->phase_proportion);

    StoreSumDeltas (&delta[k], &md->charge_balance_unknown->delta, -comp_ptr->formula_z * comp_ptr->phase_proportion);

		

    eos_list->Clear();

    parent_count = 0;

		

    CopyToTempEOSList(next_elt, 1.0);

    ChangeHydrogenInEOSList(0);

		

		for (jj = 0; jj < eos_list->Count(); jj++)

    {

			eos_p = (*eos_list)[jj];



      master_ptr = eos_p->e->primary;

      if (!master_ptr->in)

				master_ptr = master_ptr->s->secondary;



			if (master_ptr == NULL)

				return false;



      if (master_ptr->s->type == SURF)

      {

				if (!Equal(u_j->moles, u_k->moles * eos_p->coef * comp_ptr->phase_proportion, (LDBLE)5.0 * md->convergence_tolerance))

				  u_j->moles = u_k->moles * eos_p->coef * comp_ptr->phase_proportion;

      }

			

      coef = eos_p->coef;

			

      if (master_ptr->s == gd->s_hplus)

      {

				row = md->mass_hydrogen_unknown->number;

				unknown_ptr = md->mass_hydrogen_unknown;

      }

      else if (master_ptr->s == gd->s_h2o)

      {

				row = md->mass_oxygen_unknown->number;

				unknown_ptr = md->mass_oxygen_unknown;

      }

      else

      {

				row = master_ptr->u->number;

				unknown_ptr = master_ptr->u;

      }

			

      StoreJacob0(row, u_k->number, coef * comp_ptr->phase_proportion);

      StoreSumDeltas(&delta[k], &unknown_ptr->delta, -coef * comp_ptr->phase_proportion);

    }

  }

	

  return true;

}



bool ModelEngine::BuildGasPhase()

{

	//

	// Put coefficients into lists to sum iaps to test for equilibrium

	// Put coefficients into lists to build jacobian for 

	// sum of partial pressures equation and

	// mass balance equations for elements contained in gases

	//



  int i, j, k;

  int row, col;

  Master *master_ptr;

  ReactionToken *rxn_ptr;

  GasComp *gas_comp_ptr;

  Phase *phase_ptr;

  Unknown *unknown_ptr;

  LDBLE coef, coef_elt;

	ElementOfSpecies *eos_p;



  if (md->gas_unknown == NULL)

    return true;

		

	for (i = 0; i < md->use.gas_p->comps->Count(); i++)

  {

		// Determine elements in gas component

    eos_list->Clear();

    parent_count = 0;



    gas_comp_ptr = (*md->use.gas_p->comps)[i];

    phase_ptr = gas_comp_ptr->phase;

		

		if (phase_ptr->rxn_x->token_list->Count() <= 0)

      continue;

			

    CopyToTempEOSList(phase_ptr->eos_list, 1.0);

    ChangeHydrogenInEOSList(0);

		

		// Build mass balance sums for each element in gas



		// All elements in gas

    for (j = 0; j < eos_list->Count(); j++)

    {

			eos_p = (*eos_list)[j];



      unknown_ptr = NULL;

			if (eos_p->e->name == "H")

				unknown_ptr = md->mass_hydrogen_unknown;

			else if (eos_p->e->name == "O")

				unknown_ptr = md->mass_oxygen_unknown;

      else

      {

				if (eos_p->e->primary->in)

				  unknown_ptr = eos_p->e->primary->u;

				else if (eos_p->e->primary->s->secondary != NULL)

				  unknown_ptr = eos_p->e->primary->s->secondary->u;

      }

			

      if (unknown_ptr != NULL)

      {

				coef = eos_p->coef;

				StoreMB(&(gas_comp_ptr->phase->moles_x), &(unknown_ptr->f), coef);

      }

    }

		

    if (md->use.gas_p->type == PRESSURE)

    {

      // Total pressure of gases

      StoreMB(&(gas_comp_ptr->phase->p_soln_x), &(md->gas_unknown->f), 1.0);

    }

		

		// Build jacobian sums for mass balance equations

    for (j = 0; j < eos_list->Count(); j++)

    {

			eos_p = (*eos_list)[j];



      unknown_ptr = NULL;

			

			if (eos_p->e->name == "H")

				unknown_ptr = md->mass_hydrogen_unknown;

			else if (eos_p->e->name == "O")

				unknown_ptr = md->mass_oxygen_unknown;

      else

      {

				if (eos_p->e->primary->in)

				  unknown_ptr = eos_p->e->primary->u;

				else if (eos_p->e->primary->s->secondary != NULL)

				  unknown_ptr = eos_p->e->primary->s->secondary->u;

      }

			

      if (unknown_ptr == NULL)

				continue;

			

      row = unknown_ptr->number * (count_unknowns + 1);

      coef_elt = eos_p->coef;

			

			for (k = 1; k < phase_ptr->rxn_x->token_list->Count(); k++)

      {

				rxn_ptr = phase_ptr->rxn_x->token_list->Element(k);



				if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)

				  master_ptr = rxn_ptr->s->secondary;

				else

				  master_ptr = rxn_ptr->s->primary;

				

				if (master_ptr == NULL)

					return false;

				

				if (master_ptr->u == NULL)

				  continue;

				

				if (!master_ptr->in)

					return false;

				

				col = master_ptr->u->number;

				coef = coef_elt * rxn_ptr->coef;

				StoreJacob (&(gas_comp_ptr->phase->moles_x), &(arr[row + col]), coef);				

      }

			

      if (md->use.gas_p->type == PRESSURE)

      {

				// derivative wrt total moles of gas

				StoreJacob(&(gas_comp_ptr->phase->fraction_x), &(arr[row + md->gas_unknown->number]), coef_elt);

      }

    }

		

		// Build jacobian sums for sum of partial pressures equation

    if (md->use.gas_p->type != PRESSURE)

      continue;

			

    unknown_ptr = md->gas_unknown;

    row = unknown_ptr->number * (count_unknowns + 1);

		

		for (k = 1; k < phase_ptr->rxn_x->token_list->Count(); k++)

    {

			rxn_ptr = phase_ptr->rxn_x->token_list->Element(k);



			if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)

				master_ptr = rxn_ptr->s->secondary;

      else

				master_ptr = rxn_ptr->s->primary;

			

      if (master_ptr->u == NULL)

				continue;



			if (!master_ptr->in)

				return false;

			

      col = master_ptr->u->number;

      coef = rxn_ptr->coef;

      StoreJacob (&(gas_comp_ptr->phase->p_soln_x), &(arr[row + col]), coef);

    }

  }

	

  return true;

}



bool ModelEngine::BuildSSAssemblage()

{

	//

	// Put coefficients into lists to sum iaps to test for equilibrium

	// Put coefficients into lists to build jacobian for 

	//    mass action equation for component

	//    mass balance equations for elements contained in solid solutions

	//



  int i, j, k, l, m;

	bool stop;

  int row, col;

  Master *master_ptr, *m_2, *m0;

  ReactionToken *rxn_ptr;

	Unknown *u, *u2;

  SS *s_s_ptr, *s_s_ptr_old;

	//SSComp *ssc_p;

	ElementOfSpecies *eos_p;

  char token[MAX_LENGTH];

  char *ptr;



  if (md->s_s_unknown == NULL)

    return true;

		

  s_s_ptr_old = NULL;

  col = 0;

	

  for (i = 0; i < count_unknowns; i++)

  {

		u = (*unknown_list)[i];



    if (u->type != S_S_MOLES)

      continue;

			

    s_s_ptr = u->s_s;



    if (s_s_ptr != s_s_ptr_old)

    {

      col = u->number;

      s_s_ptr_old = s_s_ptr;

 …