Added code to adjust activities of master species in solution when totals are adjusted by SOLUTION_MODIFY.

git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@5270 1feff8c3-07ed-0310-ac33-dd36852eb9cd

Phreeqc.h

@@ -1082,6 +1082,7 @@ int master_delete(char *ptr);
 public:
 struct master *master_bsearch(const char *ptr);
 struct master *master_bsearch_primary(char *ptr);
+struct master *master_bsearch_secondary(char *ptr);
 private:
 struct master *master_search(char *ptr, int *n);
 struct mix *mix_bsearch(int k, int *n);

ReadClass.cxx

@@ -1155,6 +1155,7 @@ read_run_cells(void)
 	if (return_value == OPTION_KEYWORD) output_msg(OUTPUT_CHECKLINE, "\t%s\n", line);
 	return (return_value);
 }
+#ifdef SKIP
 /* ---------------------------------------------------------------------- */
 int CLASS_QUALIFIER
 read_solution_modify(void)
@@ -1244,6 +1245,105 @@ read_solution_modify(void)
 	if (return_value == OPTION_KEYWORD) output_msg(OUTPUT_CHECKLINE, "\t%s\n", line);
 	return (return_value);
 }
+#endif
+/* ---------------------------------------------------------------------- */
+int CLASS_QUALIFIER
+read_solution_modify(void)
+/* ---------------------------------------------------------------------- */
+{
+/*
+ *      Reads solution_modify data block
+ *
+ *      Arguments:
+ *         none
+ *
+ *      Returns:
+ *         KEYWORD if keyword encountered, input_error may be incremented if
+ *                    a keyword is encountered in an unexpected position
+ *         EOF     if eof encountered while reading mass balance concentrations
+ *         ERROR   if error occurred reading data
+ *
+ */
+	int return_value;
+
+	// find solution number
+	char token[MAX_LENGTH];
+	char *next;
+	int l, n_user, n;
+	next = line;
+	copy_token(token, &next, &l);
+	if (copy_token(token, &next, &l) != DIGIT)
+	{
+		input_error++;
+		sprintf(error_string, "Expected solution number following SOLUTION_MODIFY.\n%s\n", line_save);
+		error_msg(error_string, CONTINUE);
+		std::istringstream iss_in;
+		return_value = streamify_to_next_keyword(iss_in);
+		return (ERROR);
+	} 
+	else
+	{
+		sscanf(token,"%d", &n_user);
+	}	
+	/*
+	 *  Make parser
+	 */
+	std::istringstream iss_in;
+	return_value = streamify_to_next_keyword(iss_in);
+	std::ostringstream oss_out;
+	std::ostringstream oss_err;
+	CParser parser(PHREEQC_THIS_COMMA iss_in, oss_out, oss_err);
+	assert(!reading_database());
+
+	//For testing, need to read line to get started
+	parser.set_echo_file(CParser::EO_NONE);
+	std::vector < std::string > vopts;
+	std::istream::pos_type next_char;
+	parser.get_option(vopts, next_char);
+
+	if (pr.echo_input == FALSE)
+	{
+		parser.set_echo_file(CParser::EO_NONE);
+	}
+	else
+	{
+		parser.set_echo_file(CParser::EO_NOKEYWORDS);
+	}
+
+
+	if (solution_bsearch(n_user, &n, FALSE) == NULL)
+	{
+		input_error++;
+		sprintf(error_string, "Solution %d not found for SOLUTION_MODIFY.\n", n_user);
+		error_msg(error_string, CONTINUE);
+		return (ERROR);
+	}
+
+	cxxSolution sol(solution[n]);
+
+	// Clear activities so we can know what was read
+	sol.clear_master_activity();
+
+	sol.read_raw(PHREEQC_THIS_COMMA parser, false);
+
+	cxxSolution orig(solution[n]);
+
+	sol.modify_activities(PHREEQC_THIS_COMMA orig);
+
+	struct solution *soln_ptr = sol.cxxSolution2solution(PHREEQC_THIS);
+
+	/*
+	 *  This is not quite right, may not produce sort order, forced sort
+	 */
+
+	solution_free(solution[n]);
+	solution[n] = soln_ptr;
+
+	// Need to output the next keyword
+	if (return_value == OPTION_KEYWORD) output_msg(OUTPUT_CHECKLINE, "\t%s\n", line);
+	return (return_value);
+}
+
 /* ---------------------------------------------------------------------- */
 int CLASS_QUALIFIER
 read_equilibrium_phases_modify(void)

Solution.cxx

@@ -1545,7 +1545,143 @@ cxxSolution::set_master_activity(char *string, double d)
 {
 	this->master_activity[string] = d;
 }
+void
+cxxSolution::modify_activities(PHREEQC_PTR_ARG_COMMA const cxxSolution & original)
+//
+// Estimate activities after solution_modify
+//
+{
+	// Note: any master_activities in "this" have been read in SOLUTION_MODIFY
 
+	// to standardize, convert element to valence state if needed
+	// for original activity list (probably not needed)
+	cxxNameDouble orig_master_activity(original.get_master_activity());
+	cxxNameDouble::const_iterator it;
+	bool redo=true;
+	while (redo)
+	{
+		redo = false;
+		for (it = orig_master_activity.begin(); it != orig_master_activity.end(); it++)
+		{
+			struct master *master_ptr = P_INSTANCE_POINTER master_bsearch(it->first.c_str());
+			if (master_ptr != NULL)
+			{
+				if (master_ptr->primary == TRUE)
+				{
+					struct master *master_ptr_secondary = P_INSTANCE_POINTER master_bsearch_secondary(master_ptr->elt->name);
+
+					// redox element erase and replace
+					if (master_ptr_secondary != master_ptr) 
+					{
+						double d = it->second;
+						orig_master_activity.erase(orig_master_activity.find(master_ptr->elt->name));
+						orig_master_activity[master_ptr_secondary->elt->name] = d;
+						redo = true;
+						break;
+					}
+				}
+			}
+			else
+			{
+				P_INSTANCE_POINTER error_msg("Could not find master species in modify_activities.", STOP);
+			}
+		}
+	}
+
+	// also for modified activity list
+	cxxNameDouble mod_master_activity(this->master_activity);
+	redo=true;
+	while (redo)
+	{
+		redo = false;
+		for (it = mod_master_activity.begin(); it != mod_master_activity.end(); it++)
+		{
+			struct master *master_ptr = P_INSTANCE_POINTER master_bsearch(it->first.c_str());
+			if (master_ptr != NULL)
+			{
+				if (master_ptr->primary == TRUE)
+				{
+					struct master *master_ptr_secondary = P_INSTANCE_POINTER master_bsearch_secondary(master_ptr->elt->name);
+
+					// redox element erase and replace
+					if (master_ptr_secondary != master_ptr) 
+					{
+						double d = it->second;
+						mod_master_activity.erase(mod_master_activity.find(master_ptr->elt->name));
+						mod_master_activity[master_ptr_secondary->elt->name] = d;
+						redo = true;
+						break;
+					}
+				}
+			}
+			else
+			{
+				P_INSTANCE_POINTER error_msg("Could not find master species in modify_activities.", STOP);
+			}
+		}
+	}
+
+	// go through totals
+	double d = 0.0;
+	for (it = this->totals.begin(); it != this->totals.end(); ++it)
+	{
+		// find element name 
+		struct master *master_ptr = P_INSTANCE_POINTER master_bsearch(it->first.c_str());
+		char * ename = master_ptr->elt->primary->elt->name;
+		char * secondary_name;
+		if (master_ptr->primary == TRUE)
+		{
+			struct master *m_ptr = P_INSTANCE_POINTER master_bsearch_secondary(ename);
+			secondary_name = m_ptr->elt->name;
+		}
+		else
+		{
+			secondary_name = master_ptr->elt->name;
+		}
+
+		double d_mod, d_orig;
+		d_mod = this->get_total_element(ename);
+		if (d_mod <= 0) continue;
+
+		d_orig = original.get_total_element(ename);
+		if (d_orig <= 0) 
+		{
+			// add estimate for la based on concentration if not in list
+			if (mod_master_activity.find(secondary_name) == mod_master_activity.end())
+			{
+				mod_master_activity[secondary_name] = log10(d_mod) - 2.0;
+			}
+			continue;
+		}
+
+		// case where total for both orig and modified are greater than 0
+		double lratio = log10(d_mod / d_orig);
+
+		if (mod_master_activity.find(secondary_name) == mod_master_activity.end())
+		{
+			cxxNameDouble::const_iterator it1;
+			it1 = orig_master_activity.find(secondary_name);
+			if (it1 != orig_master_activity.end())
+			{
+				double d = it1->second;
+				mod_master_activity[secondary_name] = d + lratio;
+			}
+			else
+			// Has total, but no activity, should not happen
+			{
+				mod_master_activity[secondary_name] = log10(d_mod) - 2.0;
+			}
+		}
+	}
+
+	// merge activities
+	this->master_activity = orig_master_activity;
+
+	for (it = mod_master_activity.begin(); it != mod_master_activity.end(); it++)
+	{
+		this->master_activity[it->first] = it->second;
+	}
+}
 //#include "ISolution.h"
 //#include "Exchange.h"
 //#include "Surface.h"

Solution.h

@@ -132,6 +132,10 @@ class cxxSolution:public cxxNumKeyword
 	{
 		this->totals.clear();
 	}
+	void clear_master_activity()
+	{
+		this->master_activity.clear();
+	}
 
 	const cxxNameDouble & get_master_activity(void) const
 	{
@@ -154,6 +158,7 @@ class cxxSolution:public cxxNumKeyword
 
 	void read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check = true);
 	void multiply(double extensive);
+	void modify_activities(PHREEQC_PTR_ARG_COMMA const cxxSolution & original);
 #ifdef ORCHESTRA
 	void ORCH_write(std::ostream & headings, std::ostream & input_data) const;
 	void ORCH_read(std::vector < std::pair < std::string,