New solver uses VERY SMALL time steps. The circuits have some very fast...

New solver uses VERY SMALL time steps. The circuits have some very fast transients that appear to be a hard to resolve source of numerical instability. Model automatically adds resistance and capacitance to electrical lines when these are absent. Added a kludgly exictation controller.

src/ElectricalData.cpp

@@ -32,28 +32,34 @@ void ElectricalData::configureGenerator(genr_t& genr_data)
 	 */
 	if (genr_data.genr_mw * SBase > 600.0) // Big nuke plant
 	{
-		genr_data.H = 3.5*genr_data.genr_mw;
+		genr_data.H = 3.5;
+		genr_data.Lg = 0.26/(2.0*pi*freq);
 	}
 	else if (genr_data.genr_mw * SBase > 300.0) // Big fossil plant
 	{
-		genr_data.H = 4.0*genr_data.genr_mw;
+		genr_data.H = 4.0;
+		genr_data.Lg = 0.25/(2.0*pi*freq);
 	}
 
 	else if (genr_data.genr_mw * SBase > 150.0) // Smaller steam plant
 	{
-		genr_data.H = 6.0*genr_data.genr_mw;
+		genr_data.H = 6.0;
+		genr_data.Lg = 0.26/(2.0*pi*freq);
 	}
 	else 
 	{
 		genr_data.H = 2.0;
+		genr_data.Lg = 0.2/(2.0*pi*freq);
 	}
 	// Try to get close to an appropriate set of initial conditions
 	genr_data.w0 = 0.0;
 	genr_data.Tspd1 = 10.0;
 	genr_data.Tspd2 = 1.0;
 	genr_data.R = 1.0/0.05;
-	genr_data.Agc = genr_data.R/100.0;
-	genr_data.D = 1.0;
+	genr_data.Agc = genr_data.R/10.0;
+	genr_data.K = 100.0;
+	genr_data.G = 400.0;
+	genr_data.Rg = 0.0;
 }
 
 void ElectricalData::initialize_circuit()
@@ -64,12 +70,27 @@ void ElectricalData::initialize_circuit()
 	for (unsigned i = 0; i < lines.size(); i++)
 	{
 		assert(real(lines[i].B) == 0.0);
+		// Some capacitance belongs to every line
+		if (lines[i].B == Complex(0.0,0.0))
+		{
+			double L = imag(1.0/lines[i].y)/(2.0*pi*freq);
+			double C = L*1E-2;
+			lines[i].B = Complex(0.0,C*2.0*pi*freq);
+		}
+		// Some resistance belongs to every line
+		if (real(lines[i].y) == 0.0)
+		{
+			lines[i].Z = 1.0/lines[i].y;
+			lines[i].Z = Complex(imag(lines[i].Z)/(6.0*pi*freq),imag(lines[i].Z));
+			lines[i].y = 1.0/lines[i].Z;
+		}
+		double C = imag(lines[i].B)/(2.0*pi*freq);
 		for (unsigned j = 0; j < nodes.size(); j++)
 		{
 			if (lines[i].from == nodes[j].ID)
-				nodes[j].B += lines[i].B/2.0;
+				nodes[j].Cshunt += C/2.0;
 			else if (lines[i].to == nodes[j].ID)
-				nodes[j].B += lines[i].B/2.0;
+				nodes[j].Cshunt += C/2.0;
 		}
 		// Calculate line parameters
 		lines[i].B = Complex(0.0,0.0);
@@ -83,7 +104,6 @@ void ElectricalData::initialize_circuit()
 	{
 		node_to_idx[nodes[i].ID] = i;
 		nodes[i].y = nodes[i].Z = Complex(0.0,0.0);
-		nodes[i].R = nodes[i].L = nodes[i].C = nodes[i].Cshunt = 0.0;
 		if (nodes[i].B != Complex(0.0,0.0))
 		{
 			assert(real(nodes[i].B) == 0.0);
@@ -91,7 +111,7 @@ void ElectricalData::initialize_circuit()
 				nodes[i].y += nodes[i].B;
 			// Charging as capacitance
 			else
-				nodes[i].Cshunt = imag(nodes[i].B)/(2.0*pi*freq);
+				nodes[i].Cshunt += imag(nodes[i].B)/(2.0*pi*freq);
 		}
 		Complex V = polar(nodes[i].vamp,nodes[i].theta);
 		Complex S(nodes[i].load_mw,nodes[i].load_mvar);
@@ -107,10 +127,12 @@ void ElectricalData::initialize_circuit()
 			nodes[i].L = imag(nodes[i].Z)/(2.0*pi*freq);
 		else if (imag(nodes[i].Z) < 0.0)
 		{
-			nodes[i].C = -1.0/(imag(nodes[i].Z)*(2.0*pi*freq));
+			nodes[i].C = -1.0/(2.0*pi*freq*imag(nodes[i].Z));
 			assert(nodes[i].C > 0.0);
 		}
-		// Lump a purely capacitive load into the shunt capacitor
+		// Lump a purely capacitive load into the shunt capacitor or, 
+		// if there is no load, add a small load that will dissipate
+		// high frequency signals
 		if (nodes[i].R == 0.0 && nodes[i].L == 0.0)
 		{
 			nodes[i].Cshunt += nodes[i].C;
@@ -118,11 +140,8 @@ void ElectricalData::initialize_circuit()
 			nodes[i].y = Complex(0.0,0.0);
 			nodes[i].Z = 1.0/nodes[i].y;
 		}
-		// Must have a small shunt capacitance
 		if (nodes[i].Cshunt == 0.0)
-		{
 			nodes[i].Cshunt = 1E-5;
-		}
 		nodes[i].B = Complex(0,2.0*pi*freq*nodes[i].Cshunt);
 	}
 	// Find initial generator currents and voltages
@@ -216,6 +235,7 @@ void ElectricalData::initialize_circuit()
 	{
 		if (nodes[entry.second].genr.genr_mw == 0.0 && nodes[entry.second].genr.genr_mvar == 0.0)
 			continue;
+		configureGenerator(nodes[entry.second].genr);
 		Complex Igenr = nodes[entry.second].V*(nodes[entry.second].B+nodes[entry.second].y);
 		for (auto line: lines)
 		{
@@ -232,18 +252,20 @@ void ElectricalData::initialize_circuit()
 		nodes[entry.second].genr.genr_mw = real(S);
 		nodes[entry.second].genr.genr_mvar = imag(S);
 		nodes[entry.second].genr.Pm0 = 0.0;
+		nodes[entry.second].genr.Rg = 0.0;
+		// Critically damp the generator
+		Complex X(nodes[entry.second].genr.Rg,nodes[entry.second].genr.Lg*2.0*pi*freq);
+		Complex Vgenr = Igenr*X+nodes[entry.second].V;
 		for (int i = 0; i < 3; i++)
 		{
-			double ig = abs(Igenr)*cos(arg(Igenr)+phase_angle[i]);
-			double vg = abs(Vnode[entry.second])*cos(arg(Vnode[entry.second])+phase_angle[i]);
+			double ig = nodes[entry.second].genr.i[i] = abs(Igenr)*cos(arg(Igenr)+phase_angle[i]);
+			double vg = abs(Vgenr)*cos(arg(Vgenr)+phase_angle[i]);
 			nodes[entry.second].genr.Pm0 += ig*vg;
 		}
-		nodes[entry.second].genr.Iamp = abs(Igenr);
-		nodes[entry.second].genr.Itheta = arg(Igenr);
-		configureGenerator(nodes[entry.second].genr);
+		nodes[entry.second].genr.Eamp = abs(Vgenr);
+		nodes[entry.second].genr.Etheta = arg(Vgenr);
 	}
 	delete Y;
 	delete [] Vsrc;
 	delete [] Vnode;
-	
 }

src/ElectricalData.h

@@ -60,8 +60,10 @@ class ElectricalData
 			/// These are the speed control parameters
 			double Tspd1, Tspd2, R, Agc; // Speed gain
 			/// Reference output current amplitude and phase
-			double Iamp, Itheta;
-			double D; // Friction for off nominal speed
+			double Eamp, Etheta, Lg, Rg;
+			double i[3];
+			/// Voltage control
+			double K, G;
 			/// Power for plant sizing information
 			double genr_mw, genr_mvar;
 			/// Initial speed

src/Makefile

@@ -31,7 +31,7 @@ test1: app
 	#./grid2ode ../data/test3.txt
 	#./grid2ode ../data/test2.txt
 	./m2c circuit.mo
-	${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp sim.cpp poly.cpp event.cpp 
+	${CXX} ${CFLAGS} ${INCLUDE} -o sim func.cpp euler.cpp event.cpp 
 
 app: ${OBJS}
 	${CXX} -o grid2ode ${CFLAGS} ${OBJS} ${LIBS}

src/event.cpp

@@ -5,11 +5,11 @@ void event(double* const dq, double* const q, double tNow)
 {
 	// Drop the the line from bus 2 to bus 3 in the IEEE 14 bus model
 	// Setting the line current to zero is the same as removing it.
-	q[ia_2_2_3] = 0.0;
-	q[ib_2_2_3] = 0.0;
-	q[ic_2_2_3] = 0.0;
-	dq[ia_2_2_3] = 0.0;
-	dq[ib_2_2_3] = 0.0;
-	dq[ic_2_2_3] = 0.0; 
+	q[ia_3_2_4] = 0.0;
+	q[ib_3_2_4] = 0.0;
+	q[ic_3_2_4] = 0.0;
+	dq[ia_3_2_4] = 0.0;
+	dq[ib_3_2_4] = 0.0;
+	dq[ic_3_2_4] = 0.0;
 }
 

src/m2c.lex

@@ -140,6 +140,7 @@ void number(const char* num)
 
 %option yylineno
 %option noyywrap
+%option nodefault
 
 ID [A-Za-z][A-Za-z_0-9]*
 NUMBER -?[0-9]*\.?[0-9]+([e][-+]?[0-9]+)?
@@ -152,7 +153,8 @@ OP [-+*/(),]
 {OP} {mathop(yytext);}
 {ID} {identifier(yytext);}
 {NUMBER} {number(yytext);}
-[.|\n] {}
+[.|\n] { }
+[ \t]+ { }
 
 %%
 
@@ -203,7 +205,7 @@ int main(int argc, char* argv[])
 	func_out << "#include \"vars.h\"" << std::endl;
 	func_out << "#include <iostream>" << std::endl;
 	func_out << "#include <cmath>" << std::endl;
-	func_out << "using namespace std;" << std::endl << std::endl;
+	func_out << "using namespace std;" << std::endl;
 	vars_out << "#ifndef _vars_h_" << std::endl;
 	vars_out << "#define _vars_h_" << std::endl << std::endl;
 	vars_out << "static const double pi = 3.141592653589793;" << std::endl << std::endl;
@@ -220,8 +222,8 @@ int main(int argc, char* argv[])
 			angle.push_back(plot_var(state_var,var_idx));
 		vars_out << "static const int " << state_var << " = " << (k++) << "; // " << (var_idx++) << std::endl;
 	}
-	vars_out << std::endl << "// Algebraic variables " << std::endl << std::endl;
-	func_out << std::endl << "const int numStateVars = " << state_vars.size() << ";" << std::endl << std::endl;
+	vars_out << std::endl << "// Algebraic variables " << std::endl;
+	func_out << std::endl << "const int numStateVars = " << state_vars.size() << ";" << std::endl;
 	k = 0;
 	for (auto alg_var: alg_vars)
 	{
@@ -235,7 +237,7 @@ int main(int argc, char* argv[])
 			volt_amp_est.push_back(plot_var(alg_var,var_idx));
 		vars_out << "static const int " << alg_var << " = " << (k++) << "; // " << (var_idx++) << std::endl;
 	}
-	func_out << std::endl << "const int numAlgVars = " << alg_vars.size() << ";" << std::endl << std::endl;
+	func_out << "const int numAlgVars = " << alg_vars.size() << ";" << std::endl << std::endl;
 
 	func_out << "void func(double* const dq, const double* const q, double* const alg, const double time)";
 	func_out << std::endl << "{" << std::endl;

src/main.cpp

@@ -17,8 +17,6 @@ static const char* const phases[3] = { "a", "b", "c"};
 // Set of voltage amplitude measurement variables
 static set<string> voltageAmplitudeEstimators;
 
-// Default rate of change for a transformer tap
-static const double TfRate = 1E-1;
 // Default initial values for voltages and currents
 
 // Initial value data for a state variable
@@ -45,7 +43,7 @@ set<state_var_init_t> stateVars;
 string number(double val)
 {
 	ostringstream sin;
-	sin << std::setprecision(20) << val;
+	sin << setprecision(10) << val;
 	return sin.str();
 }
 
@@ -89,83 +87,54 @@ string current(string flowName, string phase)
 	return string("i")+phase+string("_")+flowName;
 }
 
-// Transformer modulus at a bus and phase to a string
-string transformerModulus(string bus_name, string phase, ElectricalData::bus_data_t& bus)
-{
-	return string("tf_")+bus_name+"_"+phase+"_"+number(bus.ID);
-}
-
-// Line transformer modulus to a string
-string transformerModulus(string bus_name, string phase, ElectricalData::line_t& line)
-{
-	return string("tf_")+bus_name+"_"+phase+"_"+number(line.to)+"_"+number(line.from);
-}
-
 // Generator phase angle to a string
 string generatorAngle(string busName)
 {
 	return string("theta_")+busName;
 }
 
-// Generator output current
-string generatorCurrent(string busName, string phase)
+// Generator output voltage
+string generatorVoltage(string busName, string phase)
 {
-	return string("ig_")+busName+"_"+phase;
+	return string("vg_")+busName+"_"+phase;
 }
 
-// Generate functions to estimate the amplitude of a sinuisoid
-string sineAmplitudeMeasurement(string busName, string phase, string genName, string& estimate)
+// Generator output current
+string generatorCurrent(string busName, string phase)
 {
-	/**
-	 * Three point estimate from Shang-Teh Wu and Jyun-Lang Hong,
-	 * Five-Point Amplitude Estimation ofSinusoidal Signals: With Applicationto LVDT Signal Conditioning,
-	 * IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 59, NO. 3, MARCH 2010
-	 *
-	 * This is almost certainly good enough if we use 60 hertz as the signal frequency. We should never be
-	 * more than 1.5% away from that because the generators would be tripping offline otherwise.
-	 * For busses attached to the generator, we'll use the generator frequency.
-	 */
-	string freq = (genName.size() > 0) ? "("+number(2.0*pi*sys_freq_hz)+"+"+omega(genName)+")" : number(2.0*pi*sys_freq_hz);
-	estimate = "Amp_"+busName+"_"+phase;
-	if (voltageAmplitudeEstimators.find(estimate) != voltageAmplitudeEstimators.end())
-		return "";
-	voltageAmplitudeEstimators.insert(estimate);
-	string y1 = "Y_"+busName+"_"+phase+"_1";
-	string y2 = "Y_"+busName+"_"+phase+"_2";
-	string y3 = "Y_"+busName+"_"+phase+"_3";
-	string y1_gen = y1+"="+voltage(busName,phase)+";";
-	string y2_gen = y2+"=delay("+voltage(busName,phase)+",1.0/(3.0*"+freq+"/(2.0*pi)),time);";
-	string y3_gen = y3+"=delay("+voltage(busName,phase)+",2.0/(3.0*"+freq+"/(2.0*pi)),time);";
-	string Aest = estimate+"=(2.0/3.0)*sqrt(fabs("+y1+"*"+y1+"+"+y2+"*"+y2+"+"+y3+"*"+y3+"-"+y1+"*"+y2+"-"+y1+"*"+y3+"-"+y2+"*"+y3+"));";
-	return y1_gen+"\n"+y2_gen+"\n"+y3_gen+"\n"+Aest+"\n";
+	return string("ig_")+busName+"_"+phase;
 }
 
-string generatorExcitor(string busName, string phase)
+string generatorExcitor(string busName)
 {
-	return "E_"+busName+"_"+phase;
+	return "E_"+busName;
 }
 
-// Generator injected current
-string generatorSourceCurrent(string busName, string phase)
+// Generator source voltage
+string generatorSourceVoltage(string busName, string phase, ElectricalData::genr_t genr)
 {
+	string ef = generatorExcitor(busName);
 	string result;
 	string freq = number(2.0*pi*sys_freq_hz);
 	string w = omega(busName);
-	string theta = generatorAngle(busName);
+	string g = number((2.0*pi*sys_freq_hz)/genr.Eamp);
+	string theta = generatorAngle(busName)+"+"+number(genr.Etheta);;
 	if (phase == phases[0])
-		result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+"+theta+");\n";
 	else if (phase == phases[1])
-		result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+2.0*pi/3.0+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+2.0*pi/3.0+"+theta+");\n";
 	else
-		result += generatorCurrent(busName,phase)+"="+generatorExcitor(busName,phase)+"*cos(("+freq+"+"+w+")*time+4.0*pi/3.0+"+theta+");\n";
+		result += generatorVoltage(busName,phase)+"=("+number(genr.Eamp)+"+"+ef+")*cos(("+freq+"+"+w+")*time+4.0*pi/3.0+"+theta+");\n";
 	return result;
 }
 
 // Create dynamic model of a generator
 string createGenerator(string busName, ElectricalData::bus_data_t data)
 {
+	int i = 0;
 	string result;
 	string varName;
+	string Ef = generatorExcitor(busName);
 	string Hstr = number(2.0*data.genr.H);
 	string w = omega(busName);
 	string Pm = mech_power(busName);
@@ -176,38 +145,38 @@ string createGenerator(string busName, ElectricalData::bus_data_t data)
 	// Power demand at the terminal is sum of demand for each phase
 	for (auto phase: phases)
 	{
-		string estimate;
-		string amp = sineAmplitudeMeasurement(busName,phase,busName,estimate);
-		if (!amp.empty())
-			result += amp;
-		string Ie = generatorExcitor(busName,phase);
-		stateVars.insert(state_var_init_t(Ie,data.genr.Iamp));
-		result += "der("+Ie+")=10.0*("+number(data.vamp)+"-"+estimate+");\n";
-		result += generatorSourceCurrent(busName,phase);
+		stateVars.insert(state_var_init_t(generatorCurrent(busName,phase),data.genr.i[i++]));
+		result += generatorSourceVoltage(busName,phase,data.genr);
+		result += "der("+generatorCurrent(busName,phase)+")=("+
+			generatorVoltage(busName,phase)+"-"+voltage(busName,phase)+"-"+
+			generatorCurrent(busName,phase)+"*"+number(data.genr.Rg)+")/"+number(data.genr.Lg)+";\n";
 		if (powerDraw.empty())
-			powerDraw = voltage(busName,phase)+"*"+generatorCurrent(busName,phase);
+			powerDraw = generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
 		else 
-			powerDraw += "+"+voltage(busName,phase)+"*"+generatorCurrent(busName,phase);
+			powerDraw += "+"+generatorVoltage(busName,phase)+"*"+generatorCurrent(busName,phase);
 	}
 	result += Pe+"="+powerDraw+";\n";
-	stateVars.insert(state_var_init_t(generatorAngle(busName),data.genr.Itheta));
+	stateVars.insert(state_var_init_t(generatorAngle(busName),0.0));
 	stateVars.insert(state_var_init_t(w,data.genr.w0));
 	stateVars.insert(state_var_init_t(Pm,data.genr.Pm0));
 	stateVars.insert(state_var_init_t(C,data.genr.Pm0));
 	stateVars.insert(state_var_init_t(Pref,data.genr.Pm0));
+	stateVars.insert(state_var_init_t(Ef,0.0));
 	// Comment these out and see the next comment block to keep the generators still
+	result += "der("+Ef+")="+number(data.genr.G)+"*("+Pm+"-"+Pe+"-"+number(data.genr.K)+"*"+Ef+");\n";
 	result += "der("+generatorAngle(busName)+")="+w+";\n";
-	result += string("der(")+w+string(")=("+Pm+
-		"-"+Pe+"-"+number(data.genr.D)+"*"+w+")/")+Hstr+";\n";
+	result += "der("+w+")=("+Pm+"-"+Pe+")/"+Hstr+";\n";
 	result += "der("+C+")="+number(data.genr.Tspd1)+"*("+Pref+"-"+number(data.genr.R)+
 		"*"+w+"-"+C+");\n";
 	result += string("der(")+Pm+")="+
 		number(data.genr.Tspd2)+"*("+C+"-"+Pm+");\n";
 	result += "der("+Pref+")=-"+number(data.genr.Agc)+"*"+w+";\n";
-/*	result += "der("+generatorAngle(busName)+")=0;\n";
+/*	result += "der("+Ef+")=0.0;\n";
+	result += "der("+generatorAngle(busName)+")=0;\n";
 	result += "der("+w+")=0;\n";
 	result += "der("+Pm+")=0;\n";
-	result += "der("+C+")=0;\n";*/ 
+	result += "der("+C+")=0;\n"; 
+	result += "der("+Pref+")=0;\n"; */
 	return result;
 }
 
@@ -243,16 +212,12 @@ string createBus(
 			else
 				isum += "+"+current(line_name(line),phase);
 			i++;
-			if (line.tside == bus.ID)
-				isum += "*"+transformerModulus(busName,phase,line);
 		}
 		// Sum of outgoing line currents
 		for (auto line: data.getLines())
 		{
 			if (line.from != bus.ID) continue;
 			isum += string("-")+current(line_name(line),phase);
-			if (line.tside == bus.ID)
-				isum += "/"+transformerModulus(busName,phase,line);
 		}
 		// Inductor, so use current through the load
 		if (bus.L > 0.0 && isfinite(bus.L))
@@ -304,31 +269,6 @@ string createLine(ElectricalData::line_t& line)
 		string va = voltage(bus_name(line.from),phase);
 		string vb = voltage(bus_name(line.to),phase);
 		string i = current(line_name(line),phase);
-		// Append the transformer modulus to the appropriate side of the line
-		// and add the tap control equations. Tap changes to maintain the attached
-		// bus voltage
-		if (line.tside == line.to)
-		{
-			string estimate;
-			string amp = sineAmplitudeMeasurement(bus_name(line.to),phase,"",estimate);
-			if (!amp.empty())
-				result += amp;
-			varName = transformerModulus(bus_name(line.tside),phase,line);
-			stateVars.insert(state_var_init_t(varName,line.turns));
-			result += string("der(")+varName+")="+number(TfRate)+"*("+estimate+"-1.0);\n";
-			vb += "*"+transformerModulus(bus_name(line.tside),phase,line);
-		}
-		else if (line.tside == line.from)
-		{
-			string estimate;
-			string amp = sineAmplitudeMeasurement(bus_name(line.from),phase,"",estimate);
-			if (!amp.empty())
-				result += amp;
-			varName = transformerModulus(bus_name(line.tside),phase,line);
-			stateVars.insert(state_var_init_t(varName,line.turns));
-			result += string("der(")+varName+")="+number(TfRate)+"*(1.0-"+estimate+");\n";
-			va += "/"+transformerModulus(bus_name(line.tside),phase,line);
-		}
 		stateVars.insert(state_var_init_t(i,line.i[phase_num]));
 		result += string("der(")+i+string(")=(")+va+"-"+vb;
 		if (line.R > 0.0)

src/poly.cpp

-#include "poly.h"
-#include <cfloat>
-#include <cmath>
-
-void poly::init(double dx)
-{
-	this->dx[0] = this->dx[1] = dx;
-}
-
-void poly::append(double dx)
-{
-	this->dx[1] = this->dx[0];
-	this->dx[0] = dx;
-}
-
-double poly::extrapolate(double h)
-{
-	return 0.5*(3.0*dx[0]-dx[1])*h;
-}
-
-double poly::time_until(double q)
-{
-	double ddx = (3.0*dx[0]-dx[1])*0.5;
-	double dt = (ddx != 0.0) ? q/fabs(ddx) : DBL_MAX;
-	return dt;
-}

src/poly.h

-#ifndef _poly_h_
-#define _poly_h_
-
-class poly
-{
-	public:
-		// Assign first estimate of derivative
-		void init(double dx);
-		// Append a new data value h units of time following the previous one
-		void append(double dx);
-		// Time to change by q
-		double time_until(double q);
-		// Change over given interval
-		double extrapolate(double h);
-	private:
-		// Two derivative approximations
-		double dx[2];
-};
-
-#endif

src/sim.cpp

-#include "sim.h"
-#include "poly.h"
-#include <cfloat>
-#include <algorithm>
-#include <vector>
-#include <list>
-#include <cassert>
-#include <iostream>
-#include <cmath>
-#include <cstring>
-using namespace std;
-
-// Current simulation time
-static double tNow = 0.0;
-// Integration quantum
-static double quant = 0.0;
-// Simulation ending time
-static double tEnd = 0.0;
-// Reporting interval
-static double cint = 0.0;
-// Delay sampling interval. Sub-milliseconds should be plenty.
-static const double dint = 0.0001;
-// Delay data stored in a circular queue
-static int dHead = 0;
-// One seconds of delay data should be more than enough
-// Delays are only used to obtain sine amplitude measurements
-static const double maxDelay = 1.0;
-static const int dSize = 1+int(maxDelay/dint);
-static double **dQueue;
-// Model variable values
-static double *q, *dq, *alg;
-static poly* p;
-
-void print()
-{
-	cout << tNow;
-	for (int i = 0; i < numStateVars; i++)
-		cout << " " << q[i];
-	for (int i = 0; i < numAlgVars; i++)
-		cout << " " << alg[i];
-	cout << endl;
-}
-
-double delay(int var, double h, double t)
-{
-	double f;
-	// If this happens, just return the oldest value. It
-	// shouldn't matter becuase this should only occur
-	// early in the simulation.
-	assert(h > 0.0);
-	if (h > maxDelay)
-	{
-		return dQueue[var][(dHead+1)%dSize];
-	}
-	// time in bin
-	double dt = t-dint*floor(t/dint);
-	// First bin below the time we want
-	int below = (int)(h/dint);
-	if (below == 0)
-	{
-		assert(dt >= h);
-		f = (dt-h)/dt;
-		return (1.0-f)*dQueue[var][dHead]+f*q[var];
-	}
-	// First bin above the time we want
-	int above = below+1;
-	f = dt/dint;
-	below = dHead-below;
-	above = dHead-above;
-	if (below < 0) below = dSize+below;
-	if (above < 0) above = dSize+above;
-	return (1.0-f)*dQueue[var][below]+f*dQueue[var][above];
-}
-
-int main(int argc, char** argv)
-{
-	int reports = 1, delays = 1;
-	double h, hReport, hEnd, hDelay;
-	if (argc != 7)
-	{
-		cout << "Must supply quantum, cint, and tend" << endl;
-		return 0;
-	}
-	for (int i = 1; i < argc; i++)
-	{
-		if (strcmp("quantum",argv[i]) == 0)
-			quant = atof(argv[++i]);
-		else if (strcmp("cint",argv[i]) == 0)
-			cint = atof(argv[++i]);
-		else if (strcmp("tend",argv[i]) == 0)
-			tEnd = atof(argv[++i]);
-	}
-	alg = new double[numAlgVars];
-	dq = new double[numStateVars];