Roughed in surface area calculation by DeClercq (sa_declercq for now).

Fixed bug where H+ was missing from list of aqueous species generated by SYS.

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

PBasic.cpp

@@ -1593,6 +1593,9 @@ listtokens(FILE * f, tokenrec * l_buf)
 		case tokdiff_c:
 			output_msg("DIFF_C");
 			break;
+		case toksa_declercq:
+			output_msg("SA_DECLERCQ");
+			break;
 		}
 		l_buf = l_buf->next;
 	}
@@ -3531,6 +3534,42 @@ factor(struct LOC_exec * LINK)
 
 		}
 		break;
+
+	case toksa_declercq:
+		{
+			double type, sa, d, m, m0, gfw;
+
+			// left parenthesis
+			require(toklp, LINK);
+
+			// first double arugument, type
+			type = realfactor(LINK);
+			require(tokcomma, LINK);
+
+			// second double arugument, Sa
+			sa = realfactor(LINK);
+			require(tokcomma, LINK);
+			
+			// third double arugument, Sa
+			d = realfactor(LINK);
+			require(tokcomma, LINK);
+
+			// fourth double arugument, m
+			m = realfactor(LINK);
+			require(tokcomma, LINK);
+
+			// fifth double arugument, m0
+			m0 = realfactor(LINK);
+			require(tokcomma, LINK);
+
+			// sixth double arugument, gfw
+			gfw = realfactor(LINK);
+			require(tokrp, LINK);
+
+			n.UU.val = (parse_all) ? 1 : PhreeqcPtr->sa_declercq(type, sa, d, m, m0, gfw);
+		}
+		break;
+
 	case tokdiff_c:
 		{
 			const char * str = stringfactor(STR1, LINK);
@@ -6830,7 +6869,8 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("eq_frac",            PBasic::tokeq_frac),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("equiv_frac",         PBasic::tokeq_frac),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("callback",           PBasic::tokcallback),
-	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("diff_c",             PBasic::tokdiff_c)
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("diff_c",             PBasic::tokdiff_c),
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sa_declercq",        PBasic::toksa_declercq)
 };
 std::map<const std::string, PBasic::BASIC_TOKEN> PBasic::command_tokens(temp_tokens, temp_tokens + sizeof temp_tokens / sizeof temp_tokens[0]);
 

PBasic.h

@@ -322,7 +322,8 @@ public:
 		tokeq_frac,
 		tokequiv_frac,
 		tokcallback,
-		tokdiff_c
+		tokdiff_c,
+		toksa_declercq
 	};
 
 #if !defined(PHREEQCI_GUI)

Phreeqc.h

@@ -93,6 +93,7 @@ public:
 	LDBLE log_activity_coefficient(const char *species_name);
 	LDBLE aqueous_vm(const char *species_name);
 	LDBLE diff_c(const char *species_name);
+	LDBLE sa_declercq(double type, double sa, double d, double m, double m0, double gfw);
 	LDBLE calc_SC(void);
 	/* VP: Density Start */
 	LDBLE calc_dens(void);

basicsubs.cpp

@@ -98,6 +98,58 @@ aqueous_vm(const char *species_name)
 	}
 	return (g);
 }
+
+/* ---------------------------------------------------------------------- */
+LDBLE Phreeqc::
+sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw)
+/* ---------------------------------------------------------------------- */
+{
+	if (sa_type == 0) 
+	{
+		// surface-area-calculation-Fixed_Surface
+		return Sa;
+	}
+	else if (sa_type == 1)
+		// surface-area-calculation-Square
+	{
+		double mass0 = m0 * gfw;
+		double V0 = mass0 / d;
+		double St0 = mass0 * Sa;        // total surface
+		double a0 = pow(V0, 1.0/3.0);   // side length
+		double Sp0 = 6.0 * a0*a0;       // surface particle
+		double np = St0 / Sp0;          // number of particles
+	    double RATS = Sa / St0;
+		double mass = m * gfw;
+		double V = mass / d;
+		double a = pow(V, 1.0/3.0); 
+		double St = 6.0 * a*a*np;
+		return St * RATS;               // total current surface
+	}
+	else if (sa_type == 2)
+	{
+		//double pi = 3.14159265359;
+		double mass0 = m0 * gfw;
+		double V0 = mass0 / d;                         // volume
+		double St0 = mass0 * Sa;                       // total surface
+		double a0 = pow(3.0 * V0/(4.0 * pi), 1.0/3.0); // ((3*V0)/(4 * 3.14159265359))^(1/3)  
+		double Sp0 = (4.0 * pi) * a0 * a0;             // surface particle
+		double np = St0 / Sp0;                         // number of particles
+		double RATS = Sa / St0;
+ 
+		double mass = m * gfw;
+		double V = mass / d;
+		double a = pow(3.0 * V/(4.0 * pi), 1.0/3.0);  //((3*V)/(4 * 3.14159265359))^(1/3)
+		double St = 4.0 * pi * a * a * np;
+		return St * RATS;                             // total current surface
+	}
+	error_string = sformatf( "Unknown surface area type in SA_DECLERCQ %d.", (int) sa_type);
+	error_msg(error_string, CONTINUE);
+	input_error++;
+	return (MISSING);
+ 
+}
+
+/* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
 diff_c(const char *species_name)
 /* ---------------------------------------------------------------------- */
@@ -2698,7 +2750,7 @@ system_total_aq(void)
 	for (i = 0; i < count_s_x; i++)
 	{
 		//if (s_x[i]->type != AQ)
-		if (s_x[i]->type > AQ)
+		if (s_x[i]->type > HPLUS)
 			continue;
 		sys[count_sys].name = string_duplicate(s_x[i]->name);
 		sys[count_sys].moles = s_x[i]->moles;