Added KDE3 version of ktechlab

git-svn-id: svn://anonsvn.kde.org/home/kde/branches/trinity/applications/ktechlab@1095338 283d02a7-25f6-0310-bc7c-ecb5cbfe19da

49 files changed, 5756 insertions, 0 deletions

diff --git a/src/electronics/simulation/Makefile.am b/src/electronics/simulation/Makefile.am
new file mode 100644
index 0000000..c45c6a0
--- /dev/null
+++ b/src/electronics/simulation/Makefile.am
@@ -0,0 +1,11 @@
+INCLUDES = -I$(top_srcdir)/src -I$(top_srcdir)/src/electronics $(all_includes)
+METASOURCES = AUTO
+noinst_LTLIBRARIES =  libelements.la
+libelements_la_SOURCES = cccs.cpp ccvs.cpp circuit.cpp currentsource.cpp \
+			diode.cpp element.cpp elementset.cpp logic.cpp matrix.cpp vccs.cpp vcvs.cpp \
+			voltagesource.cpp capacitance.cpp resistance.cpp currentsignal.cpp voltagepoint.cpp \
+			voltagesignal.cpp elementsignal.cpp nonlinear.cpp reactive.cpp vec.cpp bjt.cpp opamp.cpp \
+	inductance.cpp
+noinst_HEADERS = cccs.h ccvs.h circuit.h currentsource.h diode.h element.h \
+			elementset.h logic.h matrix.h vccs.h vcvs.h voltagesource.h capacitance.h \
+			resistance.h elementsignal.h nonlinear.h reactive.h vec.h bjt.h opamp.h inductance.h
diff --git a/src/electronics/simulation/bjt.cpp b/src/electronics/simulation/bjt.cpp
new file mode 100644
index 0000000..8daab40
--- /dev/null
+++ b/src/electronics/simulation/bjt.cpp
@@ -0,0 +1,257 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "bjt.h"
+#include "diode.h"
+#include "elementset.h"
+#include "matrix.h"
+
+#include <cmath>
+using namespace std;
+
+
+//BEGIN class BJTSettings
+BJTSettings::BJTSettings()
+{
+	I_S = 1e-16;
+	N_F = 1.0;
+	N_R = 1.0;
+	B_F = 100.0;
+	B_R = 1.0;
+}
+//END class BJTSettings
+
+
+
+//BEGIN class BJTState
+BJTState::BJTState()
+{
+	reset();
+}
+
+
+void BJTState::reset()
+{
+	for ( unsigned i = 0; i < 3; ++i )
+	{
+		for ( unsigned j = 0; j < 3; ++j )
+			A[i][j] = 0.0;
+			
+		I[i] = 0.0;
+	}
+}
+
+
+BJTState BJTState::operator-( const BJTState & s ) const
+{
+	BJTState newState( *this );
+	
+	for ( unsigned i = 0; i < 3; ++i )
+	{
+		for ( unsigned j = 0; j < 3; ++j )
+			newState.A[i][j] -= s.A[i][j];
+			
+		newState.I[i] -= s.I[i];
+	}
+	
+	return newState;
+}
+//END class BJTState
+
+
+
+//BEGIN class BJT
+BJT::BJT( const bool isNPN )
+	: NonLinear()
+{
+	V_BE_prev = 0.0;
+	V_BC_prev = 0.0;
+	m_pol = isNPN ? 1 : -1;
+	m_numCNodes = 3;
+}
+
+
+BJT::~BJT()
+{
+}
+
+
+void BJT::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[1]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+	}
+	if ( !p_cnode[2]->isGround )
+	{
+		p_A->setUse( p_cnode[2]->n(), p_cnode[2]->n(), Map::et_unstable, false );
+	}
+	
+	if ( !p_cnode[0]->isGround && !p_cnode[2]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[2]->n(), Map::et_unstable, false );
+		p_A->setUse( p_cnode[2]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+	if ( !p_cnode[1]->isGround && !p_cnode[2]->isGround )
+	{
+		p_A->setUse( p_cnode[2]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+		p_A->setUse( p_cnode[1]->n(), p_cnode[2]->n(), Map::et_unstable, false );
+	}
+}
+
+
+void BJT::add_initial_dc()
+{
+	V_BE_prev = 0.0;
+	V_BC_prev = 0.0;
+	m_os.reset();
+	update_dc();
+}
+
+
+void BJT::updateCurrents()
+{
+	if (!b_status)
+		return;
+	
+	double V_B = p_cnode[0]->v;
+	double V_C = p_cnode[1]->v;
+	double V_E = p_cnode[2]->v;
+	
+	double V_BE = (V_B - V_E) * m_pol;
+	double V_BC = (V_B - V_C) * m_pol;
+	
+	double I_BE, I_BC, I_T, g_BE, g_BC, g_IF, g_IR;
+	calcIg( V_BE, V_BC, & I_BE, & I_BC, & I_T, & g_BE, & g_BC, & g_IF, & g_IR );
+	
+	m_cnodeI[1] = I_BC - I_T;
+	m_cnodeI[2] = I_BE + I_T;
+	m_cnodeI[0] = -(m_cnodeI[1] + m_cnodeI[2]);
+}
+
+
+void BJT::update_dc()
+{
+	if (!b_status)
+		return;
+	
+	calc_eq();
+	
+	BJTState diff = m_ns - m_os;
+	for ( unsigned i = 0; i < 3; ++i )
+	{
+		for ( unsigned j = 0 ; j < 3; ++j )
+			A_g( i, j ) += diff.A[i][j];
+		
+		b_i( i ) += diff.I[i];
+	}
+	
+	m_os = m_ns;
+}
+
+
+void BJT::calc_eq()
+{
+	double V_B = p_cnode[0]->v;
+	double V_C = p_cnode[1]->v;
+	double V_E = p_cnode[2]->v;
+	
+	double V_BE = (V_B - V_E) * m_pol;
+	double V_BC = (V_B - V_C) * m_pol;
+	
+	double I_S = m_bjtSettings.I_S;
+	double N_F = m_bjtSettings.N_F;
+	double N_R = m_bjtSettings.N_R;
+	
+	// adjust voltage to help convergence
+	double V_BEcrit = diodeCriticalVoltage( I_S, N_F * V_T );
+	double V_BCcrit = diodeCriticalVoltage( I_S, N_R * V_T );
+	V_BE_prev = V_BE = diodeVoltage( V_BE, V_BE_prev, V_T * N_F, V_BEcrit );
+	V_BC_prev = V_BC = diodeVoltage( V_BC, V_BC_prev, V_T * N_R, V_BCcrit );
+	
+	double I_BE, I_BC, I_T, g_BE, g_BC, g_IF, g_IR;
+	calcIg( V_BE, V_BC, & I_BE, & I_BC, & I_T, & g_BE, & g_BC, & g_IF, & g_IR );
+	
+	double I_eq_B = I_BE - V_BE * g_BE;
+	double I_eq_C = I_BC - V_BC * g_BC;
+	double I_eq_E = I_T - V_BE * g_IF + V_BC * g_IR;
+	
+	m_ns.A[0][0] = g_BC + g_BE;
+	m_ns.A[0][1] = -g_BC;
+	m_ns.A[0][2] = -g_BE;
+	
+	m_ns.A[1][0] = -g_BC + (g_IF - g_IR);
+	m_ns.A[1][1] = g_IR + g_BC;
+	m_ns.A[1][2] = -g_IF;
+	
+	m_ns.A[2][0] = -g_BE - (g_IF - g_IR);
+	m_ns.A[2][1] = -g_IR;
+	m_ns.A[2][2] = g_BE + g_IF;
+	
+	m_ns.I[0] = (-I_eq_B - I_eq_C) * m_pol;
+	m_ns.I[1] = (+I_eq_C - I_eq_E) * m_pol;
+	m_ns.I[2] = (+I_eq_B + I_eq_E) * m_pol;
+}
+
+
+void BJT::calcIg( double V_BE, double V_BC, double * I_BE, double * I_BC, double * I_T, double * g_BE, double * g_BC, double * g_IF, double * g_IR )
+{
+	double I_S = m_bjtSettings.I_S;
+	double N_F = m_bjtSettings.N_F;
+	double N_R = m_bjtSettings.N_R;
+	double B_F = m_bjtSettings.B_F;
+	double B_R = m_bjtSettings.B_R;
+	
+	// base-emitter diodes
+	double g_tiny = (V_BE < (-10 * V_T * N_F)) ? I_S : 0;
+  
+	double I_F;
+	diodeJunction( V_BE, I_S, V_T * N_F, & I_F, g_IF );
+	
+	double I_BEI = I_F / B_F;
+	double g_BEI = *g_IF / B_F;
+	double I_BEN = g_tiny * V_BE;
+	double g_BEN = g_tiny;
+	*I_BE = I_BEI + I_BEN;
+	*g_BE = g_BEI + g_BEN;
+	
+	// base-collector diodes
+	g_tiny = (V_BC < (-10 * V_T * N_R)) ? I_S : 0;
+	
+	double I_R;
+	diodeJunction( V_BC, I_S, V_T * N_R, & I_R, g_IR );
+	
+	double I_BCI = I_R / B_R;
+	double g_BCI = *g_IR / B_R;
+	double I_BCN = g_tiny * V_BC;
+	double g_BCN = g_tiny;
+	*I_BC = I_BCI + I_BCN;
+	*g_BC = g_BCI + g_BCN;
+	
+	*I_T = I_F - I_R;
+}
+
+
+void BJT::setBJTSettings( const BJTSettings & settings )
+{
+	m_bjtSettings = settings;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+}
+//END class BJT
+
+
diff --git a/src/electronics/simulation/bjt.h b/src/electronics/simulation/bjt.h
new file mode 100644
index 0000000..52022aa
--- /dev/null
+++ b/src/electronics/simulation/bjt.h
@@ -0,0 +1,75 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef BJT_H
+#define BJT_H
+
+#include "matrix.h"
+#include "nonlinear.h"
+
+class BJTState
+{
+	public:
+		BJTState();
+		void reset();
+		
+		BJTState operator-( const BJTState & s ) const;
+	
+		double A[3][3];
+		double I[3];
+};
+
+
+class BJTSettings
+{
+	public:
+		BJTSettings();
+		
+		double I_S; ///< saturation current
+		double N_F; ///< forward emission coefficient
+		double N_R; ///< reverse emission coefficient
+		double B_F; ///< forward beta
+		double B_R; ///< reverse beta
+};
+
+
+/**
+@author David Saxton
+*/
+class BJT : public NonLinear
+{
+	public:
+		BJT( bool isNPN );
+		virtual ~BJT();
+	
+		virtual Type type() const { return Element_BJT; }
+		virtual void update_dc();
+		virtual void add_initial_dc();
+		virtual void add_map();
+		BJTSettings settings() const { return m_bjtSettings; }
+		void setBJTSettings( const BJTSettings & settings );
+	
+	protected:
+		virtual void updateCurrents();
+		/**
+		 * Calculates the new BJTState from the voltages on the nodes.
+		 */
+		void calc_eq();
+		
+		void calcIg( double V_BE, double V_BC, double * I_BE, double * I_BC, double * I_T, double * g_BE, double * g_BC, double * g_IF, double * g_IR );
+
+		BJTState m_os;
+		BJTState m_ns;
+		int m_pol;
+		double V_BE_prev, V_BC_prev;
+		BJTSettings m_bjtSettings;
+};
+
+#endif
diff --git a/src/electronics/simulation/capacitance.cpp b/src/electronics/simulation/capacitance.cpp
new file mode 100644
index 0000000..9087c7f
--- /dev/null
+++ b/src/electronics/simulation/capacitance.cpp
@@ -0,0 +1,117 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "capacitance.h"
+#include "matrix.h"
+
+Capacitance::Capacitance( const double capacitance, const double delta )
+	: Reactive(delta)
+{
+	m_cap = capacitance;
+	g_eq_old = i_eq_old = 0.;
+	m_numCNodes = 2;
+	setMethod( Capacitance::m_euler );
+}
+
+Capacitance::~Capacitance()
+{
+}
+
+void Capacitance::setCapacitance( const double c )
+{
+	m_cap = c;
+}
+
+void Capacitance::add_initial_dc()
+{
+	// We don't need to do anything here, as time_step() will do that for us,
+	// apart from to make sure our old values are 0
+	g_eq_old = i_eq_old = 0.;
+}
+
+void Capacitance::updateCurrents()
+{
+	if (!b_status) return;
+	const double r_i = (p_cnode[0]->v-p_cnode[1]->v)*g_eq_old;
+	m_cnodeI[0] = -i_eq_old-r_i;
+	m_cnodeI[1] = -m_cnodeI[0];
+}
+
+
+void Capacitance::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[1]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+	}
+	
+	if ( !p_cnode[0]->isGround && !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+		p_A->setUse( p_cnode[1]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+}
+
+
+void Capacitance::time_step()
+{
+	if (!b_status) return;
+	
+	double v = p_cnode[0]->v-p_cnode[1]->v;
+	double i_eq_new = 0.0, g_eq_new = 0.0;
+	
+	if ( m_method == Capacitance::m_euler )
+	{
+		g_eq_new = m_cap/m_delta;
+		i_eq_new = -v*g_eq_new;
+	}
+	else if ( m_method == Capacitance::m_trap ) {
+		// TODO Implement + test trapezoidal method
+		g_eq_new = 2.*m_cap/m_delta;
+	}
+	
+	if ( g_eq_old != g_eq_new )
+	{
+		A_g( 0, 0 ) += g_eq_new-g_eq_old;
+		A_g( 1, 1 ) += g_eq_new-g_eq_old;
+		A_g( 0, 1 ) -= g_eq_new-g_eq_old;
+		A_g( 1, 0 ) -= g_eq_new-g_eq_old;
+	}
+	
+	if ( i_eq_new != i_eq_old )
+	{
+		b_i( 0 ) -= i_eq_new-i_eq_old;
+		b_i( 1 ) += i_eq_new-i_eq_old;
+	}
+	
+	g_eq_old = g_eq_new;
+	i_eq_old = i_eq_new;
+}
+
+bool Capacitance::updateStatus()
+{
+	b_status = Reactive::updateStatus();
+	if ( m_method == Capacitance::m_none ) b_status = false;
+	return b_status;
+}
+
+void Capacitance::setMethod( Method m )
+{
+	m_method = m;
+	updateStatus();
+}
+
+
diff --git a/src/electronics/simulation/capacitance.h b/src/electronics/simulation/capacitance.h
new file mode 100644
index 0000000..ccc083d
--- /dev/null
+++ b/src/electronics/simulation/capacitance.h
@@ -0,0 +1,55 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CAPACITANCE_H
+#define CAPACITANCE_H
+
+#include "reactive.h"
+
+/**
+@author David Saxton
+@short Capacitance
+*/
+class Capacitance : public Reactive
+{
+public:
+	enum Method
+	{
+		m_none, // None
+		m_euler, // Backward Euler
+		m_trap // Trapezoidal (currently unimplemented)
+	};
+	Capacitance( const double capacitance, const double delta );
+	virtual ~Capacitance();
+	
+	virtual Type type() const { return Element_Capacitance; }
+	/**
+	 * Set the stepping use for numerical integration of capacitance,
+	 * and the interval between successive updates
+	 */
+	void setMethod( Method m );
+	virtual void time_step();
+	virtual void add_initial_dc();
+	void setCapacitance( const double c );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual bool updateStatus();
+	
+private:
+	double m_cap; // Capacitance
+	Method m_method; // Method of integration
+	
+	double g_eq_old;
+	double i_eq_old;
+};
+
+#endif
diff --git a/src/electronics/simulation/cccs.cpp b/src/electronics/simulation/cccs.cpp
new file mode 100644
index 0000000..9725735
--- /dev/null
+++ b/src/electronics/simulation/cccs.cpp
@@ -0,0 +1,89 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "cccs.h"
+
+CCCS::CCCS( const double gain )
+	: Element::Element()
+{
+	m_g = gain;
+	m_numCBranches = 2;
+	m_numCNodes = 4;
+}
+
+CCCS::~CCCS()
+{
+}
+
+void CCCS::setGain( const double g )
+{
+	if ( m_g == g )
+		return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	m_g = g;
+	add_initial_dc();
+}
+
+
+void CCCS::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[2]->isGround )
+	{
+		p_A->setUse_b( p_cnode[2]->n(), p_cbranch[1]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[3]->isGround )
+	{
+		p_A->setUse_b( p_cnode[3]->n(), p_cbranch[1]->n(), Map::et_constant, true );
+	}
+	p_A->setUse_d( p_cbranch[1]->n(), p_cbranch[0]->n(), Map::et_stable, true );
+	p_A->setUse_d( p_cbranch[1]->n(), p_cbranch[1]->n(), Map::et_constant, true );
+}
+
+void CCCS::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_b( 0, 0 ) = 1;
+	A_c( 0, 0 ) = 1;
+	A_b( 1, 0 ) = -1;
+	A_c( 0, 1 ) = -1;
+	A_b( 2, 1 ) = 1;
+	A_b( 3, 1 ) = -1;
+	A_d( 1, 0 ) = -m_g;
+	A_d( 1, 1 ) = 1;
+}
+
+void CCCS::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[1] = p_cbranch[0]->i;
+	m_cnodeI[0] = -m_cnodeI[1];
+	m_cnodeI[3] = p_cbranch[1]->i;
+	m_cnodeI[2] = -m_cnodeI[3];
+}
+
diff --git a/src/electronics/simulation/cccs.h b/src/electronics/simulation/cccs.h
new file mode 100644
index 0000000..ba86e9e
--- /dev/null
+++ b/src/electronics/simulation/cccs.h
@@ -0,0 +1,41 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CCCS_H
+#define CCCS_H
+
+#include "element.h"
+
+/**
+CNodes n0 and n1 are used for the current control.
+CNodes n2 and n3 are used for the current output.
+Branches b0 and b1 are for control and output
+@short Current Controlled Current Source
+@author David Saxton
+*/
+class CCCS : public Element
+{
+public:
+	CCCS( const double gain );
+	virtual ~CCCS();
+	
+	virtual Type type() const { return Element_CCCS; }
+	void setGain( const double g );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_g; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/ccvs.cpp b/src/electronics/simulation/ccvs.cpp
new file mode 100644
index 0000000..fc12bf5
--- /dev/null
+++ b/src/electronics/simulation/ccvs.cpp
@@ -0,0 +1,91 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "ccvs.h"
+
+CCVS::CCVS( const double gain )
+	: Element::Element()
+{
+	m_g = gain;
+	m_numCBranches = 2;
+	m_numCNodes = 4;
+}
+
+CCVS::~CCVS()
+{
+}
+
+void CCVS::setGain( const double g )
+{
+	if ( m_g == g )
+		return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	m_g = g;
+	add_initial_dc();
+}
+
+
+void CCVS::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[2]->isGround )
+	{
+		p_A->setUse_b( p_cnode[2]->n(), p_cbranch[1]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[1]->n(), p_cnode[2]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[3]->isGround )
+	{
+		p_A->setUse_b( p_cnode[3]->n(), p_cbranch[1]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[1]->n(), p_cnode[3]->n(), Map::et_constant, true );
+	}
+	p_A->setUse_d( p_cbranch[1]->n(), p_cbranch[0]->n(), Map::et_stable, true );
+}
+
+
+void CCVS::add_initial_dc()
+{
+	if (!b_status) return;
+	
+	A_b( 0, 0 ) = 1;
+	A_c( 0, 0 ) = 1;
+	A_b( 1, 0 ) = -1;
+	A_c( 0, 1 ) = -1;
+	A_b( 2, 1 ) = 1;
+	A_c( 1, 2 ) = 1;
+	A_b( 3, 1 ) = -1;
+	A_c( 1, 3 ) = -1;
+	A_d( 1, 0 ) = -m_g;
+}
+
+void CCVS::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[1] = p_cbranch[0]->i;
+	m_cnodeI[0] = -m_cnodeI[1];
+	m_cnodeI[3] = p_cbranch[0]->i;
+	m_cnodeI[2] = -m_cnodeI[3];
+}
+
diff --git a/src/electronics/simulation/ccvs.h b/src/electronics/simulation/ccvs.h
new file mode 100644
index 0000000..bcb1ac0
--- /dev/null
+++ b/src/electronics/simulation/ccvs.h
@@ -0,0 +1,41 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CCVS_H
+#define CCVS_H
+
+#include "element.h"
+
+/**
+CNodes n0 and n1 are used for the current control.
+CNodes n2 and n3 are used for the voltage output.
+Branches b0 and b1 are for control and output
+@short Current Controlled Voltage Source
+@author David Saxton
+*/
+class CCVS : public Element
+{
+public:
+	CCVS( const double gain );
+	virtual ~CCVS();
+	
+	virtual Type type() const { return Element_CCVS; }
+	void setGain( const double g );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_g; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/circuit.cpp b/src/electronics/simulation/circuit.cpp
new file mode 100644
index 0000000..c152756
--- /dev/null
+++ b/src/electronics/simulation/circuit.cpp
@@ -0,0 +1,550 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include <vector>
+#include "circuit.h"
+#include "circuitdocument.h"
+#include "element.h"
+#include "elementset.h"
+#include "logic.h"
+#include "matrix.h"
+#include "nonlinear.h"
+#include "pin.h"
+#include "reactive.h"
+#include "wire.h"
+
+
+#include <cmath>
+#include <map>
+
+typedef std::multimap<int, PinList> PinListMap;
+
+//BEGIN class Circuit
+Circuit::Circuit()
+{
+	m_bCanAddChanged = true;
+	m_pNextChanged[0] = m_pNextChanged[1] = 0l;
+	m_logicOutCount = 0;
+	m_bCanCache = false;
+	m_pLogicOut = 0l;
+	m_elementSet = new ElementSet( this, 0, 0 );
+	m_cnodeCount = m_branchCount = -1;
+	m_prepNLCount = 0;
+	m_pLogicCacheBase = new LogicCacheNode;
+}
+
+Circuit::~Circuit()
+{
+	delete m_elementSet;
+	delete m_pLogicCacheBase;
+	delete[] m_pLogicOut;
+}
+
+
+void Circuit::addPin( Pin *node )
+{
+	if ( m_pinList.contains(node) ) return;
+	m_pinList.append(node);
+}
+
+void Circuit::addElement( Element *element )
+{
+	if ( m_elementList.contains(element) ) return;
+	m_elementList.append(element);
+}
+
+bool Circuit::contains( Pin *node )
+{
+	return m_pinList.contains(node);
+}
+
+
+// static function
+int Circuit::identifyGround( PinList nodeList, int *highest )
+{
+	// What this function does:
+	// We are given a list of pins. First, we divide them into groups of pins
+	// that are directly connected to each other (e.g. through wires or
+	// switches). Then, each group of connected pins is looked at to find the
+	// pin with the highest "ground priority", and this is taken to be
+	// the priority of the group. The highest ground priority from all the
+	// groups is recorded. If the highest ground priority found is the maximum,
+	// then all the pins in groups with this priority are marked as ground
+	// (their eq-id is set to -1). Otherwise, the first group of pins with the
+	// highest ground priority found is marked as ground, and all others are
+	// marked as non ground (their eq-id is set to 0).
+	
+	int temp_highest;
+	if (!highest)
+		highest = &temp_highest;
+	
+	// Now to give all the Pins ids
+	PinListMap eqs;
+	while ( !nodeList.isEmpty() )
+	{
+		PinList associated;
+		PinList nodes;
+		Pin *node = *nodeList.begin();
+		recursivePinAdd( node, &nodeList, &associated, &nodes );
+		if ( nodes.size() > 0 )
+		{
+			eqs.insert( std::make_pair( associated.size(), nodes ) );
+		}
+	}
+	
+	
+	// Now, we want to look through the associated Pins, 
+	// to find the ones with the highest "Ground Priority". Anything with a lower
+	// priority than Pin::gt_never will not be considered
+	*highest = Pin::gt_never; // The highest priority found so far
+	int numGround = 0; // The number of node groups found with that priority
+	const PinListMap::iterator eqsEnd = eqs.end();
+	for ( PinListMap::iterator it = eqs.begin(); it != eqsEnd; ++it )
+	{
+		int highPri = Pin::gt_never; // The highest priority found in these group of nodes
+		const PinList::iterator send = it->second.end();
+		for ( PinList::iterator sit = it->second.begin(); sit != send; ++sit )
+		{
+			if ( (*sit)->groundType() < highPri )
+				highPri = (*sit)->groundType();
+		}
+		
+		if ( highPri == *highest )
+			numGround++;
+		
+		else if ( highPri < *highest ) 
+		{
+			numGround = 1;
+			*highest = highPri;
+		}
+	}
+	
+	if ( *highest == Pin::gt_never )
+	{
+		(*highest)--;
+		numGround=0;
+	}
+	// If there are no Always Ground nodes, then we only want to set one of the nodes as ground
+	else if ( *highest > Pin::gt_always )
+		numGround = 1;
+	
+	
+	// Now, we can give the nodes their cnode ids, or tell them they are ground
+	bool foundGround = false; // This is only used when we don't have a Always ground node
+	for ( PinListMap::iterator it = eqs.begin(); it != eqsEnd; ++it )
+	{
+		bool ground = false;
+		const PinList::iterator send = it->second.end();
+		for ( PinList::iterator sit = it->second.begin(); sit != send; ++sit )
+		{
+			ground |= (*sit)->groundType() <= (*highest);
+		}
+		if ( ground && (!foundGround || *highest == Pin::gt_always ) )
+		{
+			for ( PinList::iterator sit = it->second.begin(); sit != send; ++sit )
+			{
+				(*sit)->setEqId(-1);
+			}
+			foundGround = true;
+		}
+		else
+		{
+			for ( PinList::iterator sit = it->second.begin(); sit != send; ++sit )
+			{
+				(*sit)->setEqId(0);
+			}
+		}
+	}
+	
+	return numGround;
+}
+
+
+void Circuit::init()
+{
+	m_branchCount = 0;
+	
+	const ElementList::iterator listEnd = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != listEnd; ++it )
+	{
+		m_branchCount += (*it)->numCBranches();
+	}
+	
+	// Now to give all the Pins ids
+	int groundCount = 0;
+	PinListMap eqs;
+	PinList unassignedNodes = m_pinList;
+	while ( !unassignedNodes.isEmpty() )
+	{
+		PinList associated;
+		PinList nodes;
+		Pin *node = *unassignedNodes.begin();
+		if ( recursivePinAdd( node, &unassignedNodes, &associated, &nodes ) ) {
+			groundCount++;
+		}
+		if ( nodes.size() > 0 ) {
+			eqs.insert( std::make_pair( associated.size(), nodes ) );
+		}
+	}
+	
+	m_cnodeCount = eqs.size() - groundCount;
+	
+	delete m_pLogicCacheBase;
+	m_pLogicCacheBase = 0l;
+	
+	delete m_elementSet;
+	m_elementSet = new ElementSet( this, m_cnodeCount, m_branchCount );
+	
+	// Now, we can give the nodes their cnode ids, or tell them they are ground
+	Vector *x = m_elementSet->x();
+	int i=0;
+	const PinListMap::iterator eqsEnd = eqs.end();
+	for ( PinListMap::iterator it = eqs.begin(); it != eqsEnd; ++it )
+	{
+		bool foundGround = false;
+		
+		const PinList::iterator sEnd = it->second.end();
+		for ( PinList::iterator sit = it->second.begin(); sit != sEnd; ++sit )
+			foundGround |= (*sit)->eqId() == -1;
+			
+		if ( foundGround )
+			continue;
+		
+		bool foundEnergyStoragePin = false;
+		
+		for ( PinList::iterator sit = it->second.begin(); sit != sEnd; ++sit )
+		{
+			(*sit)->setEqId(i);
+			
+			bool energyStorage = false;
+			const ElementList elements = (*sit)->elements();
+			ElementList::const_iterator elementsEnd = elements.end();
+			for ( ElementList::const_iterator it = elements.begin(); it != elementsEnd; ++it )
+			{
+				if ( !*it )
+					continue;
+				
+				if ( ((*it)->type() == Element::Element_Capacitance)
+									|| ((*it)->type() == Element::Element_Inductance) )
+				{
+					energyStorage = true;
+					break;
+				}
+			}
+			
+			// A pin attached to an energy storage pin overrides one that doesn't.
+			// If the two pins have equal status with in this regard, we pick the
+			// one with the highest absolute voltage on it.
+			
+			if ( foundEnergyStoragePin && !energyStorage )
+				continue;
+			
+			double v = (*sit)->voltage();
+			
+			if ( energyStorage && !foundEnergyStoragePin )
+			{
+				foundEnergyStoragePin = true;
+				(*x)[i] = v;
+				continue;
+			}
+			
+			if ( std::abs(v) > std::abs( (*x)[i] ) )
+				(*x)[i] = v;
+		}
+		i++;
+	}
+	
+	
+	// And add the elements to the elementSet
+	for ( ElementList::iterator it = m_elementList.begin(); it != listEnd; ++it )
+	{
+		// We don't want the element to prematurely try to do anything,
+		// as it doesn't know its actual cnode ids yet
+		(*it)->setCNodes();
+		(*it)->setCBranches();
+		m_elementSet->addElement(*it);
+	}	
+	// And give the branch ids to the elements
+	i=0;
+	for ( ElementList::iterator it = m_elementList.begin(); it != listEnd; ++it )
+	{
+		switch ( (*it)->numCBranches() )
+		{
+			case 0:
+				break;
+			case 1:
+				(*it)->setCBranches( i );
+				i += 1;
+				break;
+			case 2:
+				(*it)->setCBranches( i, i+1 );
+				i += 2;
+				break;
+			case 3:
+				(*it)->setCBranches( i, i+1, i+2 );
+				i += 3;
+				break;
+			default:
+				// What the?!
+				break;
+		}
+	}
+}
+
+
+void Circuit::initCache()
+{
+	m_elementSet->updateInfo();
+	
+	m_bCanCache = true;
+	m_logicOutCount = 0;
+	
+	delete[] m_pLogicOut;
+	m_pLogicOut = 0l;
+	
+	delete m_pLogicCacheBase;
+	m_pLogicCacheBase = 0l;
+	
+	const ElementList::iterator end = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != end && m_bCanCache; ++it )
+	{
+		switch ( (*it)->type() )
+		{
+			case Element::Element_BJT:
+			case Element::Element_CCCS:
+			case Element::Element_CCVS:
+			case Element::Element_CurrentSource:
+			case Element::Element_Diode:
+			case Element::Element_LogicIn:
+			case Element::Element_OpAmp:
+			case Element::Element_Resistance:
+			case Element::Element_VCCS:
+			case Element::Element_VCVS:
+			case Element::Element_VoltagePoint:
+			case Element::Element_VoltageSource:
+			{
+				break;
+			}
+				
+			case Element::Element_LogicOut:
+			{
+				m_logicOutCount++;
+				break;
+			}
+				
+			case Element::Element_CurrentSignal:
+			case Element::Element_VoltageSignal:
+			case Element::Element_Capacitance:
+			case Element::Element_Inductance:
+			{
+				m_bCanCache = false;
+				break;
+			}
+		}
+	}
+	
+	if ( !m_bCanCache )
+		return;
+	
+	m_pLogicOut = new LogicOut*[m_logicOutCount];
+	unsigned i = 0;
+	for ( ElementList::iterator it = m_elementList.begin(); it != end && m_bCanCache; ++it )
+	{
+		if ( (*it)->type() == Element::Element_LogicOut )
+			m_pLogicOut[i++] = static_cast<LogicOut*>(*it);
+	}
+	
+	m_pLogicCacheBase = new LogicCacheNode;
+}
+
+
+void Circuit::setCacheInvalidated()
+{
+	if (m_pLogicCacheBase)
+	{
+		delete m_pLogicCacheBase->high;
+		m_pLogicCacheBase->high = 0l;
+	
+		delete m_pLogicCacheBase->low;
+		m_pLogicCacheBase->low = 0l;
+	
+		delete m_pLogicCacheBase->data;
+		m_pLogicCacheBase->data = 0l;
+	}
+}
+
+
+void Circuit::cacheAndUpdate()
+{
+	LogicCacheNode * node = m_pLogicCacheBase;
+	for ( unsigned i = 0; i < m_logicOutCount; i++ )
+	{
+		if ( m_pLogicOut[i]->outputState() )
+		{
+			if (!node->high)
+				node->high = new LogicCacheNode;
+			
+			node = node->high;
+		}
+		else
+		{
+			if (!node->low)
+				node->low = new LogicCacheNode;
+			
+			node = node->low;
+		}
+	}
+	
+	if ( node->data )
+	{
+		(*m_elementSet->x()) = *node->data;
+		m_elementSet->updateInfo();
+		return;
+	}
+	
+	if ( m_elementSet->containsNonLinear() )
+		m_elementSet->doNonLinear( 150, 1e-10, 1e-13 );
+	else
+		m_elementSet->doLinear(true);
+	
+	node->data = new Vector( m_elementSet->x()->size() );
+	*node->data = *m_elementSet->x();
+}
+
+
+void Circuit::createMatrixMap()
+{
+	m_elementSet->createMatrixMap();
+}
+
+
+bool Circuit::recursivePinAdd( Pin *node, PinList *unassignedNodes, PinList *associated, PinList *nodes )
+{
+	if ( !unassignedNodes->contains(node) )
+		return false;
+	unassignedNodes->remove(node);
+	
+	bool foundGround = node->eqId() == -1;
+	
+	const PinList circuitDependentPins = node->circuitDependentPins();
+	const PinList::const_iterator dEnd = circuitDependentPins.end();
+	for ( PinList::const_iterator it = circuitDependentPins.begin(); it != dEnd; ++it )
+	{
+		if ( !associated->contains(*it) )
+			associated->append(*it);
+	}
+	
+	nodes->append(node);
+	
+	const PinList localConnectedPins = node->localConnectedPins();
+	const PinList::const_iterator end = localConnectedPins.end();
+	for ( PinList::const_iterator it = localConnectedPins.begin(); it != end; ++it )
+		foundGround |= recursivePinAdd( *it, unassignedNodes, associated, nodes );
+	
+	return foundGround;
+}
+
+
+void Circuit::doNonLogic()
+{
+	if ( !m_elementSet || m_cnodeCount+m_branchCount <= 0 )
+		return;
+	
+	if (m_bCanCache)
+	{
+		if ( !m_elementSet->b()->isChanged() && !m_elementSet->matrix()->isChanged() )
+			return;
+		cacheAndUpdate();
+		updateNodalVoltages();
+		m_elementSet->b()->setUnchanged();
+		return;
+	}
+	
+	stepReactive();
+	if ( m_elementSet->containsNonLinear() )
+	{
+		m_elementSet->doNonLinear( 10, 1e-9, 1e-12 );
+		updateNodalVoltages();
+	}
+	else
+	{
+		if ( m_elementSet->doLinear(true) )
+			updateNodalVoltages();
+	}
+}
+
+
+void Circuit::stepReactive()
+{
+	ElementList::iterator listEnd = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != listEnd; ++it )
+	{
+		Element * const e = *it;
+		if ( e && e->isReactive() )
+			(static_cast<Reactive*>(e))->time_step();
+	}
+}
+
+
+void Circuit::updateNodalVoltages()
+{
+	CNode **_cnodes = m_elementSet->cnodes();
+	
+	const PinList::iterator endIt = m_pinList.end();
+	for ( PinList::iterator it = m_pinList.begin(); it != endIt; ++it )
+	{
+		Pin * const node = *it;
+		int i = node->eqId();
+		if ( i == -1 )
+			node->setVoltage(0.);
+		else
+		{
+			const double v = _cnodes[i]->v;
+			node->setVoltage( std::isfinite(v)?v:0. );
+		}
+	}
+}
+
+
+void Circuit::updateCurrents()
+{
+	ElementList::iterator listEnd = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != listEnd; ++it )
+	{
+		(*it)->updateCurrents();
+	}
+}
+
+void Circuit::displayEquations()
+{
+	m_elementSet->displayEquations();
+}
+//END class Circuit
+
+
+
+//BEGIN class LogicCacheNode
+LogicCacheNode::LogicCacheNode()
+{
+	low = 0l;
+	high = 0l;
+	data = 0l;
+}
+
+
+LogicCacheNode::~LogicCacheNode()
+{
+	delete low;
+	delete high;
+	delete data;
+}
+//END class LogicCacheNode
+
+
diff --git a/src/electronics/simulation/circuit.h b/src/electronics/simulation/circuit.h
new file mode 100644
index 0000000..2455edc
--- /dev/null
+++ b/src/electronics/simulation/circuit.h
@@ -0,0 +1,137 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CIRCUIT_H
+#define CIRCUIT_H
+
+#include <qguardedptr.h>
+#include "qstringlist.h"
+#include "qvaluelist.h"
+
+#include "elementset.h"
+
+class CircuitDocument;
+class Wire;
+class Pin;
+class Element;
+class LogicOut;
+
+typedef QValueList<QGuardedPtr<Pin> > PinList;
+typedef QValueList<Element*> ElementList;
+
+
+class LogicCacheNode
+{
+	public:
+		LogicCacheNode();
+		~LogicCacheNode();
+	
+		LogicCacheNode * high;
+		LogicCacheNode * low;
+		Vector * data;
+};
+
+
+/**
+Usage of this class (usually invoked from CircuitDocument):
+(1) Add Wires, Pins and Elements to the class as appropriate
+(2) Call init to initialize the simulation
+(3) Control the simulation with step()
+
+This class can be considered a bridge between the gui-tainted CircuitDocument - specific
+to this implementation, and the pure untainted ElementSet. Please keep it that way.
+
+@short Simulates a collection of components
+@author David Saxton
+*/
+class Circuit
+{
+	public:
+		Circuit();
+		~Circuit();
+	
+		void addPin( Pin *node );
+		void addElement( Element *element );
+	
+		bool contains( Pin *node );
+		bool containsNonLinear() const { return m_elementSet->containsNonLinear(); }
+	
+		void init();
+		/**
+		 * Called after everything else has been setup - before doNonLogic or
+		 * doLogic are called for the first time. Preps the circuit.
+		 */
+		void initCache();
+		/**
+		 * Marks all cached results as invalidated and removes them.
+		 */
+		void setCacheInvalidated();
+		/**
+		 * Solves for non-logic elements
+		 */
+		void doNonLogic();
+		/**
+		 * Solves for logic elements (i.e just does fbSub)
+		 */
+		void doLogic() { m_elementSet->doLinear(false); }
+	
+		void displayEquations();
+		void updateCurrents();
+	
+		void createMatrixMap();
+		/**
+		 * This will identify the ground node and non-ground nodes in the given set.
+		 * Ground will be given the eqId -1, non-ground of 0.
+		 * @param highest The highest ground type of the groundnodes found. If no
+			ground nodes were found, this will be (gt_never-1).
+		 * @returns the number of ground nodes. If all nodes are at or below the
+		 * 			gt_never threshold, then this will be zero.
+		 */
+		static int identifyGround( PinList nodeList, int *highest = 0l );
+		
+		void setNextChanged( Circuit * circuit, unsigned char chain ) { m_pNextChanged[chain] = circuit; }
+		Circuit * nextChanged( unsigned char chain ) const { return m_pNextChanged[chain]; }
+		void setCanAddChanged( bool canAdd ) { m_bCanAddChanged = canAdd; }
+		bool canAddChanged() const { return m_bCanAddChanged; }
+		
+	protected:
+		void cacheAndUpdate();
+		/**
+		 * Update the nodal voltages from those calculated in ElementSet
+		 */
+		void updateNodalVoltages();
+		/**
+		 * Step the reactive elements.
+		 */
+		void stepReactive();
+		/**
+		 * Returns true if any of the nodes are ground
+		 */
+		static bool recursivePinAdd( Pin *node, PinList *unassignedNodes, PinList *associated, PinList *nodes );
+	
+		int m_cnodeCount;
+		int m_branchCount;
+		int m_prepNLCount; // Count until next m_elementSet->prepareNonLinear() is called
+	
+		PinList m_pinList;
+		ElementList m_elementList;
+		ElementSet *m_elementSet;
+		
+		//Stuff for caching
+		bool m_bCanCache;
+		LogicCacheNode * m_pLogicCacheBase;
+		unsigned m_logicOutCount;
+		LogicOut ** m_pLogicOut;
+		
+		bool m_bCanAddChanged;
+		Circuit * m_pNextChanged[2];
+};
+
+#endif
diff --git a/src/electronics/simulation/currentsignal.cpp b/src/electronics/simulation/currentsignal.cpp
new file mode 100644
index 0000000..5e5388d
--- /dev/null
+++ b/src/electronics/simulation/currentsignal.cpp
@@ -0,0 +1,67 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "currentsignal.h"
+#include "element.h"
+#include "matrix.h"
+
+CurrentSignal::CurrentSignal( const double delta, const double current )
+	: Reactive::Reactive(delta)
+{
+	m_current = current;
+	m_oldCurrent = m_newCurrent = 0.;
+	m_numCNodes = 2;
+}
+
+CurrentSignal::~CurrentSignal()
+{
+}
+
+void CurrentSignal::setCurrent( const double i )
+{
+	if ( m_oldCurrent != m_newCurrent ) add_initial_dc();
+	m_newCurrent *= i/m_current; // Instead of calling step again, we can just "adjust" what the current should be
+	m_current = i;
+	add_initial_dc();
+}
+
+
+void CurrentSignal::add_map()
+{
+	// We don't need a map for current signal :-)
+}
+
+
+void CurrentSignal::add_initial_dc()
+{
+	if ( !b_status )
+		return;
+	
+	if ( m_newCurrent == m_oldCurrent )
+		return;
+
+	b_i( 0 ) -= m_newCurrent-m_oldCurrent;
+	b_i( 1 ) += m_newCurrent-m_oldCurrent;
+	
+	m_oldCurrent = m_newCurrent;
+}
+
+void CurrentSignal::updateCurrents()
+{
+	m_cnodeI[1] = m_newCurrent;
+	m_cnodeI[0] = -m_newCurrent;
+}
+
+void CurrentSignal::time_step()
+{
+	add_initial_dc(); // Make sure our old and new are synced
+	m_newCurrent = m_current*advance();
+	add_initial_dc();
+}
diff --git a/src/electronics/simulation/currentsignal.h b/src/electronics/simulation/currentsignal.h
new file mode 100644
index 0000000..b217b79
--- /dev/null
+++ b/src/electronics/simulation/currentsignal.h
@@ -0,0 +1,43 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CURRENTSIGNAL_H
+#define CURRENTSIGNAL_H
+
+#include "reactive.h"
+#include "elementsignal.h"
+
+/**
+@short CurrentSignal
+@author David saxton
+*/
+class CurrentSignal : public Reactive, public ElementSignal
+{
+public:
+	CurrentSignal(  const double delta, const double current );
+	virtual ~CurrentSignal();
+	
+	virtual Element::Type type() const { return Element_CurrentSignal; }
+	void setCurrent( const double current );
+	double current() { return m_current; }
+	virtual void time_step();
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_current; // Current
+	double m_oldCurrent; // Old calculated current
+	double m_newCurrent; // New calculated current
+};
+
+#endif
diff --git a/src/electronics/simulation/currentsource.cpp b/src/electronics/simulation/currentsource.cpp
new file mode 100644
index 0000000..675b0b7
--- /dev/null
+++ b/src/electronics/simulation/currentsource.cpp
@@ -0,0 +1,64 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "currentsource.h"
+#include "elementset.h"
+#include "matrix.h"
+
+CurrentSource::CurrentSource( const double current )
+	: Element::Element()
+{
+	m_i = current;
+	m_numCNodes = 2;
+}
+
+
+CurrentSource::~CurrentSource()
+{
+}
+
+
+void CurrentSource::setCurrent( const double i )
+{
+	if ( i == m_i ) return;
+
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	// Remove the old current
+	m_i = -m_i;
+	add_initial_dc();
+	
+	m_i = i;
+	add_initial_dc();
+}
+
+
+void CurrentSource::add_map()
+{
+	// We don't need a map for current source :-)
+}
+
+
+void CurrentSource::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	b_i( 0 ) -= m_i;
+	b_i( 1 ) += m_i;
+}
+
+
+void CurrentSource::updateCurrents()
+{
+	m_cnodeI[0] = -m_i;
+	m_cnodeI[1] = m_i;
+}
diff --git a/src/electronics/simulation/currentsource.h b/src/electronics/simulation/currentsource.h
new file mode 100644
index 0000000..9b72e0d
--- /dev/null
+++ b/src/electronics/simulation/currentsource.h
@@ -0,0 +1,39 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef CURRENTSOURCE_H
+#define CURRENTSOURCE_H
+
+#include "element.h"
+
+/**
+cnode n0 has current flowing otu of it, cnode n1 has current flowing into it
+@author David Saxton
+@short Current Source
+*/
+class CurrentSource : public Element
+{
+public:
+	CurrentSource( const double current );
+	virtual ~CurrentSource();
+	
+	virtual Type type() const { return Element_CurrentSource; }
+	void setCurrent( const double i );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_i; // Current
+};
+
+#endif
diff --git a/src/electronics/simulation/diode.cpp b/src/electronics/simulation/diode.cpp
new file mode 100644
index 0000000..e13d478
--- /dev/null
+++ b/src/electronics/simulation/diode.cpp
@@ -0,0 +1,198 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+ 
+#include <vector>
+#include "diode.h"
+#include "elementset.h"
+#include "matrix.h"
+
+#include <cmath>
+
+
+//BEGIN class Diode Settings
+DiodeSettings::DiodeSettings()
+{
+	reset();
+}
+
+
+void DiodeSettings::reset()
+{
+	I_S = 1e-15;
+	N = 1.0;
+	V_B = 4.7;
+// 	R = 0.001;
+}
+//END class Diode Settings
+
+
+
+//BEGIN class Diode
+Diode::Diode()
+	: NonLinear()
+{
+	m_numCNodes = 2;
+	g_new = g_old = I_new = I_old = V_prev = 0.0;
+}
+
+
+Diode::~Diode()
+{
+}
+
+
+void Diode::add_map()
+{
+	if (!b_status)
+		return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[1]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+	}
+	
+	if ( !p_cnode[0]->isGround && !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[1]->n(), Map::et_unstable, false );
+		p_A->setUse( p_cnode[1]->n(), p_cnode[0]->n(), Map::et_unstable, false );
+	}
+}
+
+
+void Diode::add_initial_dc()
+{
+	g_new = g_old = I_new = I_old = V_prev = 0.0;
+	update_dc();
+}
+
+
+double Diode::current() const
+{
+	if (!b_status)
+		return 0.0;
+	
+	double I;
+	calcIg( p_cnode[0]->v - p_cnode[1]->v, & I, 0 );
+	
+	return I;
+}
+
+
+void Diode::updateCurrents()
+{
+	if (!b_status)
+		return;
+	
+	m_cnodeI[1] = current();
+	m_cnodeI[0] = -m_cnodeI[1];
+}
+
+
+void Diode::update_dc()
+{
+	if (!b_status)
+		return;
+	
+	calc_eq();
+	
+	A_g( 0, 0 ) += g_new - g_old;
+	A_g( 1, 1 ) += g_new - g_old;
+	A_g( 0, 1 ) -= g_new - g_old;
+	A_g( 1, 0 ) -= g_new - g_old;
+	
+	b_i( 0 ) -= I_new - I_old;
+	b_i( 1 ) += I_new - I_old;
+	
+	g_old = g_new;
+	I_old = I_new;
+}
+
+
+
+#ifndef MIN
+# define MIN(x,y) (((x) < (y)) ? (x) : (y))
+#endif
+
+
+
+void Diode::calc_eq()
+{
+	double I_S = m_diodeSettings.I_S;
+	double N = m_diodeSettings.N;
+	double V_B = m_diodeSettings.V_B;
+// 	double R = m_diodeSettings.R;
+	
+	double v = p_cnode[0]->v - p_cnode[1]->v;
+	
+	// adjust voltage to help convergence
+	double V_crit = diodeCriticalVoltage( I_S, N * V_T );
+	if (V_B != 0 && v < MIN (0, -V_B + 10 * N * V_T))
+	{
+		double V = -(v + V_B);
+		V = diodeVoltage( V, -(V_prev + V_B), V_T * N, V_crit );
+		v = -(V + V_B);
+	}
+	else
+		v = diodeVoltage( v, V_prev, V_T * N, V_crit );
+	
+	V_prev = v;
+	
+	double I_D;
+	calcIg( v, & I_D, & g_new );
+	
+	I_new = I_D - (v * g_new);
+}
+
+
+void Diode::calcIg( double V, double * I_D, double * g ) const
+{
+	double I_S = m_diodeSettings.I_S;
+	double N = m_diodeSettings.N;
+	double V_B = m_diodeSettings.V_B;
+// 	double R = m_diodeSettings.R;
+	
+	double gtiny = (V < - 10 * V_T * N && V_B != 0) ? I_S : 0;
+	
+	if ( V >= (-3 * N * V_T) )
+	{
+		if ( g )
+			*g = diodeConductance( V, I_S, V_T * N ) + gtiny;
+		*I_D = diodeCurrent( V, I_S, V_T * N ) + (gtiny * V);
+	}
+	else if ( V_B == 0 || V >= -V_B )
+	{
+		double a = (3 * N * V_T) / (V * M_E);
+		a = a * a * a;
+		*I_D = (-I_S * (1 + a)) + (gtiny * V);
+		if ( g )
+			*g = ((I_S * 3 * a) / V) + gtiny;
+	}
+	else
+	{
+		double a = exp( -(V_B + V) / N / V_T );
+		*I_D = (-I_S * a) + (gtiny * V);
+		if ( g )
+			*g = I_S * a / V_T / N + gtiny;
+	}
+}
+
+
+void Diode::setDiodeSettings( const DiodeSettings & settings )
+{
+	m_diodeSettings = settings;
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+}
+//END class Diode
+
diff --git a/src/electronics/simulation/diode.h b/src/electronics/simulation/diode.h
new file mode 100644
index 0000000..0b13946
--- /dev/null
+++ b/src/electronics/simulation/diode.h
@@ -0,0 +1,73 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef DIODE_H
+#define DIODE_H
+
+#include "nonlinear.h"
+
+class DiodeSettings
+{
+	public:
+		DiodeSettings();
+		void reset();
+		
+		double I_S;	///< Diode saturation current
+		double N;	///< Emission coefficient
+		double V_B;	///< Reverse breakdown
+// 		double R;	///< Series resistance
+};
+
+
+/**
+This simulates a diode. The simulated diode characteristics are:
+
+@li I_s: Diode saturation current
+@li V_T: Thermal voltage = kT/4 = 25mV at 20 C
+@li n: Emission coefficient, typically between 1 and 2
+@li V_RB: Reverse breakdown (large negative voltage)
+@li G_RB: Reverse breakdown conductance
+@li R_D: Base resistance of diode
+
+@short Simulates the electrical property of diode-ness
+@author David Saxton
+*/
+class Diode : public NonLinear
+{
+	public:
+		Diode();
+		virtual ~Diode();
+	
+		virtual void update_dc();
+		virtual void add_initial_dc();
+		virtual void add_map();
+		virtual Element::Type type() const { return Element_Diode; }
+		DiodeSettings settings() const { return m_diodeSettings; }
+		void setDiodeSettings( const DiodeSettings & settings );
+		/**
+		 * Returns the current flowing through the diode
+		 */
+		double current() const;
+	
+	protected:
+		virtual void updateCurrents();
+		void calc_eq();
+		
+		void calcIg( double V, double * I, double * g ) const;
+		
+		double g_new, g_old;
+		double I_new, I_old;
+		double V_prev;
+	
+		DiodeSettings m_diodeSettings;
+};
+
+#endif
+
diff --git a/src/electronics/simulation/element.cpp b/src/electronics/simulation/element.cpp
new file mode 100644
index 0000000..2411897
--- /dev/null
+++ b/src/electronics/simulation/element.cpp
@@ -0,0 +1,193 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "element.h"
+#include "elementset.h"
+
+#include <assert.h>
+
+#include <kdebug.h>
+
+Element::Element()
+{
+	b_status = false;
+	p_A = 0l;
+	p_eSet = 0l;
+	p_b = 0l;
+	b_componentDeleted = false;
+	b_eSetDeleted = true;
+	
+	for ( int i=0; i<8; i++ )
+		p_cnode[i] = 0l;
+	
+	resetCurrents();
+	
+	for ( int i=0; i<4; i++ )
+		p_cbranch[i] = 0l;
+	
+	m_numCBranches = 0;
+	m_numCNodes = 0;
+}
+
+Element::~ Element()
+{
+}
+
+void Element::resetCurrents()
+{
+	for ( int i=0; i<8; i++ )
+		m_cnodeI[i] = 0.0;
+}
+
+void Element::setElementSet( ElementSet *c )
+{
+	assert(!b_componentDeleted);
+	assert(!p_eSet);
+	if (!c) return elementSetDeleted();
+	b_eSetDeleted = false;
+	p_eSet = c;
+	p_A = p_eSet->matrix();
+	p_b = p_eSet->b();
+	updateStatus();
+}
+
+void Element::componentDeleted()
+{
+// 	assert(!b_componentDeleted);
+	if (b_componentDeleted)
+	{
+		// Something strange happened here....
+	}
+	if (b_eSetDeleted) return delete this;
+	b_componentDeleted = true;
+	b_status = false;
+// 	kdDebug() << "Element::componentDeleted(): Setting b_status to false, this="<<this<<endl;
+	
+	p_eSet = 0l;
+	p_A = 0l;
+	p_b = 0l;
+	setCNodes();
+	setCBranches();
+}
+
+void Element::elementSetDeleted()
+{
+// 	assert(!b_eSetDeleted);
+	if (b_eSetDeleted)
+	{
+		// Something strange happened here....
+	}
+	if (b_componentDeleted) return delete this;
+	b_eSetDeleted = true;
+	b_status = false;
+// 	kdDebug() << "Element::elementSetDeleted(): Setting b_status to false, this="<<this<<endl;
+	
+	p_eSet = 0l;
+	p_A = 0l;
+	p_b = 0l;
+	setCNodes();
+	setCBranches();
+}
+
+
+void Element::setCNodes( const int n0, const int n1, const int n2, const int n3 )
+{
+	if ( !p_eSet )
+	{
+// 		cerr << "Element::setCNodes: can't set nodes without circuit!"<<endl;
+		for ( int i=0; i<8; i++ )
+			p_cnode[i] = 0l;
+		return;
+	}
+	
+	// MAX_CNODES-1 should match the last array index below.
+	assert( MAX_CNODES == 4 );
+	p_cnode[0] = (n0>-1)?p_eSet->cnodes()[n0]:(n0==-1?p_eSet->ground():0l);
+	p_cnode[1] = (n1>-1)?p_eSet->cnodes()[n1]:(n1==-1?p_eSet->ground():0l);
+	p_cnode[2] = (n2>-1)?p_eSet->cnodes()[n2]:(n2==-1?p_eSet->ground():0l);
+	p_cnode[3] = (n3>-1)?p_eSet->cnodes()[n3]:(n3==-1?p_eSet->ground():0l);
+	updateStatus();
+}
+
+void Element::setCBranches( const int b0, const int b1, const int b2, const int b3 )
+{
+	if ( !p_eSet )
+	{
+// 		cerr << "Element::setCBranches: can't set branches without circuit!"<<endl;
+		for ( int i=0; i<4; i++ ) p_cbranch[i] = 0l;
+		return;
+	}
+	p_cbranch[0] = (b0>-1)?p_eSet->cbranches()[b0]:0l;
+	p_cbranch[1] = (b1>-1)?p_eSet->cbranches()[b1]:0l;
+	p_cbranch[2] = (b2>-1)?p_eSet->cbranches()[b2]:0l;
+	p_cbranch[3] = (b3>-1)?p_eSet->cbranches()[b3]:0l;
+	updateStatus();
+}
+
+bool Element::updateStatus()
+{
+	// First, set status to false if all nodes in use are ground
+	b_status = false;
+	for ( int i=0; i<m_numCNodes; i++ )
+	{
+		b_status |= p_cnode[i]?!p_cnode[i]->isGround:false;
+	}
+	
+	// Set status to false if any of the nodes are not set
+	for ( int i=0; i<m_numCNodes; i++ )
+	{
+		if (!p_cnode[i]) b_status = false;
+	}
+	
+	// Finally, set status to false if not all the required branches are set
+	for ( int i=0; i<m_numCBranches; i++ )
+	{
+		if (!p_cbranch[i]) b_status = false;
+	}
+	
+	// Finally, check for various pointers
+	if ( !p_eSet || !p_A || !p_b ) b_status = false;
+	
+	if (!b_status)
+	{
+		resetCurrents();
+	}
+	// And return the status :-)
+// 	kdDebug() << "Element::updateStatus(): Setting b_status to "<<(b_status?"true":"false")<<" this="<<this<<endl;
+	return b_status;
+}
+
+double Element::cbranchCurrent( const int branch )
+{
+	if ( !b_status || branch<0 || branch>=m_numCBranches ) return 0.;
+	return (*p_cbranch)[branch].i;
+}
+
+double Element::cnodeVoltage( const int node )
+{
+	if ( !b_status || node<0 || node>=m_numCNodes ) return 0.;
+	return (*p_cnode)[node].v;
+}
+
+
+CNode::CNode()
+{
+	m_n = 0;
+	v = 0.0;
+	isGround = false;
+}
+
+CBranch::CBranch()
+{
+	m_n = 0;
+	i = 0.0;
+}
+
+
diff --git a/src/electronics/simulation/element.h b/src/electronics/simulation/element.h
new file mode 100644
index 0000000..e05de46
--- /dev/null
+++ b/src/electronics/simulation/element.h
@@ -0,0 +1,255 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef ELEMENT_H
+#define ELEMENT_H
+
+#include "elementset.h"
+#include "matrix.h"
+
+class ElementSet;
+class Vector;
+typedef unsigned int uint;
+
+const double T = 300.; // Temperature in Kelvin
+const double K = 1.38e-23; // Boltzmann's constant
+const double q = 1.602e-19; // Charge on an electron
+const double V_T = K*T/q; // Thermal voltage
+const double gmin = 1e-12; // Minimum parallel conductance used in dc domain
+
+class CNode
+{
+public:
+	CNode();
+	void set_n( const uint n ) { m_n=n; }
+	uint n() { return m_n; }
+	double v; // Voltage on node. This is set from the last calculated voltage.
+	bool isGround; // True for ground nodes. Obviously, you should ignore n and v if this is true
+private:
+	uint m_n; // CNode number
+};
+
+class CBranch
+{
+public:
+	CBranch();
+	void set_n( const uint n ) { m_n=n; }
+	uint n() { return m_n; }
+	double i; // Current flowing through branch. This is set from the last calculated current.
+private:
+	uint m_n; // CBranch number
+};
+
+const int MAX_CNODES = 4;
+
+// Default node number that represents no node (remember that
+// Ground node is -1, and the rest are numbered from 0 to n-1
+const int noCNode = -2;
+// Likewise for branch (although there is no "ground" branch;
+// it is merely -2 for likeness with noCNode)
+const int noBranch = -2;
+
+/**
+@short Represents a circuit element (such as resistance)
+@author David Saxton
+*/
+class Element
+{
+public:
+	enum Type
+	{
+		Element_BJT,
+		Element_Capacitance,
+		Element_CCCS,
+		Element_CCVS,
+		Element_CurrentSignal,
+		Element_CurrentSource,
+		Element_Diode,
+		Element_Inductance,
+		Element_LogicIn,
+		Element_LogicOut,
+		Element_OpAmp,
+		Element_Resistance,
+		Element_VCCS,
+		Element_VCVS,
+		Element_VoltagePoint,
+		Element_VoltageSignal,
+		Element_VoltageSource
+	};
+	
+	Element();
+	virtual ~Element();
+	/**
+	 * This must be called when the circuit is changed. The function will get
+	 * all the required pointers from ElementSet
+	 */
+	virtual void setElementSet( ElementSet *c );
+	/**
+	 * Returns a pointer to the current element set
+	 */
+	ElementSet *elementSet() { return p_eSet; }
+	/**
+	 * Tells the element which nodes to use. Remember that -1 is ground. You
+	 * should refer to the individual elements for which nodes are used for what.
+	 */
+	void setCNodes( const int n0 = noCNode, const int n1 = noCNode, const int n2 = noCNode, const int n3 = noCNode );
+	/**
+	 * Tells the element it's branch numbers (if it should have one). Not
+	 * all elements use this.
+	 */
+	void setCBranches( const int b0 = noBranch, const int b1 = noBranch, const int b2 = noBranch, const int b3 = noBranch );
+	/**
+	 * Returns a pointer to the given CNode
+	 */
+	CNode *cnode( const uint num ) { return p_cnode[num]; }
+	/**
+	 * Returns a pointer to the given CNode
+	 */
+	CBranch *cbranch( const uint num ) { return p_cbranch[num]; }
+	/**
+	 * Returns the number of branches used by the element
+	 */
+	int numCBranches() { return m_numCBranches; }
+	/**
+	 * Returns the number of circuit nodes used by the element
+	 */
+	int numCNodes() { return m_numCNodes; }
+	/**
+	 * Call this function to tell the element to calculate the
+	 * current flowing *into* it's cnodes *from* the element. You
+	 * can get the currents with m_cnodeI. Child class must implement this function.
+	 */
+	virtual void updateCurrents() = 0;
+	/**
+	 * Returns true for reactive elements that need stepping for numerical-integration
+	 * (such as capacitors)
+	 */
+	virtual bool isReactive() { return false; }
+	/**
+	 * Returns true for NonLinear elements that need iteration to converge to a solution
+	 * as the matrix A is a function of x.
+	 */
+	virtual bool isNonLinear() { return false; }
+	/**
+	 * Returns the type of element
+	 */
+	virtual Type type() const = 0;
+	/**
+	 * Call this function to tell the element to add its map to the matrix in use
+	 */
+	virtual void add_map() {};
+	/**
+	 * Does the required MNA stuff. This should be called from ElementSet when necessary.
+	 */
+	virtual void add_initial_dc() = 0;
+	/**
+	 * This is called from the Component destructor. When elementSetDeleted has
+	 * also been called, this class will delete itself.
+	 */
+	 void componentDeleted();
+	 void elementSetDeleted();
+	
+	double m_cnodeI[8]; ///< Current flowing into the cnodes from the element
+	double cbranchCurrent( const int branch );
+	double cnodeVoltage( const int node );
+	
+protected:
+	/**
+	 * Resets all calculated currents in the nodes to 0
+	 */
+	void resetCurrents();
+	
+	inline double & A_g( uint i, uint j );
+	inline double & A_b( uint i, uint j );
+	inline double & A_c( uint i, uint j );
+	inline double & A_d( uint i, uint j );
+	
+	inline double & b_i( uint i );
+	inline double & b_v( uint i );
+	
+	ElementSet *p_eSet;
+	Matrix *p_A;
+	Vector *p_b;
+	CNode *p_cnode[MAX_CNODES];
+	CBranch *p_cbranch[4];
+	
+	/**
+	 * True when the element can do add_initial_dc(), i.e. when it has
+	 * pointers to the circuit, and at least one of its nodes is not ground.
+	 */
+	bool b_status;
+	/**
+	 * Update the status, returning b_status
+	 */
+	virtual bool updateStatus();
+	/**
+	 * Set by child class - the number of branches that the element uses
+	 * Typically, this is 0, but could be 1 (e.g. independent voltage source)
+	 * or 2 (e.g. cccs)
+	 */
+	int m_numCBranches;
+	/**
+	 * Set by child class - the number of circuit nodes that the element uses
+	 */
+	int m_numCNodes;
+	
+private:
+	bool b_componentDeleted;
+	bool b_eSetDeleted; 
+	double m_temp;
+};
+
+
+double & Element::A_g( uint i, uint j )
+{
+	if ( p_cnode[i]->isGround || p_cnode[j]->isGround )
+		return m_temp;
+	return p_A->g( p_cnode[i]->n(), p_cnode[j]->n() );
+}
+
+
+double & Element::A_b( uint i, uint j )
+{
+	if ( p_cnode[i]->isGround )
+		return m_temp;
+	return p_A->b( p_cnode[i]->n(), p_cbranch[j]->n() );
+}
+
+
+double & Element::A_c( uint i, uint j )
+{
+	if ( p_cnode[j]->isGround )
+		return m_temp;
+	return p_A->c( p_cbranch[i]->n(), p_cnode[j]->n() );
+}
+
+
+double & Element::A_d( uint i, uint j )
+{
+	return p_A->d( p_cbranch[i]->n(), p_cbranch[j]->n() );
+}
+
+
+
+double & Element::b_i( uint i )
+{
+	if ( p_cnode[i]->isGround )
+		return m_temp;
+	
+	return (*p_b)[ p_cnode[i]->n() ];
+}
+
+
+double & Element::b_v( uint i )
+{
+	return (*p_b)[ p_eSet->cnodeCount() + p_cbranch[i]->n() ];
+}
+
+#endif
diff --git a/src/electronics/simulation/elementset.cpp b/src/electronics/simulation/elementset.cpp
new file mode 100644
index 0000000..25057c2
--- /dev/null
+++ b/src/electronics/simulation/elementset.cpp
@@ -0,0 +1,253 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "bjt.h"
+#include "circuit.h"
+#include "elementset.h"
+#include "element.h"
+#include "logic.h"
+#include "matrix.h"
+#include "nonlinear.h"
+#include "vec.h"
+
+#include <kdebug.h>
+
+#include <cmath>
+#include <iostream>
+#include <assert.h>
+
+ElementSet::ElementSet( Circuit * circuit, const int n, const int m )
+{
+	m_pCircuit = circuit;
+	m_cn = n;
+	m_cb = m;
+	p_logicIn = 0l;
+	p_A = new Matrix( m_cn, m_cb );
+	p_b = new Vector(m_cn+m_cb);
+	p_x = new Vector(m_cn+m_cb);
+	p_x_prev = new Vector(m_cn+m_cb);
+	m_cbranches = new CBranch*[m_cb];
+	m_cnodes = new CNode*[m_cn];
+	for ( uint i=0; i<m_cn; i++ )
+	{
+		m_cnodes[i] = new CNode();
+		m_cnodes[i]->set_n(i);
+	}
+	for ( uint i=0; i<m_cb; i++ )
+	{
+		m_cbranches[i] = new CBranch();
+		m_cbranches[i]->set_n(i);
+	}
+	m_ground = new CNode();
+	m_ground->isGround = true;
+	b_containsNonLinear = false;
+}
+
+ElementSet::~ElementSet()
+{
+	const ElementList::iterator end = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != end; ++it )
+	{
+		// Note: By calling setElementSet(0l), we might have deleted it (the Element will commit
+		// suicide when both the the ElementSet and Component to which it belongs have deleted
+		// themselves). So be very careful it you plan to do anything with the (*it) pointer
+		if (*it) (*it)->elementSetDeleted();
+	}
+	for ( uint i=0; i<m_cn; i++ )
+	{
+		delete m_cnodes[i];
+	}
+	for ( uint i=0; i<m_cb; i++ )
+	{
+		delete m_cbranches[i];
+	}
+	delete[] m_cbranches;
+	delete[] m_cnodes;
+	delete[] p_logicIn;
+	delete m_ground;
+	delete p_A;
+	delete p_b;
+	delete p_x;
+	delete p_x_prev;
+}
+
+
+void ElementSet::setCacheInvalidated()
+{
+	m_pCircuit->setCacheInvalidated();
+}
+
+
+void ElementSet::addElement( Element *e )
+{
+	if ( !e || m_elementList.contains(e) ) return;
+	e->setElementSet(this);
+	m_elementList.append(e);
+	if ( e->isNonLinear() )
+	{
+		b_containsNonLinear = true;
+		m_cnonLinearList.append( static_cast<NonLinear*>(e) );
+	}
+}
+
+
+void ElementSet::createMatrixMap()
+{
+	p_A->createMap();
+	
+	
+	// And do our logic as well...
+	
+	m_clogic = 0;
+	ElementList::iterator end = m_elementList.end();
+	for ( ElementList::iterator it = m_elementList.begin(); it != end; ++it )
+	{
+		if ( dynamic_cast<LogicIn*>(*it) )
+			m_clogic++;
+	}
+	
+	p_logicIn = new LogicIn*[m_clogic];
+	int i=0;
+	for ( ElementList::iterator it = m_elementList.begin(); it != end; ++it )
+	{
+		if ( LogicIn * in = dynamic_cast<LogicIn*>(*it) )
+			p_logicIn[i++] = in;
+	}
+}
+
+
+void ElementSet::doNonLinear( int maxIterations, double maxErrorV, double maxErrorI )
+{
+// 	p_x_prev->reset();
+	
+	// And now tell the cnodes and cbranches about their new voltages & currents
+	updateInfo();
+	
+	const NonLinearList::iterator end = m_cnonLinearList.end();
+	
+	int k = 0;
+	do
+	{
+		// Tell the nonlinear elements to update its J, A and b from the newly calculated x
+		for ( NonLinearList::iterator it = m_cnonLinearList.begin(); it != end; ++it )
+			(*it)->update_dc();
+		
+		*p_x = *p_b;
+		p_A->performLU();
+		p_A->fbSub(p_x);
+		updateInfo();
+		
+		// Now, check for convergence
+		bool converged = true;
+		for ( unsigned i = 0; i < m_cn; ++i )
+		{
+			double diff = std::abs( (*p_x_prev)[i] - (*p_x)[i] );
+			if ( diff > maxErrorI )
+			{
+				converged = false;
+				break;
+			}
+		}
+		if ( converged )
+		{
+			for ( unsigned i = m_cn; i < m_cn+m_cb; ++i )
+			{
+				double diff = std::abs( (*p_x_prev)[i] - (*p_x)[i] );
+				if ( diff > maxErrorV )
+				{
+					converged = false;
+					break;
+				}
+			}
+		}
+		
+		*p_x_prev = *p_x;
+		
+		if ( converged )
+			break;
+	}
+	while ( ++k < maxIterations );
+}
+
+
+bool ElementSet::doLinear( bool performLU )
+{
+	if ( b_containsNonLinear || (!p_b->isChanged() && ((performLU && !p_A->isChanged()) || !performLU)) )
+		return false;
+	
+	if (performLU)
+		p_A->performLU();
+	
+	*p_x = *p_b;
+	p_A->fbSub(p_x);
+	updateInfo();
+	p_b->setUnchanged();
+	
+	return true;
+}
+
+
+void ElementSet::updateInfo()
+{
+	for ( uint i=0; i<m_cn; i++ )
+	{
+		const double v = (*p_x)[i];
+		if (std::isfinite(v)) {
+			m_cnodes[i]->v = v;
+		} else {
+			(*p_x)[i] = 0.;
+			m_cnodes[i]->v = 0.;
+		}
+	}
+	for ( uint i=0; i<m_cb; i++ )
+	{
+		// NOTE: I've used lowercase and uppercase "I" here, so be careful!
+		const double I = (*p_x)[i+m_cn];
+		if (std::isfinite(I)) {
+			m_cbranches[i]->i = I;
+		} else {
+			(*p_x)[i+m_cn] = 0.;
+			m_cbranches[i]->i = 0.;
+		}
+	}
+	
+	// Tell logic to check themselves
+	for ( uint i=0; i<m_clogic; ++i )
+	{
+		p_logicIn[i]->check();
+	}
+}
+
+
+void ElementSet::displayEquations()
+{
+	std::cout.setf(std::ios_base::fixed);
+	std::cout.precision(5);
+	std::cout.setf(std::ios_base::showpoint);
+	std::cout << "A x = b :"<<std::endl;
+	for ( uint i=0; i<m_cn+m_cb; i++ )
+	{
+		std::cout << "( ";
+		for ( uint j=0; j<m_cn+m_cb; j++ )
+		{
+			const double value = p_A->g(i,j);
+// 			if		( value > 0 ) cout <<"+";
+// 			else if ( value == 0 ) cout <<" ";
+			std::cout.width(10);
+			std::cout << value<<" ";
+		}
+		std::cout << ") ( "<<(*p_x)[i]<<" ) = ( ";
+		std::cout<<(*p_b)[i]<<" )"<<std::endl;
+	}
+	std::cout << "A_LU:"<<std::endl;
+	p_A->displayLU();
+}
+
+
diff --git a/src/electronics/simulation/elementset.h b/src/electronics/simulation/elementset.h
new file mode 100644
index 0000000..c7ef7ca
--- /dev/null
+++ b/src/electronics/simulation/elementset.h
@@ -0,0 +1,131 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef ELEMENTSET_H
+#define ELEMENTSET_H
+
+#include <vector>
+
+#include <qvaluelist.h>
+
+class CBranch;
+class Circuit;
+class CNode;
+class Element;
+class ElementSet;
+class LogicIn;
+class Matrix;
+class NonLinear;
+class Vector;
+
+typedef QValueList<Element*> ElementList;
+typedef QValueList<NonLinear*> NonLinearList;
+
+/**
+Steps in simulation of a set of elements:
+(1) Create this class with given number of nodes "n" and voltage sources "m"
+(2) Add the various elements with addElement.
+(3) Call performDC()
+(4) Get the nodal voltages and voltage currents with x()
+(5) Repeat steps 3 and 4 if necessary for further transient analysis.
+
+This class shouldn't be confused with the Circuit class, but considered a helper class to Circuit.
+Circuit will handle the simulation of a set of components over time. This just finds the DC-operating
+point of the circuit for a given set of elements.
+
+@short Handles a set of circuit elements
+@author David Saxton
+*/
+class ElementSet
+{
+public:
+	/**
+	 * Create a new circuit, with "n" nodes and "m" voltage sources.
+	 * After creating the circuit, you must call setGround to specify
+	 * the ground nodes, before adding any elements.
+	 */
+	ElementSet( Circuit * circuit, const int n, const int m );
+	/**
+	 * Destructor. Note that only the matrix and supporting data is deleted.
+	 * i.e. Any elements added to the circuit will not be deleted.
+	 */
+	~ElementSet();
+	Circuit * circuit() const { return m_pCircuit; }
+	void addElement( Element *e );
+	void setCacheInvalidated();
+	/**
+	 * Returns the matrix in use. This is created once on the creation of the ElementSet
+	 * class, and deleted in the destructor, so the pointer returned will never change.
+	 */
+	Matrix *matrix() const { return p_A; }
+	/**
+	 * Returns the vector for b (i.e. the independent currents & voltages)
+	 */
+	Vector *b() const { return p_b; }
+	/**
+	 * Returns the vector for x (i.e. the currents & voltages at the branches and nodes)
+	 */
+	Vector *x() const { return p_x; }
+	/**
+	 * @return if we have any nonlinear elements (e.g. diodes, tranaistors).
+	 */
+	bool containsNonLinear() const { return b_containsNonLinear; }
+	/**
+	 * Solves for nonlinear elements, or just does linear if it doesn't contain
+	 * any nonlinear.
+	 */
+	void doNonLinear( int maxIterations, double maxErrorV = 1e-9, double maxErrorI = 1e-12 );
+	/**
+	 * Solves for linear and logic elements.
+	 * @returns true if anything changed
+	 */
+	bool doLinear( bool performLU );
+	CBranch **cbranches() const { return m_cbranches; }
+	CNode **cnodes() const { return m_cnodes; }
+	CNode *ground() const { return m_ground; }
+	/**
+	 * Returns the number of nodes in the circuit (excluding ground 'nodes')
+	 */
+	int cnodeCount() const { return m_cn; }
+	/**
+	 * Returns the number of voltage sources in the circuit
+	 */
+	int cbranchCount() const { return m_cb; }
+	
+	void createMatrixMap();
+	/**
+	 * Displays the matrix equations Ax=b and J(dx)=-r
+	 */
+	void displayEquations();
+	/**
+	 * Update the nodal voltages and branch currents from the x vector
+	 */
+	void updateInfo();
+	
+private:
+	Matrix *p_A;
+	Vector *p_x;
+	Vector *p_x_prev;
+	Vector *p_b;
+	uint m_cn;
+	uint m_cb;
+	ElementList m_elementList;
+	NonLinearList m_cnonLinearList;
+	CBranch **m_cbranches; // Pointer to an array of cbranches
+	CNode **m_cnodes; // Pointer to an array of cnodes
+	CNode *m_ground;
+	uint m_clogic;
+	LogicIn **p_logicIn;
+	bool b_containsNonLinear;
+	Circuit * m_pCircuit;
+};
+
+#endif
+
diff --git a/src/electronics/simulation/elementsignal.cpp b/src/electronics/simulation/elementsignal.cpp
new file mode 100644
index 0000000..31d7d78
--- /dev/null
+++ b/src/electronics/simulation/elementsignal.cpp
@@ -0,0 +1,64 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementsignal.h"
+#include <cmath>
+
+ElementSignal::ElementSignal()
+{
+	m_type = ElementSignal::st_sinusoidal;
+	m_time = 0.;
+	m_frequency = 0.;
+}
+
+ElementSignal::~ElementSignal()
+{
+}
+
+void ElementSignal::setStep( double delta, Type type, double frequency )
+{
+	m_type = type;
+	m_delta = delta;
+	m_frequency = frequency;
+	m_omega = 6.283185307179586*m_frequency;
+	m_time = 1./(4.*m_frequency);
+}
+
+double ElementSignal::advance()
+{
+	m_time += m_delta;
+	if ( m_time >= 1./m_frequency ) m_time -= 1./m_frequency;
+	
+	switch (m_type)
+	{
+		case ElementSignal::st_sawtooth:
+		{
+			// TODO Sawtooth signal
+			return 0.;
+		}
+		case ElementSignal::st_square:
+		{
+			return (sin(m_time*m_omega)>=0)?1:-1;
+		}
+		case ElementSignal::st_triangular:
+		{
+			// TODO Triangular signal
+			return 0.;
+		}
+		case ElementSignal::st_sinusoidal:
+		default:
+		{
+			return sin(m_time*m_omega);
+		}
+	}
+}
+
+
+
diff --git a/src/electronics/simulation/elementsignal.h b/src/electronics/simulation/elementsignal.h
new file mode 100644
index 0000000..c26bea1
--- /dev/null
+++ b/src/electronics/simulation/elementsignal.h
@@ -0,0 +1,45 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef ELEMENTSIGNAL_H
+#define ELEMENTSIGNAL_H
+
+/**
+@short Provides different signals
+@author David Saxton
+*/
+class ElementSignal
+{
+public:
+	enum Type
+	{
+		st_sinusoidal,
+		st_square,
+		st_sawtooth,
+		st_triangular
+	};
+	ElementSignal();
+	~ElementSignal();
+	
+	void setStep( double delta, Type type, double frequency );
+	/**
+	 * Advances the timer, returns amplitude (between -1 and 1)
+	 */
+	double advance();
+
+protected:
+	Type m_type;
+	double m_time;
+	double m_frequency;
+	double m_delta;
+	double m_omega; // Used for sinusoidal signal
+};
+
+#endif
diff --git a/src/electronics/simulation/inductance.cpp b/src/electronics/simulation/inductance.cpp
new file mode 100644
index 0000000..22c5f9e
--- /dev/null
+++ b/src/electronics/simulation/inductance.cpp
@@ -0,0 +1,126 @@
+/***************************************************************************
+ *   Copyright (C) 2005 by David Saxton                                    *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "inductance.h"
+#include "matrix.h"
+
+Inductance::Inductance( double inductance, double delta )
+	: Reactive(delta)
+{
+	m_inductance = inductance;
+	r_eq_old = v_eq_old = 0.0;
+	m_numCNodes = 2;
+	m_numCBranches = 1;
+	setMethod( Inductance::m_euler );
+}
+
+
+Inductance::~Inductance()
+{
+}
+
+
+void Inductance::setInductance( double i )
+{
+	m_inductance = i;
+}
+
+
+void Inductance::add_initial_dc()
+{
+	A_c( 0, 0 ) = 1;
+	A_b( 0, 0 ) = 1;
+	A_c( 0, 1 ) = -1;
+	A_b( 1, 0 ) = -1;
+	
+	// The adding of r_eg and v_eq will be done for us by time_step.
+	// So for now, just reset the constants used.
+	r_eq_old = v_eq_old = 0.0;
+}
+
+
+void Inductance::updateCurrents()
+{
+	if (!b_status)
+		return;
+	m_cnodeI[0] = p_cbranch[0]->i;
+	m_cnodeI[1] = -m_cnodeI[0];
+}
+
+
+void Inductance::add_map()
+{
+	if (!b_status)
+		return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+	}
+	
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_constant, true );
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+	}
+	
+	p_A->setUse_d( p_cbranch[0]->n(), p_cbranch[0]->n(), Map::et_unstable, false );
+}
+
+
+void Inductance::time_step()
+{
+	if (!b_status) return;
+	
+	double i = p_cbranch[0]->i;
+	double v_eq_new = 0.0, r_eq_new = 0.0;
+	
+	if ( m_method == Inductance::m_euler )
+	{
+		r_eq_new = m_inductance/m_delta;
+		v_eq_new = -i*r_eq_new;
+	}
+	else if ( m_method == Inductance::m_trap ) {
+		// TODO Implement + test trapezoidal method
+		r_eq_new = 2.0*m_inductance/m_delta;
+	}
+	
+	if ( r_eq_old != r_eq_new )
+	{
+		A_d( 0, 0 ) -= r_eq_new - r_eq_old;
+	}
+	
+	if ( v_eq_new != v_eq_old )
+	{
+		b_v( 0 ) += v_eq_new - v_eq_old;
+	}
+	
+	r_eq_old = r_eq_new;
+	v_eq_old = v_eq_new;
+}
+
+
+bool Inductance::updateStatus()
+{
+	b_status = Reactive::updateStatus();
+	if ( m_method == Inductance::m_none )
+		b_status = false;
+	return b_status;
+}
+
+
+void Inductance::setMethod( Method m )
+{
+	m_method = m;
+	updateStatus();
+}
+
diff --git a/src/electronics/simulation/inductance.h b/src/electronics/simulation/inductance.h
new file mode 100644
index 0000000..46ccb09
--- /dev/null
+++ b/src/electronics/simulation/inductance.h
@@ -0,0 +1,55 @@
+/***************************************************************************
+ *   Copyright (C) 2005 by David Saxton                                    *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef INDUCTANCE_H
+#define INDUCTANCE_H
+
+#include "reactive.h"
+
+/**
+
+@author David Saxton
+*/
+class Inductance : public Reactive
+{
+	public:
+		enum Method
+		{
+			m_none, // None
+			m_euler, // Backward Euler
+			m_trap // Trapezoidal (currently unimplemented)
+		};
+		Inductance( double capacitance, double delta );
+		virtual ~Inductance();
+	
+		virtual Type type() const { return Element_Inductance; }
+		/**
+		 * Set the stepping use for numerical integration of inductance, and the
+		 * interval between successive updates.
+		 */
+		void setMethod( Method m );
+		virtual void time_step();
+		virtual void add_initial_dc();
+		void setInductance( double i );
+		virtual void add_map();
+	
+	protected:
+		virtual void updateCurrents();
+		virtual bool updateStatus();
+	
+	private:
+		double m_inductance; // Inductance
+		Method m_method; // Method of integration
+		
+		double r_eq_old;
+		double v_eq_old;
+};
+
+#endif
diff --git a/src/electronics/simulation/logic.cpp b/src/electronics/simulation/logic.cpp
new file mode 100644
index 0000000..031dd2e
--- /dev/null
+++ b/src/electronics/simulation/logic.cpp
@@ -0,0 +1,319 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include <vector>
+#include "circuit.h"
+#include "elementset.h"
+#include "logic.h"
+#include "matrix.h"
+#include "simulator.h"
+#include "src/core/ktlconfig.h"
+
+LogicIn::LogicIn( LogicConfig config )
+	: Element::Element()
+{
+	m_config = config;
+	m_pCallbackFunction = 0l;
+	m_numCNodes = 1;
+	m_bLastState = false;
+	m_pNextLogic = 0l;
+	setLogic(getConfig());
+}
+
+LogicIn::~LogicIn()
+{
+	Simulator::self()->removeLogicInReferences(this);
+}
+
+
+void LogicIn::setCallback( CallbackClass * object, CallbackPtr func )
+{
+	m_pCallbackFunction = func;
+	m_pCallbackObject = object;
+}
+
+
+void LogicIn::check()
+{
+	if (!b_status)
+		return;
+	
+	bool newState;
+	if (m_bLastState)
+	{
+		// Was high, will still be high unless voltage is less than falling trigger
+		newState = p_cnode[0]->v > m_config.fallingTrigger;
+	}
+	else
+	{
+		// Was low, will still be low unless voltage is more than rising trigger
+		newState = p_cnode[0]->v > m_config.risingTrigger;
+	}
+	
+	if ( m_pCallbackFunction && (newState != m_bLastState) )
+	{
+		m_bLastState = newState;
+		(m_pCallbackObject->*m_pCallbackFunction)(newState);
+	}
+	m_bLastState = newState;
+}
+
+
+void LogicIn::setLogic( LogicConfig config )
+{
+	m_config = config;
+	check();
+}
+
+
+void LogicIn::setElementSet( ElementSet *c )
+{
+	if (c)
+		m_pNextLogic = 0l;
+	else
+		m_cnodeI[0] = 0.;
+	
+	Element::setElementSet(c);
+}
+
+
+void LogicIn::add_initial_dc()
+{
+}
+
+
+void LogicIn::updateCurrents()
+{
+}
+
+
+LogicConfig LogicIn::getConfig()
+{
+	LogicConfig c;
+	c.risingTrigger = KTLConfig::logicRisingTrigger();
+	c.fallingTrigger = KTLConfig::logicFallingTrigger();
+	c.output = KTLConfig::logicOutputHigh();
+	c.highImpedance = KTLConfig::logicOutputHighImpedance();
+	c.lowImpedance = KTLConfig::logicOutputLowImpedance();
+	return c;
+}
+
+
+LogicOut::LogicOut( LogicConfig config, bool _high )
+	: LogicIn(config)
+{
+	m_bCanAddChanged = true;
+	m_bOutputHighConductanceConst = false;
+	m_bOutputLowConductanceConst = false;
+	m_bOutputHighVoltageConst = false;
+	m_pNextChanged[0] = m_pNextChanged[1] = 0l;
+	m_pSimulator = 0l;
+	m_bUseLogicChain = false;
+	b_state = false;
+	m_numCNodes = 1;
+	m_vHigh = m_gHigh = m_gLow = 0.0;
+	m_old_g_out = m_g_out = 0.0;
+	m_old_v_out = m_v_out = 0.0;
+	setHigh(_high);
+	
+	// Although we already call this function in LogicIn's constructor, our
+	// virtual function will not have got called, so we have to call it again.
+	setLogic(getConfig());
+}
+
+LogicOut::~LogicOut()
+{
+	if (!m_pSimulator)
+		m_pSimulator = Simulator::self();
+	
+	// Note that although this function will get called in the destructor of
+	// LogicIn, we must call it here as well as it needs to be called before
+	// removeLogicOutReferences(this) is called.
+	m_pSimulator->removeLogicInReferences(this);
+	
+	m_pSimulator->removeLogicOutReferences(this);
+}
+
+
+void LogicOut::setUseLogicChain( bool use )
+{
+	if (!m_pSimulator)
+		m_pSimulator = Simulator::self();
+	
+	m_bUseLogicChain = use;
+	if (use)
+		setElementSet(0l);
+}
+
+
+void LogicOut::setElementSet( ElementSet *c )
+{
+	if (!m_pSimulator)
+		m_pSimulator = Simulator::self();
+	
+	if (c)
+	{
+		m_bUseLogicChain = false;
+		m_pNextChanged[0] = m_pNextChanged[1] = 0l;
+	}
+	
+	// NOTE Make sure that the next two lines are the same as those in setHigh and setLogic
+	m_g_out = b_state ? m_gHigh : m_gLow;
+	m_v_out = b_state ? m_vHigh : 0.0;
+	
+	LogicIn::setElementSet(c);
+}
+
+
+void LogicOut::setOutputHighConductance( double g )
+{
+	m_bOutputHighConductanceConst = true;
+	if ( g == m_gHigh )
+		return;
+	m_gHigh = g;
+	configChanged();
+}
+
+
+void LogicOut::setOutputLowConductance( double g )
+{
+	m_bOutputLowConductanceConst = true;
+	if ( g == m_gLow )
+		return;
+	m_gLow = g;
+	configChanged();
+}
+
+
+void LogicOut::setOutputHighVoltage( double v )
+{
+	m_bOutputHighVoltageConst = true;
+	if ( v == m_vHigh )
+		return;
+	m_vHigh = v;
+	configChanged();
+}
+
+
+void LogicOut::setLogic( LogicConfig config )
+{
+	m_config = config;
+	
+	if (!m_bOutputHighConductanceConst)
+		m_gHigh = 1.0/config.highImpedance;
+	
+	if (!m_bOutputLowConductanceConst)
+		m_gLow = (config.lowImpedance == 0.0) ? 0.0 : 1.0/config.lowImpedance;
+	
+	if (!m_bOutputHighVoltageConst)
+		m_vHigh = config.output;
+	
+	configChanged();
+}
+
+
+void LogicOut::configChanged()
+{
+	if (m_bUseLogicChain)
+		return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	// Re-add the DC stuff using the new values
+	
+	m_old_g_out = m_g_out;
+	m_old_v_out = m_v_out;
+	
+	// NOTE Make sure that the next two lines are the same as those in setElementSet and setHigh
+	m_g_out = b_state ? m_gHigh : m_gLow;
+	m_v_out = b_state ? m_vHigh : 0.0;
+	
+	add_initial_dc();
+	
+	m_old_g_out = 0.;
+	m_old_v_out = 0.;
+	
+	check();
+}
+
+
+void LogicOut::add_map()
+{
+	if (!b_status) return;
+	
+	p_A->setUse( p_cnode[0]->n(), p_cnode[0]->n(), Map::et_variable, false );
+}
+
+
+void LogicOut::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_g( 0, 0 ) += m_g_out-m_old_g_out;
+	b_i( 0 ) += m_g_out*m_v_out-m_old_g_out*m_old_v_out;
+}
+
+void LogicOut::updateCurrents()
+{
+	if (m_bUseLogicChain)
+	{
+		m_cnodeI[0] = 0.;
+		return;
+	}
+	if (!b_status) {
+		return;
+	}
+	m_cnodeI[0] = (p_cnode[0]->v-m_v_out)*m_g_out;
+}
+
+void LogicOut::setHigh( bool high )
+{
+	if ( high == b_state )
+		return;
+	
+	if (m_bUseLogicChain)
+	{
+		b_state = high;
+		
+		for ( LogicIn * logic = this; logic; logic = logic->nextLogic() )
+			logic->setLastState(high);
+		
+		if (m_bCanAddChanged)
+		{
+			m_pSimulator->addChangedLogic(this);
+			m_bCanAddChanged = false;
+		}
+	
+		return;
+	}
+	
+	m_old_g_out = m_g_out;
+	m_old_v_out = m_v_out;
+	
+	// NOTE Make sure that the next two lines are the same as those in setElementSet and setLogic
+	m_g_out = high ? m_gHigh : m_gLow;
+	m_v_out = high ? m_vHigh : 0.0;
+	
+	add_initial_dc();
+	
+	m_old_g_out = 0.;
+	m_old_v_out = 0.;
+	
+	b_state = high;
+	
+	if ( p_eSet && p_eSet->circuit()->canAddChanged() )
+	{
+		m_pSimulator->addChangedCircuit( p_eSet->circuit() );
+		p_eSet->circuit()->setCanAddChanged(false);
+	}
+}
+
diff --git a/src/electronics/simulation/logic.h b/src/electronics/simulation/logic.h
new file mode 100644
index 0000000..be8374f
--- /dev/null
+++ b/src/electronics/simulation/logic.h
@@ -0,0 +1,211 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef LOGIC_H
+#define LOGIC_H
+
+#include "element.h"
+
+#include <qguardedptr.h>
+#include <qvaluelist.h>
+
+class Component;
+class Pin;
+class Simulator;
+
+typedef QValueList<QGuardedPtr<Pin> > PinList;
+
+typedef struct
+{
+	float risingTrigger; // Trigger on rising edge
+	float fallingTrigger; // Trigger on falling edge
+	float output; // Output voltage
+	float highImpedance; // Output impedance when high
+	float lowImpedance; // Output impedance when low
+} LogicConfig;
+
+class CallbackClass {};
+typedef void(CallbackClass::*CallbackPtr)( bool isHigh );
+
+/**
+Use this class for Logic Inputs - this will have infinite impedance.
+Use isHigh() will return whether the voltage level at the pin
+is high than the predetermined voltage threshold, and setHigh() will make the
+output high/low, also according to the predetermined logic type / voltages.
+
+@short Boolean Logic input
+*/
+class LogicIn : public Element
+{
+	public:
+		LogicIn( LogicConfig config );
+		virtual ~LogicIn();
+	
+		virtual Type type() const { return Element_LogicIn; }
+		virtual void setElementSet( ElementSet *c );
+	
+		/**
+		 * Set logic values from the LogicConfig.
+		 */
+		virtual void setLogic( LogicConfig config );
+		/**
+		 * Check if the input state has changed, to see if we need to callback.
+		 */
+		void check();
+		/**
+		 * Returns whether the pin is 'high', as defined for the logic type
+		 * Note: this is defined as the voltage on the pin, as opposed to what the 
+		 * state was set to (the two are not necessarily the same).
+		 */
+		bool isHigh() const { return m_bLastState; }
+		/**
+		 * When the logic state on this LogicIn changes, the function passed in this
+		 * function will be called. At most one Callback can be added per LogicIn.
+		 */
+		void setCallback( CallbackClass * object, CallbackPtr func );
+		/**
+		 * Reads the LogicConfig values in from KTLConfig, and returns them in a
+		 * nice object form.
+		 */
+		static LogicConfig getConfig();
+		/**
+		 * If this belongs to a logic chain, then this will be called from the chain.
+		 */
+		void setLastState( bool state ) { m_bLastState = state; }
+		/**
+		 * Returns a pointer to the next LogicIn in the chain.
+		 */
+		LogicIn * nextLogic() const { return m_pNextLogic; }
+		/**
+		 * Sets the next LogicIn in the chain.
+		 */
+		void setNextLogic( LogicIn * next ) { m_pNextLogic = next; }
+		/**
+		 * Calls the callback function, if there is one.
+		 */
+		void callCallback()
+		{
+			if (m_pCallbackFunction)
+				(m_pCallbackObject->*m_pCallbackFunction)(m_bLastState);
+		}
+	
+	protected:
+		virtual void updateCurrents();
+		virtual void add_initial_dc();
+	
+		CallbackPtr m_pCallbackFunction;
+		CallbackClass * m_pCallbackObject;
+		bool m_bLastState;
+		LogicIn * m_pNextLogic;
+		LogicConfig m_config;
+};
+
+
+/**
+@short Logic output/input
+*/
+class LogicOut : public LogicIn
+{
+	public:
+		LogicOut( LogicConfig config, bool _high );
+		virtual ~LogicOut();
+	
+		virtual void setLogic( LogicConfig config );
+		virtual void setElementSet( ElementSet *c );
+		virtual void add_map();
+		virtual Type type() const { return Element_LogicOut; }
+	
+		/**
+		 * Call this function to override the logic-high output impedance as set by
+		 * the user. Once set, the impedance will not be changed by the user
+		 * updating the config; only by subsequent calls to this function.
+		 */
+		void setOutputHighConductance( double g );
+		/**
+		 * Call this function to override the logic-low output impedance as set by
+		 * the user. Once set, the impedance will not be changed by the user
+		 * updating the config; only by subsequent calls to this function.
+		 */
+		void setOutputLowConductance( double g );
+		/**
+		 * Call this function to override the logic out voltage as set by the
+		 * user. Once set, the impedance will not be changed by the user
+		 * updating the config; only by subsequent calls to this function.
+		 */
+		void setOutputHighVoltage( double v );
+		/**
+		 * Returns the voltage that this will output when high.
+		 */
+		double outputHighVoltage() const { return m_vHigh; }
+		/**
+		 * Sets the pin to be high/low
+		 */
+		void setHigh( bool high );
+		/**
+		 * @returns the state that this is outputting (regardless of voltage level on logic)
+		 */
+		bool outputState() const { return b_state; }
+		/**
+		 * Set whether or not this LogicOut is the head of a LogicChain (controls
+		 * itself and a bunch of LogicIns).
+		 */
+		void setUseLogicChain( bool use );
+		/**
+		 * When a LogicOut configured as the start of a LogicChain changes start, it
+		 * appends a pointer to itself to the list of change LogicOut, starting from
+		 * the Simulator. This functions enables appending the next changed LogicOut
+		 * to this one.
+		 */
+		void setNextChanged( LogicOut * logicOut, unsigned char chain ) { m_pNextChanged[chain] = logicOut; }
+		/**
+		 * To avoid a pointer to this LogicOut being added twice in one
+		 * iteration due to the state changing twice, this LogicOut sets an
+		 * added flag to true after adding it to the list of changed. The flag
+		 * must be reset to false with this function (done by Simulator).
+		 */
+		void setCanAddChanged( bool canAdd ) { m_bCanAddChanged = canAdd; }
+		/**
+		 * Returns the next LogicOut that has changed, when configured as the start
+		 * of a LogicChain.
+		 * @see setNextChanged
+		 */
+		LogicOut * nextChanged( unsigned char chain ) const { return m_pNextChanged[chain]; }
+		PinList pinList;
+		PinList::iterator pinListBegin;
+		PinList::iterator pinListEnd;
+	
+	protected:
+		void configChanged();
+		virtual void updateCurrents();
+		virtual void add_initial_dc();
+	
+		// Pre-initalized levels from config
+		double m_gHigh;
+		double m_gLow;
+		double m_vHigh;
+		
+		// Whether to use the user-defined logic values
+		bool m_bOutputHighConductanceConst;
+		bool m_bOutputLowConductanceConst;
+		bool m_bOutputHighVoltageConst;
+	
+		double m_g_out;
+		double m_v_out;
+		double m_old_g_out;
+		double m_old_v_out;
+		bool b_state;
+		bool m_bCanAddChanged;
+		LogicOut * m_pNextChanged[2];
+		Simulator * m_pSimulator;
+		bool m_bUseLogicChain;
+};
+
+#endif
+
diff --git a/src/electronics/simulation/matrix.cpp b/src/electronics/simulation/matrix.cpp
new file mode 100644
index 0000000..fb1248f
--- /dev/null
+++ b/src/electronics/simulation/matrix.cpp
@@ -0,0 +1,546 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "matrix.h"
+
+#include <kdebug.h>
+
+#include <assert.h>
+
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+/// Minimum value before an entry is deemed "zero"
+const double epsilon = 1e-50;
+
+Matrix::Matrix( uint n, uint m )
+{
+	m_n = n;
+	m_m = m;
+	m_size = m_n+m_m;
+	
+	m_mat = new matrix(m_size);
+	m_lu = new matrix(m_size);
+	m_y = new double[m_size];
+	m_inMap = new int[m_size];
+// 	m_outMap = new int[m_size];
+	m_map = new Map(m_size);
+	zero();
+}
+
+
+Matrix::~Matrix()
+{
+	delete m_map;
+	delete m_mat;
+	delete m_lu;
+	delete [] m_y;
+	delete [] m_inMap;
+// 	delete [] m_outMap;
+}
+
+
+void Matrix::zero()
+{
+	for ( uint i=0; i<m_size; i++ )
+	{
+		for ( uint j=0; j<m_size; j++ )
+		{
+			(*m_mat)[i][j] = 0.;
+			(*m_lu)[i][j] = 0.;
+		}
+		m_inMap[i] = i;
+// 		m_outMap[i] = i;
+	}
+	
+	max_k = 0;
+}
+
+
+void Matrix::setUse( const uint i, const uint j, Map::e_type type, bool big )
+{
+	m_map->setUse( i, j, type, big );
+}
+
+
+void Matrix::createMap()
+{
+	int newMap[m_size];
+	m_map->createMap(newMap);
+	for ( uint i=0; i<m_size; i++ )
+	{
+		const int nu = newMap[i];
+		if ( nu != m_inMap[i] )
+		{
+			int old = -1;
+			for ( uint j=0; j<m_size && old == -1; j++ )
+			{
+				if ( m_inMap[j] == nu ) {
+					old = j;
+				}
+			}
+			assert( old != -1 );
+			swapRows( old, i );
+		}
+	}
+}
+
+
+void Matrix::swapRows( const uint a, const uint b )
+{
+	if ( a == b ) return;
+	m_mat->swapRows( a, b );
+	
+	const int old = m_inMap[a];
+	m_inMap[a] = m_inMap[b];
+	m_inMap[b] = old;
+	
+	max_k = 0;
+}
+
+
+/*void Matrix::genOutMap()
+{
+	for ( uint i=0; i<m_size; i++ )
+	{
+		m_outMap[ m_inMap[i] ] = i;
+	}
+}*/
+
+
+void Matrix::operator=( Matrix *const m )
+{
+	for ( uint _i=0; _i<m_size; _i++ )
+	{
+		uint i = m_inMap[_i];
+		for ( uint j=0; j<m_size; j++ )
+		{
+			(*m_mat)[i][j] = m->m(i,j);
+		}
+	}
+	
+	max_k = 0;
+}
+
+void Matrix::operator+=( Matrix *const m )
+{
+	for ( uint _i=0; _i<m_size; _i++ )
+	{
+		uint i = m_inMap[_i];
+		for ( uint j=0; j<m_size; j++ )
+		{
+			(*m_mat)[i][j] += m->m(i,j);
+		}
+	}
+	
+	max_k = 0;
+}
+
+void Matrix::performLU()
+{
+// 	max_k = 0;
+	uint n = m_size;
+	if ( n == 0 ) return;
+	
+	// Copy the affected segment to LU
+	for ( uint i=max_k; i<n; i++ )
+	{
+		for ( uint j=max_k; j<n; j++ )
+		{
+			(*m_lu)[i][j] = (*m_mat)[i][j];
+		}
+	}
+	
+	// LU decompose the matrix, and store result back in matrix
+	for ( uint k=0; k<n-1; k++ )
+	{
+		double * const lu_K_K = &(*m_lu)[k][k];
+		if ( std::abs(*lu_K_K) < 1e-10 )
+		{
+			if ( *lu_K_K < 0. ) *lu_K_K = -1e-10;
+			else *lu_K_K = 1e-10;
+		}
+		for ( uint i=std::max(k,max_k)+1; i<n; i++ )
+		{
+			(*m_lu)[i][k] /= *lu_K_K;
+		}
+		for ( uint i=std::max(k,max_k)+1; i<n; i++ )
+		{
+			const double lu_I_K = (*m_lu)[i][k];
+			if ( std::abs(lu_I_K) > 1e-12 )
+			{
+				for ( uint j=std::max(k,max_k)+1; j<n; j++ )
+				{
+					(*m_lu)[i][j] -= lu_I_K*(*m_lu)[k][j];
+				}
+			}
+		}
+	}
+	
+	max_k = n;
+}
+
+void Matrix::fbSub( Vector* b )
+{
+	if ( m_size == 0 ) return;
+	
+	for ( uint i=0; i<m_size; i++ )
+	{
+		m_y[m_inMap[i]] = (*b)[i];
+	}
+	
+	// Forward substitution
+	for ( uint i=1; i<m_size; i++ )
+	{
+		double sum = 0;
+		for ( uint j=0; j<i; j++ )
+		{
+			sum += (*m_lu)[i][j]*m_y[j];
+		}
+		m_y[i] -= sum;
+	}
+	
+	// Back substitution
+	m_y[m_size-1] /= (*m_lu)[m_size-1][m_size-1];
+	for ( int i=m_size-2; i>=0; i-- )
+	{
+		double sum = 0;
+		for ( uint j=i+1; j<m_size; j++ )
+		{
+			sum += (*m_lu)[i][j]*m_y[j];
+		}
+		m_y[i] -= sum;
+		m_y[i] /= (*m_lu)[i][i];
+	}
+	
+	for ( uint i=0; i<m_size; i++ )
+		(*b)[i] = m_y[i];
+}
+
+
+void Matrix::multiply( Vector *rhs, Vector *result )
+{
+	if ( !rhs || !result ) return;
+	result->reset();
+	for ( uint _i=0; _i<m_size; _i++ )
+	{
+		uint i = m_inMap[_i];
+		for ( uint j=0; j<m_size; j++ )
+		{
+			(*result)[_i] += (*m_mat)[i][j] * (*rhs)[j];
+		}
+	}
+}
+
+
+void Matrix::displayMatrix()
+{
+	uint n = m_size;
+	for ( uint _i=0; _i<n; _i++ )
+	{
+		uint i = m_inMap[_i];
+		for ( uint j=0; j<n; j++ )
+		{
+			if ( j > 0 && (*m_mat)[i][j] >= 0 ) kdDebug() << "+";
+			kdDebug() << (*m_mat)[i][j] << "("<<j<<")";
+		}
+		kdDebug()  << endl;
+	}
+}
+
+void Matrix::displayLU()
+{
+	uint n = m_size;
+	for ( uint _i=0; _i<n; _i++ )
+	{
+		uint i = m_inMap[_i];
+// 		uint i = _i;
+		for ( uint j=0; j<n; j++ )
+		{
+			if ( j > 0 && (*m_lu)[i][j] >= 0 ) std::cout << "+";
+			std::cout << (*m_lu)[i][j] << "("<<j<<")";
+		}
+		std::cout << std::endl;
+	}
+	std::cout << "m_inMap:    ";
+	for ( uint i=0; i<n; i++ )
+	{
+		std::cout << i<<"->"<<m_inMap[i]<<"  ";
+	}
+	std::cout << std::endl;
+	/*cout << "m_outMap:   ";
+	for ( uint i=0; i<n; i++ )
+	{
+		cout << i<<"->"<<m_outMap[i]<<"  ";
+	}
+	cout << endl;*/
+}
+
+
+Map::Map( const uint size )
+{
+	m_size = size;
+	m_map = new ETMap( m_size );
+	reset();
+}
+
+
+Map::~Map()
+{
+	delete m_map;
+}
+
+
+void Map::reset()
+{
+	for ( uint i=0; i<m_size; i++ )
+	{
+		for ( uint j=0; j<m_size; j++ )
+		{
+			(*m_map)[i][j] = 0;
+		}
+	}
+}
+
+
+void Map::setUse( const uint i, const uint j, Map::e_type type, bool big )
+{
+	if ( type == Map::et_none ) {
+		(*m_map)[i][j] = Map::et_none;
+	} else {
+		(*m_map)[i][j] = type | (big)?Map::et_big:0;
+	}
+}
+
+
+void Map::createMap( int *map )
+{
+	assert(map);
+	
+	// In this function, the passes through that we make want to be done from
+	// top left to bottom right, to minimise fill-in
+	
+	// available[i] is true if an external-row can be mapped to internal-row "i"
+	// map[i] gives the internal-row for external-row i
+	bool available[m_size];
+	for ( uint i=0; i<m_size; i++ )
+	{
+		available[i] = true;
+		map[i] = -1;
+	}
+	
+	// This loop looks through columns and rows to find any swaps that are necessary
+	// (e.g. only one matrix-element in that row/column), and if no necessary swaps
+	// were found, then it will swap two rows according to criteria given below
+	bool badMap = false;
+	bool changed;
+	do
+	{
+		changed = false;
+		
+		// Pass through columns
+		int E,N;
+		uint highest = 0;
+		for ( uint j=0; j<m_size; j++ )
+		{
+			if ( map[j] == -1 ) // If we haven't mapped this column yet
+			{
+				int count = 0; // Number of "spare" elements
+				int element; // Last element that is "spare", only applicable if count=1
+				for ( uint i=0; i<m_size; i++ )
+				{
+					if ( available[i] && (*m_map)[i][j] )
+					{
+						count++;
+						element = i;
+					}
+				}
+				if ( count == 0 ) {
+					badMap = true;
+				}
+				else if ( count == 1 )
+				{
+					const uint newType = (*m_map)[element][j];
+					if ( typeCmp( newType, highest) )
+					{
+						E=element;
+						N=j;
+						highest=newType;
+					}
+				}
+			}
+		}
+		// Pass through rows
+		for ( uint i=0; i<m_size; i++ )
+		{
+			if ( map[i] == -1 ) // If we haven't mapped this row yet
+			{
+				int count = 0; // Number of "spare" elements
+				int element; // Last element that is "spare", only applicable if count=1
+				for ( uint j=0; j<m_size; j++ )
+				{
+					if ( available[j] && (*m_map)[i][j] )
+					{
+						count++;
+						element = j;
+					}
+				}
+				if ( count == 0 ) {
+					badMap = true;
+				}
+				else if ( count == 1 )
+				{
+					const uint newType = (*m_map)[i][element];
+					if ( typeCmp( newType, highest) )
+					{
+						E=element;
+						N=i;
+						highest=newType;
+					}
+				}
+			}
+		}
+		if (highest)
+		{
+			available[E] = false;
+			map[N] = E;
+			changed = true;
+		}
+		if (!changed)
+		{
+			int next = -1; // next is the row to mapped to (interally)
+			uint j=0;
+			
+			/// TODO We want to change this search to one that finds a swap, taking into acocunt the priorities given below
+			while ( next == -1 && j<m_size )
+			{
+				if ( available[j] ) next=j;
+				j++;
+			}
+			uint i=0;
+			while ( i<m_size && map[i] != -1 ) i++;
+			if ( next != -1 && i < m_size )
+			{
+				available[next] = false;
+				map[i] = next;
+				changed = true;
+			}
+		}
+	} while (changed);
+	
+	if (badMap)
+	{
+// 		cerr << "Map::createMap: unable to create decent mapping; do not trust the matrix, Neo!"<<endl;
+	}
+	
+	for ( int i = 0; i < int(m_size); ++i )
+	{
+		assert( map[i] >= 0 && map[i] < int(m_size) );
+	}
+	
+	// Ignore this, for now:
+	
+	// Now, we want to order the matrix, with the following priorities:
+	//	(1) How often values change
+	//	(2) How few values there are
+	//	(3) How large the values are
+	// For each value in the column, 
+}
+
+
+bool Map::typeCmp( const uint t1, const uint t2 )
+{
+	if (!t2) return true;
+	if (!t1) return false;
+	
+	int t1_score = 1;
+	if		( t1 | Map::et_constant )	t1_score += 64;
+	else if ( t1 | Map::et_stable )		t1_score += 16;
+	else if ( t1 | Map::et_variable )	t1_score += 4;
+	
+	int t2_score = 1;
+	if		( t2 | Map::et_constant )	t2_score += 64;
+	else if ( t2 | Map::et_stable )		t2_score += 16;
+	else if ( t2 | Map::et_variable )	t2_score += 4;
+	
+	if ( t1 | Map::et_big ) t1_score *= 2;
+	if ( t2 | Map::et_big ) t2_score *= 2;
+	
+	return ( t1_score >= t2_score );
+}
+
+
+Matrix22::Matrix22()
+{
+	reset();
+}
+
+bool Matrix22::solve()
+{
+	const double old_x1 = m_x1;
+	const double old_x2 = m_x2;
+	
+	const bool e11 = std::abs((m_a11))<epsilon;
+	const bool e12 = std::abs((m_a12))<epsilon;
+	const bool e21 = std::abs((m_a21))<epsilon;
+	const bool e22 = std::abs((m_a22))<epsilon;
+	
+	if (e11)
+	{
+		if ( e12||e21 )
+			return false;
+		m_x2 = m_b1/m_a12;
+		m_x1 = (m_b2-(m_a22*m_x2))/m_a21;
+	}
+	else if (e12)
+	{
+		if ( e11||e22 )
+			return false;
+		m_x1 = m_b1/m_a11;
+		m_x2 = (m_b2-(m_a21*m_x1))/m_a22;
+	}
+	else if (e21)
+	{
+		if ( e11||e22 )
+			return false;
+		m_x2 = m_b2/m_a22;
+		m_x1 = (m_b1-(m_a12*m_x2))/m_a11;
+	}
+	else if (e22)
+	{
+		if ( e12||e21 )
+			return false;
+		m_x1 = m_b2/m_a21;
+		m_x2 = (m_b1-(m_a11*m_x1))/m_a12;
+	}
+	else
+	{
+		m_x2 = (m_b2-(m_a21*m_b1/m_a11))/(m_a22-(m_a21*m_a12/m_a11));
+		m_x1 = (m_b1-(m_a12*m_x2))/m_a11;
+	}
+	if ( !std::isfinite(m_x1) || !std::isfinite(m_x2) )
+	{
+		m_x1 = old_x1;
+		m_x2 = old_x2;
+		return false;
+	}
+	return true;
+}
+
+void Matrix22::reset()
+{
+	m_a11=m_a12=m_a21=m_a22=0.;
+	m_b1=m_b2=0.;
+	m_x1=m_x2=0.;
+}
+
+
+
diff --git a/src/electronics/simulation/matrix.h b/src/electronics/simulation/matrix.h
new file mode 100644
index 0000000..4c3e518
--- /dev/null
+++ b/src/electronics/simulation/matrix.h
@@ -0,0 +1,248 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include "vec.h"
+
+#include <vector>
+
+class Map;
+
+typedef std::vector<std::vector<uint> > ETMap;
+
+/**
+@short Handles row-wise permutation of matrixes
+*/
+class Map
+{
+public:
+	enum e_type
+	{
+		et_none = 0x0, // Default value, none
+		et_constant = 0x1, // Never changes value in lifetime of matrix
+		et_stable = 0x2, // Value changes occasionally, e.g. user changing resistance value
+		et_variable = 0x3, // Rate of changing is unknown, probably average - e.g. logic
+		et_unstable = 0x4, // Value changes practically every call of performLU
+		et_big = 0x8 // Ignore this :-) It is used to OR with one of the others, but it should only ever be accessed by the class
+	};
+	Map( const uint size );
+	~Map();
+
+	void setUse( const uint i, const uint j, Map::e_type type, bool big );
+	/**
+	 * Generates an optimal permutation, returned in the given array
+	 */
+	void createMap( int *map );
+	/**
+	 * Resets the info to a blank pattern
+	 */
+	void reset();
+	
+protected:
+	/**
+	 * Compares two "types", returning true if t1 >= t2, else false
+	 */
+	bool typeCmp( const uint t1, const uint t2 );
+	
+private:
+	ETMap *m_map;
+	uint m_size;
+};
+
+/**
+This class performs matrix storage, lu decomposition, forward and backward
+substitution, and a few other useful operations. Steps in using class:
+(1) Create an instance of this class with the correct size
+(2) Define the matrix pattern as neccessary:
+	(1) Call zero (unnecessary after initial ceration) to reset the pattern
+		& matrix
+	(2) Call setUse to set the use of each element in the matrix
+	(3) Call createMap to generate the row-wise permutation mapping for use
+		in partial pivoting
+(3) Add the values to the matrix
+(4) Call performLU, and get the results with fbSub
+(5) Repeat 2, 3, 4 or 5 as necessary.
+@todo We need to allow createMap to work while the matrix has already been initalised
+@short Matrix manipulation class tailored for circuit equations
+@author David Saxton
+*/
+class Matrix
+{
+public:
+	/**
+	 * Creates a size x size square matrix m, with all values zero,
+	 * and a right side vector x of size m+n
+	 */
+	Matrix( uint n, uint m );
+	~Matrix();
+	/**
+	 * Sets all elements to zero
+	 */
+	void zero();
+	/** 
+	 * Sets the type of (matrix-) element at the position i, j.
+	 * @param type Describes how often the value changes
+	 * @param big Set this true if the value is likely to be >= 1, else false
+	 */
+	void setUse( const uint i, const uint j, Map::e_type type, bool big );
+	void setUse_b( const uint i, const uint j, Map::e_type type, bool big )
+	{
+		setUse( i, j+m_n, type, big );
+	}
+	void setUse_c( const uint i, const uint j, Map::e_type type, bool big )
+	{
+		setUse( i+m_n, j, type, big );
+	}
+	void setUse_d( const uint i, const uint j, Map::e_type type, bool big )
+	{
+		setUse( i+m_n, j+m_n, type, big );
+	}
+	/**
+	 * Generates the row-wise permutation mapping from the values set by setUse
+	 */
+	void createMap();
+	/**
+	 * Returns true if the matrix is changed since last calling performLU()
+	 * - i.e. if we do need to call performLU again.
+	 */
+	inline bool isChanged() const { return max_k < m_size; }
+	/**
+	 * Performs LU decomposition. Going along the rows,
+	 * the value of the decomposed LU matrix depends only on
+	 * the previous values.
+	 */
+	void performLU();
+	/**
+	 * Applies the right side vector (x) to the decomposed matrix,
+	 * with the solution returned in x.
+	 */
+	void fbSub( Vector* x );
+	/**
+	 * Prints the matrix to stdout
+	 */
+	void displayMatrix();
+	/**
+	 * Prints the LU-decomposed matrix to stdout
+	 */
+	void displayLU();
+	/**
+	 * Sets the element matrix at row i, col j to value x
+	 */
+	double& g( uint i, uint j )
+	{
+		i = m_inMap[i];
+		if ( i<max_k ) max_k=i;
+		if ( j<max_k ) max_k=j;
+		
+		// I think I need the next line...
+		if ( max_k>0 ) max_k--;
+		
+		return (*m_mat)[i][j];
+	}
+	double& b( uint i, uint j ) { return g( i, j+m_n ); }
+	double& c( uint i, uint j ) { return g( i+m_n, j ); }
+	double& d( uint i, uint j ) { return g( i+m_n, j+m_n ); }
+	const double g( uint i, uint j ) const { return (*m_mat)[m_inMap[i]][j]; }
+	const double b( uint i, uint j ) const { return g( i, j+m_n ); }
+	const double c( uint i, uint j ) const { return g( i+m_n, j ); }
+	const double d( uint i, uint j ) const { return g( i+m_n, j+m_n ); }
+	/**
+	 * Returns the value of matrix at row i, col j.
+	 */
+	const double m( uint i, uint j ) const
+	{
+		return (*m_mat)[m_inMap[i]][j];
+	}
+	/**
+	 * Multiplies this matrix by the Vector rhs, and places the result
+	 * in the vector pointed to by result. Will fail if wrong size.
+	 */
+	void multiply( Vector *rhs, Vector *result );
+	/**
+	 * Sets the values of this matrix to that of the given matrix
+	 */
+	void operator=( Matrix *const m );
+	/**
+	 * Adds the values of the given matrix to this matrix
+	 */
+	void operator+=( Matrix *const m );
+	
+private:
+	/**
+	 * Swaps around the rows in the (a) the matrix; and (b) the mappings
+	 */
+	void swapRows( const uint a, const uint b );
+	
+// 	// Generates m_outMap from m_inMap
+// 	void genOutMap();
+	
+	uint m_n;
+	uint m_m;
+	uint m_size;
+	uint max_k;
+	
+	int *m_inMap; // Rowwise permutation mapping from external reference to internal storage
+// 	int *m_outMap; // Opposite of m_inMap
+	
+	matrix *m_mat;
+	matrix *m_lu;
+	double *m_y; // Avoids recreating it lots of times
+	Map *m_map;
+};
+
+
+/**
+This class provides a very simple, lightweight, 2x2 matrix solver.
+It's fast and reliable. Set the values for the entries of A and b:
+
+A x = b
+
+call solve() (which returns true if successful - i.e. exactly one solution to the
+matrix), and get the values of x with the appropriate functions.
+
+@short 2x2 Matrix
+@author David Saxton
+*/
+class Matrix22
+{
+public:
+	Matrix22();
+	
+	double &a11() { return m_a11; }
+	double &a12() { return m_a12; }
+	double &a21() { return m_a21; }
+	double &a22() { return m_a22; }
+	
+	double &b1() { return m_b1; }
+	double &b2() { return m_b2; }
+	
+	/**
+	 * Solve the matrix. Returns true if successful (i.e. non-singular), else
+	 * false. Get the solution with x1() and x2().
+	 */
+	bool solve();
+	/**
+	 * Resets all entries to zero
+	 */
+	void reset();
+	
+	double x1() const { return m_x1; }
+	double x2() const { return m_x2; }
+	
+private:
+	double m_a11, m_a12, m_a21, m_a22;
+	double m_b1, m_b2;
+	double m_x1, m_x2;
+};
+
+
+#endif
diff --git a/src/electronics/simulation/nonlinear.cpp b/src/electronics/simulation/nonlinear.cpp
new file mode 100644
index 0000000..c975f16
--- /dev/null
+++ b/src/electronics/simulation/nonlinear.cpp
@@ -0,0 +1,97 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "matrix.h"
+#include "nonlinear.h"
+
+#include <cmath>
+using namespace std;
+
+const double KTL_MAX_DOUBLE = 1.7976931348623157e+308; ///< 7fefffff ffffffff
+const int KTL_MAX_EXPONENT = int( log( KTL_MAX_DOUBLE ) );
+
+
+NonLinear::NonLinear()
+	: Element()
+{
+}
+
+
+#ifndef MIN
+# define MIN(x,y) (((x) < (y)) ? (x) : (y))
+#endif
+
+
+// The function computes the exponential pn-junction current.
+double NonLinear::diodeCurrent( double v, double I_S, double Vte ) const
+{
+	return I_S * (exp( MIN( v / Vte, KTL_MAX_EXPONENT ) ) - 1);
+}
+
+
+
+double NonLinear::diodeConductance( double v, double I_S, double Vte ) const
+{
+	return I_S * exp( MIN( v / Vte, KTL_MAX_EXPONENT ) ) / Vte;
+}
+
+
+
+double NonLinear::diodeVoltage( double V, double V_prev, double V_T, double Vcrit ) const
+{
+	if ( V > Vcrit && fabs( V - V_prev ) > 2 * V_T )
+	{
+		if ( V_prev > 0 )
+		{
+			double arg = (V - V_prev) / V_T;
+			if (arg > 0)
+				V = V_prev + V_T * (2 + log( arg - 2 ));
+			else
+				V = V_prev - V_T * (2 + log( 2 - arg ));
+		}
+		else
+			V = (V_prev < 0) ? (V_T * log (V / V_T)) : Vcrit;
+	}
+	else
+	{
+		if ( V < 0 )
+		{
+			double arg = (V_prev > 0) ? (-1 - V_prev) : (2 * V_prev - 1);
+			if (V < arg)
+				V = arg;
+		}
+	}
+	return V;
+}
+
+
+double NonLinear::diodeCriticalVoltage( double I_S, double V_Te ) const
+{
+	return V_Te * log( V_Te / M_SQRT2 / I_S );
+}
+
+
+void NonLinear::diodeJunction( double V, double I_S, double V_Te, double * I, double * g ) const
+{
+	if (V < -3 * V_Te)
+	{
+		double a = 3 * V_Te / (V * M_E);
+		a = a * a * a;
+		*I = -I_S * (1 + a);
+		*g = +I_S * 3 * a / V;
+	}
+	else
+	{ 
+		double e = exp( MIN( V / V_Te, KTL_MAX_EXPONENT ) );
+		*I = I_S * (e - 1);
+		*g = I_S * e / V_Te;
+	}
+}
+
diff --git a/src/electronics/simulation/nonlinear.h b/src/electronics/simulation/nonlinear.h
new file mode 100644
index 0000000..422f3a4
--- /dev/null
+++ b/src/electronics/simulation/nonlinear.h
@@ -0,0 +1,48 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2005 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef NONLINEAR_H
+#define NONLINEAR_H
+
+#include "element.h"
+
+/**
+@short Represents a non-linear circuit element (such as a diode)
+@author David Saxton
+*/
+class NonLinear : public Element
+{
+	public:
+		NonLinear();
+	
+		virtual bool isNonLinear() { return true; }
+		/**
+		 * Newton-Raphson iteration: Update equation system.
+		 */
+		virtual void update_dc() = 0;
+		
+	protected:
+		double diodeCurrent( double v, double I_S, double Vte ) const;
+		/**
+		 * Conductance of the diode - the derivative of Schockley's
+		 * approximation.
+		 */
+		double diodeConductance( double v, double I_S, double Vte ) const;
+		/**
+		 * Limits the diode voltage to prevent divergence in the nonlinear
+		 * iterations.
+		 */
+		double diodeVoltage( double v, double V_prev, double Vt, double V_crit ) const;
+		void diodeJunction( double v, double I_S, double Vte, double * I, double * g ) const;
+		
+		double diodeCriticalVoltage( double I_S, double Vte ) const;
+};
+
+#endif
diff --git a/src/electronics/simulation/opamp.cpp b/src/electronics/simulation/opamp.cpp
new file mode 100644
index 0000000..45fbf02
--- /dev/null
+++ b/src/electronics/simulation/opamp.cpp
@@ -0,0 +1,77 @@
+/***************************************************************************
+ *   Copyright (C) 2005 by David Saxton                                    *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "opamp.h"
+
+OpAmp::OpAmp()
+	: Element::Element()
+{
+	m_numCBranches = 1;
+	m_numCNodes = 3;
+}
+
+
+OpAmp::~OpAmp()
+{
+}
+
+
+void OpAmp::add_map()
+{
+	if (!b_status)
+		return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		// Non-inverting input
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+	
+	if ( !p_cnode[2]->isGround )
+	{
+		// Inverting input
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[2]->n(), Map::et_constant, true );
+	}
+	
+	if ( !p_cnode[1]->isGround )
+	{
+		// Output
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+	}
+}
+
+
+void OpAmp::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	// Non-inverting input
+	A_c( 0, 0 ) = 1;
+	
+	// Inverting input
+	A_c( 0, 2 ) = -1;
+	
+	// Output
+	A_b( 1, 0 ) = 1;
+}
+
+
+void OpAmp::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[0] = m_cnodeI[2] = 0.0;
+	m_cnodeI[1] = p_cbranch[0]->i;
+}
+
+
+
diff --git a/src/electronics/simulation/opamp.h b/src/electronics/simulation/opamp.h
new file mode 100644
index 0000000..22fe843
--- /dev/null
+++ b/src/electronics/simulation/opamp.h
@@ -0,0 +1,36 @@
+/***************************************************************************
+ *   Copyright (C) 2005 by David Saxton                                    *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef OPAMP_H
+#define OPAMP_H
+
+#include <element.h>
+
+/**
+node 0: non-inverting input
+node 1: output
+node 2: inverting input
+@author David Saxton
+*/
+class OpAmp : public Element
+{
+	public:
+		OpAmp();
+		virtual ~OpAmp();
+	
+		virtual Type type() const { return Element_OpAmp; }
+		virtual void add_map();
+	
+	protected:
+		virtual void updateCurrents();
+		virtual void add_initial_dc();
+};
+
+#endif
diff --git a/src/electronics/simulation/reactive.cpp b/src/electronics/simulation/reactive.cpp
new file mode 100644
index 0000000..83fcfd4
--- /dev/null
+++ b/src/electronics/simulation/reactive.cpp
@@ -0,0 +1,36 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+
+#include "reactive.h"
+
+Reactive::Reactive( const double delta )
+	: Element()
+{
+	m_delta = delta;
+}
+
+
+Reactive::~Reactive()
+{
+}
+
+
+void Reactive::setDelta( double delta )
+{
+	m_delta = delta;
+	updateStatus();
+}
+
+
+bool Reactive::updateStatus()
+{
+	return Element::updateStatus();
+}
diff --git a/src/electronics/simulation/reactive.h b/src/electronics/simulation/reactive.h
new file mode 100644
index 0000000..1142b34
--- /dev/null
+++ b/src/electronics/simulation/reactive.h
@@ -0,0 +1,42 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef REACTIVE_H
+#define REACTIVE_H
+
+#include "element.h"
+
+/**
+@short Represents a reactive element (such as a capacitor)
+@author David Saxton
+*/
+class Reactive : public Element
+{
+public:
+	Reactive( const double delta );
+	virtual ~Reactive();
+	
+	virtual bool isReactive() { return true; }
+	/**
+	 * Call this function to set the time period (in seconds)
+	 */
+	void setDelta( double delta );
+	/**
+	 * Called on every time step for the element to update itself
+	 */
+	virtual void time_step() = 0;
+	
+protected:
+	virtual bool updateStatus();
+	
+	double m_delta; // Delta time interval
+};
+
+#endif
diff --git a/src/electronics/simulation/resistance.cpp b/src/electronics/simulation/resistance.cpp
new file mode 100644
index 0000000..9c2d25d
--- /dev/null
+++ b/src/electronics/simulation/resistance.cpp
@@ -0,0 +1,95 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "resistance.h"
+
+// #include <kdebug.h>
+
+Resistance::Resistance( const double resistance )
+	: Element::Element()
+{
+	m_g = resistance < 1e-9 ? 1e9 : 1./resistance;
+	m_numCNodes = 2;
+// 	kdDebug() << k_funcinfo << endl;
+}
+
+Resistance::~Resistance()
+{
+// 	kdDebug() << k_funcinfo << endl;
+}
+
+void Resistance::setConductance( const double g )
+{
+	if ( g == m_g )
+		return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	// Remove old resistance
+	m_g = -m_g;
+	add_initial_dc();
+
+	m_g = g;
+	add_initial_dc();
+}
+
+
+void Resistance::setResistance( const double r )
+{
+	setConductance( r < 1e-9 ? 1e9 : 1./r );
+}
+
+
+void Resistance::add_map()
+{
+	if (!b_status)
+		return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[0]->n(), Map::et_stable, false );
+	}
+	if ( !p_cnode[1]->isGround ) {
+		p_A->setUse( p_cnode[1]->n(), p_cnode[1]->n(), Map::et_stable, false );
+	}
+	
+	if ( !p_cnode[0]->isGround && !p_cnode[1]->isGround )
+	{
+		p_A->setUse( p_cnode[0]->n(), p_cnode[1]->n(), Map::et_stable, false );
+		p_A->setUse( p_cnode[1]->n(), p_cnode[0]->n(), Map::et_stable, false );
+	}
+}
+
+
+void Resistance::add_initial_dc()
+{
+	if (!b_status) return;
+	
+	A_g( 0, 0 ) += m_g;
+	A_g( 1, 1 ) += m_g;
+	A_g( 0, 1 ) -= m_g;
+	A_g( 1, 0 ) -= m_g;
+}
+
+
+void Resistance::updateCurrents()
+{
+	if (!b_status) return;
+	const double v=p_cnode[0]->v-p_cnode[1]->v; 
+	m_cnodeI[1] = v*m_g;
+	m_cnodeI[0] = -m_cnodeI[1];
+}
+
+
+
+
diff --git a/src/electronics/simulation/resistance.h b/src/electronics/simulation/resistance.h
new file mode 100644
index 0000000..7e5a62e
--- /dev/null
+++ b/src/electronics/simulation/resistance.h
@@ -0,0 +1,43 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef RESISTANCE_H
+#define RESISTANCE_H
+
+#include "element.h"
+
+/**
+@short Resistance
+@author David saxton
+*/
+class Resistance : public Element
+{
+public:
+	Resistance( const double resistance );
+	virtual ~Resistance();
+	
+	virtual Type type() const { return Element_Resistance; }
+	
+	void setConductance( const double g );
+	void setResistance( const double r );
+	
+	double resistance() { return 1/m_g; }
+	double conductance() { return m_g; }
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_g; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/vccs.cpp b/src/electronics/simulation/vccs.cpp
new file mode 100644
index 0000000..7ff1bc1
--- /dev/null
+++ b/src/electronics/simulation/vccs.cpp
@@ -0,0 +1,93 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "vccs.h"
+
+VCCS::VCCS( const double gain )
+	: Element::Element()
+{
+	m_g = gain;
+	m_numCBranches = 1;
+	m_numCNodes = 4;
+}
+
+
+VCCS::~VCCS()
+{
+}
+
+
+void VCCS::setGain( const double g )
+{
+	if ( g == m_g ) return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	// Remove old values
+	m_g = -m_g;
+	add_initial_dc();
+	
+	// Add new values
+	m_g = g;
+	add_initial_dc();
+}
+
+
+void VCCS::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_stable, false );
+	}
+	if ( !p_cnode[1]->isGround  )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_stable, false );
+	}
+	if ( !p_cnode[2]->isGround  )
+	{
+		p_A->setUse_b( p_cnode[2]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[2]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[3]->isGround  )
+	{
+		p_A->setUse_b( p_cnode[3]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[3]->n(), Map::et_constant, true );
+	}
+}
+
+
+void VCCS::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_g( 2, 0 ) += m_g;
+	A_g( 3, 0 ) -= m_g;
+	A_g( 2, 1 ) -= m_g;
+	A_g( 3, 1 ) += m_g;
+}
+
+
+void VCCS::updateCurrents()
+{
+	if (!b_status)
+		return;
+	
+	m_cnodeI[0] = m_cnodeI[1] = 0.;
+	m_cnodeI[3] = (p_cnode[0]->v-p_cnode[1]->v)*m_g;
+	m_cnodeI[2] = -m_cnodeI[3];
+}
+
+
diff --git a/src/electronics/simulation/vccs.h b/src/electronics/simulation/vccs.h
new file mode 100644
index 0000000..26f1101
--- /dev/null
+++ b/src/electronics/simulation/vccs.h
@@ -0,0 +1,40 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VCCS_H
+#define VCCS_H
+
+#include "element.h"
+
+/**
+CNodes n0 and n1 are used for the voltage control.
+CNodes n2 and n3 are used for the current output.
+@short Voltage Controlled Current Source
+@author David Saxton
+*/
+class VCCS : public Element
+{
+public:
+	VCCS( const double gain );
+	virtual ~VCCS();
+	
+	virtual Type type() const { return Element_VCCS; }
+	void setGain( const double g );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_g; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/vcvs.cpp b/src/electronics/simulation/vcvs.cpp
new file mode 100644
index 0000000..68604df
--- /dev/null
+++ b/src/electronics/simulation/vcvs.cpp
@@ -0,0 +1,93 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "elementset.h"
+#include "matrix.h"
+#include "vcvs.h"
+
+VCVS::VCVS( const double gain )
+	: Element::Element()
+{
+	m_g = gain;
+	m_numCBranches = 1;
+	m_numCNodes = 4;
+}
+
+
+VCVS::~VCVS()
+{
+}
+
+
+void VCVS::setGain( const double g )
+{
+	if ( g == m_g )
+		return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	m_g = -m_g;
+	add_initial_dc();
+
+	m_g = g;
+	add_initial_dc();
+}
+
+
+void VCVS::add_map()
+{
+	if (!b_status)
+		return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_stable, false );
+	}
+	if ( !p_cnode[1]->isGround  )
+	{
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_stable, false );
+	}
+	if ( !p_cnode[2]->isGround  )
+	{
+		p_A->setUse_b( p_cnode[2]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[2]->n(), Map::et_constant, true );
+	}
+	if ( !p_cnode[3]->isGround  )
+	{
+		p_A->setUse_b( p_cnode[3]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[3]->n(), Map::et_constant, true );
+	}
+}
+
+
+void VCVS::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_c( 0, 0 ) -= m_g;
+	A_c( 0, 1 ) += m_g;
+	A_b( 2, 0 ) = 1;
+	A_c( 0, 2 ) = 1;
+	A_b( 3, 0 ) = -1;
+	A_c( 0, 3 ) = -1;
+}
+
+
+void VCVS::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[0] = m_cnodeI[1] = 0.;
+	m_cnodeI[3] = p_cbranch[0]->i;
+	m_cnodeI[2] = -m_cnodeI[3];
+}
+
+
diff --git a/src/electronics/simulation/vcvs.h b/src/electronics/simulation/vcvs.h
new file mode 100644
index 0000000..862dea9
--- /dev/null
+++ b/src/electronics/simulation/vcvs.h
@@ -0,0 +1,40 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VCVS_H
+#define VCVS_H
+
+#include "element.h"
+
+/**
+Voltage source between nodes c2 and c3
+Controlling voltage between nodes c0 and c1
+@short Voltage Controlled Voltage  Source
+@author David Saxton
+*/
+class VCVS : public Element
+{
+public:
+	VCVS( const double gain );
+	virtual ~VCVS();
+	
+	virtual Type type() const { return Element_VCVS; }
+	void setGain( const double g );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_g; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/vec.cpp b/src/electronics/simulation/vec.cpp
new file mode 100644
index 0000000..d45f505
--- /dev/null
+++ b/src/electronics/simulation/vec.cpp
@@ -0,0 +1,166 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "vec.h"
+
+#include <assert.h>
+#include <cmath>
+#include <string.h>
+using namespace std;
+
+Vector::Vector( const int size )
+{
+	m_size = size;
+	m_vec = new double[m_size];
+	reset();
+	b_changed = true;
+}
+
+
+Vector::~Vector()
+{
+	// Hmm...this looks like it's the correct format, although valgrind complains
+	// about improper memory use. Interesting. "delete m_vec" definitely leaks
+	// memory, so this seems like the lesser of two evils. I miss C memory allocation
+	// somtimes, with a nice simple free :p
+	delete [] m_vec;
+}
+
+
+void Vector::reset()
+{
+	for ( int i=0; i<m_size; i++ )
+	{
+		m_vec[i] = 0.;
+	}
+	b_changed = true;
+}
+
+
+void Vector::limitTo( double scaleMax, Vector * limitVector )
+{
+	assert( limitVector );
+	assert( limitVector->size() == size() );
+	
+	for ( int i = 0; i < m_size; ++i )
+	{
+		double limitAbs = std::abs( limitVector->m_vec[i] );
+		if ( limitAbs < 1e-6 )
+			limitAbs = 1e-6;
+		
+		double thisAbs = std::abs( m_vec[i] );
+		if ( thisAbs < 1e-6 )
+			thisAbs = 1e-6;
+		
+		if ( (thisAbs / limitAbs) > scaleMax )
+			m_vec[i] /= (thisAbs / limitAbs);
+		
+		else if ( (limitAbs / thisAbs) > scaleMax )
+			m_vec[i] /= (limitAbs / thisAbs);
+	}
+	b_changed = true;
+}
+
+
+void Vector::operator += ( Vector *rhs )
+{
+	if (!rhs) return;
+	for ( int i=0; i<m_size; i++ )
+	{
+		m_vec[i] += (*rhs)[i];
+	}
+	b_changed = true;
+}
+
+
+void Vector::operator -= ( Vector *rhs )
+{
+	if (!rhs) return;
+	for ( int i=0; i<m_size; i++ )
+	{
+		m_vec[i] -= (*rhs)[i];
+	}
+	b_changed = true;
+}
+
+
+void Vector::operator *=( double s )
+{
+	for ( int i=0; i<m_size; i++ )
+	{
+		m_vec[i] *= s;
+	}
+	b_changed = true;
+}
+
+	
+void Vector::operator = ( Vector& v )
+{
+	assert( size() == v.size() );
+	memcpy( m_vec, v.m_vec, m_size * sizeof( double ) );
+	b_changed = true;
+}
+
+
+void Vector::negative( Vector *rhs )
+{
+	if (!rhs) return;
+	for ( int i=0; i<m_size; i++ )
+	{
+		m_vec[i] = -(*rhs)[i];
+	}
+	b_changed = true;
+}
+
+
+double Vector::abs() const
+{
+	double s=0;
+	for ( int i=0; i<m_size; i++ )
+	{
+		s += m_vec[i]*m_vec[i];
+	}
+	return sqrt(s);
+}
+
+
+
+// matrix stuff...
+
+matrix::matrix( const uint size )
+{
+	m_size = size;
+	m_mat = new Vector*[m_size];
+	for ( uint i=0; i<m_size; ++i )
+	{
+		m_mat[i] = new Vector(m_size);
+	}
+}
+
+matrix::~matrix()
+{
+	for ( uint i=0; i<m_size; ++i )
+	{
+		delete m_mat[i];
+	}
+	delete [] m_mat;
+}
+
+
+void matrix::swapRows( const uint a, const uint b )
+{
+	if ( a == b ) return;
+	Vector *v = m_mat[a];
+	m_mat[a] = m_mat[b];
+	m_mat[b] = v;
+}
+
+
+
diff --git a/src/electronics/simulation/vec.h b/src/electronics/simulation/vec.h
new file mode 100644
index 0000000..7192bb3
--- /dev/null
+++ b/src/electronics/simulation/vec.h
@@ -0,0 +1,109 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VEC_H
+#define VEC_H
+
+typedef unsigned uint;
+
+/**
+@short Vector of doubles, faster than STL vector
+@author David Saxton
+*/
+class Vector
+{
+public:
+	Vector( const int size );
+	~Vector();
+	
+	double & operator[]( const int i ) { b_changed=true; return m_vec[i]; }
+	const double operator[]( const int i ) const { return m_vec[i]; }
+	int size() const { return m_size; }
+	/**
+	 * Resets all values to 0
+	 */
+	void reset();
+	/**
+	 * Limits the absolute value of each component of this vector to scaleMax
+	 * times bigger or smaller than the components of limitVector.
+	 */
+	void limitTo( double scaleMax, Vector * limitVector );
+	/**
+	 * Adds the Vector rhs to this
+	 */
+	void operator+=( Vector *rhs );
+	/**
+	 * Subtracts the Vector rhs from this
+	 */
+	void operator-=( Vector *rhs );
+	/**
+	 * Multiplies this Vector by the given scaler constant
+	 */
+	void operator*=( double s );
+	/**
+	 * Sets this vector equal to the given vector
+	 */
+	void operator=( Vector& v );
+	/**
+	 * Copies the negative values of the given vector to this vector.
+	 * (i.e. sets this = -rhs )
+	 */
+	void negative( Vector *rhs );
+	/**
+	 * Returns the absolute value of this vector, defined as the squareroot
+	 * of the sum of the square of each element of the vector
+	 */
+	double abs() const;
+	/**
+	 * Returns true if the vector has changed since setUnchanged was last called
+	 * Note that this will return true if the vector has just been read, due to
+	 * limitations with the [] operator.
+	 */
+	inline bool isChanged() const { return b_changed; }
+	/**
+	 * Sets the changed status to false.
+	 */
+	inline void setUnchanged() { b_changed=false; }
+
+private:
+	Vector( const Vector & );
+	Vector & operator= ( const Vector & );
+	
+	bool b_changed;
+	double *m_vec;
+	int m_size;
+};
+
+/**
+@short Container for Vector of Vectors
+@author David Saxton
+*/
+class matrix
+{
+public:
+	matrix( const uint size );
+	~matrix();
+	
+	Vector & operator[]( const uint i ) { return *(m_mat[i]); }
+	const Vector & operator[]( const uint i ) const { return *(m_mat[i]); }
+	/**
+	 * Swaps the pointers to the given rows
+	 */
+	void swapRows( const uint a, const uint b );
+	
+private:
+	matrix( const matrix & );
+	matrix & operator= ( const matrix & );
+	
+	Vector **m_mat;
+	uint m_size;
+};
+
+#endif
diff --git a/src/electronics/simulation/voltagepoint.cpp b/src/electronics/simulation/voltagepoint.cpp
new file mode 100644
index 0000000..44d3b01
--- /dev/null
+++ b/src/electronics/simulation/voltagepoint.cpp
@@ -0,0 +1,73 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "matrix.h"
+#include "voltagepoint.h"
+#include "elementset.h"
+
+VoltagePoint::VoltagePoint( const double voltage )
+	: Element::Element()
+{
+	m_voltage = -voltage;
+	m_numCBranches = 1;
+	m_numCNodes = 1;
+}
+
+
+VoltagePoint::~VoltagePoint()
+{
+}
+
+
+void VoltagePoint::setVoltage( const double v )
+{
+	if ( -v == m_voltage ) return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	m_voltage = -v;
+	add_initial_dc();
+}
+
+
+void VoltagePoint::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+}
+
+
+void VoltagePoint::add_initial_dc()
+{
+	if (!b_status) return;
+	
+	A_b( 0, 0 ) = -1;
+	A_c( 0, 0 ) = -1;
+	
+	b_v( 0 ) = m_voltage;
+}
+
+
+void VoltagePoint::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[0] = p_cbranch[0]->i;
+}
+
+
+
+
+
diff --git a/src/electronics/simulation/voltagepoint.h b/src/electronics/simulation/voltagepoint.h
new file mode 100644
index 0000000..924b0af
--- /dev/null
+++ b/src/electronics/simulation/voltagepoint.h
@@ -0,0 +1,39 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VOLTAGEPOINT_H
+#define VOLTAGEPOINT_H
+
+#include "element.h"
+
+/**
+@short VoltagePoint
+@author David saxton
+*/
+class VoltagePoint : public Element
+{
+public:
+	VoltagePoint( const double voltage );
+	virtual ~VoltagePoint();
+	
+	virtual Type type() const { return Element_VoltagePoint; }
+	void setVoltage( const double voltage );
+	double voltage() { return m_voltage; }
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_voltage; // Conductance
+};
+
+#endif
diff --git a/src/electronics/simulation/voltagesignal.cpp b/src/electronics/simulation/voltagesignal.cpp
new file mode 100644
index 0000000..d7a5839
--- /dev/null
+++ b/src/electronics/simulation/voltagesignal.cpp
@@ -0,0 +1,79 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "matrix.h"
+#include "voltagesignal.h"
+#include "elementset.h"
+
+VoltageSignal::VoltageSignal( const double delta, const double voltage )
+	: Reactive::Reactive(delta)
+{
+	m_voltage = voltage;
+	m_numCNodes = 2;
+	m_numCBranches = 1;
+}
+
+
+VoltageSignal::~VoltageSignal()
+{
+}
+
+
+void VoltageSignal::setVoltage( const double v )
+{
+	m_voltage = v;
+}
+
+
+void VoltageSignal::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+	
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_constant, true );
+	}
+}
+
+
+void VoltageSignal::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_b( 0, 0 ) = -1;
+	A_c( 0, 0 ) = -1;
+	A_b( 1, 0 ) = 1;
+	A_c( 0, 1 ) = 1;
+}
+
+
+void VoltageSignal::time_step()
+{
+	if (!b_status) return;
+	b_v( 0 ) =  m_voltage*advance();
+}
+
+
+void VoltageSignal::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[1] = p_cbranch[0]->i;
+	m_cnodeI[0] = -m_cnodeI[1];
+}
+
+
diff --git a/src/electronics/simulation/voltagesignal.h b/src/electronics/simulation/voltagesignal.h
new file mode 100644
index 0000000..5e35e72
--- /dev/null
+++ b/src/electronics/simulation/voltagesignal.h
@@ -0,0 +1,41 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VOLTAGESIGNAL_H
+#define VOLTAGESIGNAL_H
+
+#include "reactive.h"
+#include "elementsignal.h"
+
+/**
+@short VoltageSignal
+@author David saxton
+*/
+class VoltageSignal : public Reactive, public ElementSignal
+{
+public:
+	VoltageSignal( const double delta, const double voltage );
+	virtual ~VoltageSignal();
+	
+	virtual Element::Type type() const { return Element_VoltageSignal; }
+	void setVoltage( const double voltage );
+	double voltage() { return m_voltage; }
+	virtual void time_step();
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_voltage; // Voltage
+};
+
+#endif
diff --git a/src/electronics/simulation/voltagesource.cpp b/src/electronics/simulation/voltagesource.cpp
new file mode 100644
index 0000000..b4fca64
--- /dev/null
+++ b/src/electronics/simulation/voltagesource.cpp
@@ -0,0 +1,77 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#include "matrix.h"
+#include "voltagesource.h"
+#include "elementset.h"
+
+VoltageSource::VoltageSource( const double voltage )
+	: Element::Element()
+{
+	m_v = voltage;
+	m_numCBranches = 1;
+	m_numCNodes = 2;
+}
+
+VoltageSource::~VoltageSource()
+{
+}
+
+void VoltageSource::setVoltage( const double v )
+{
+	if ( m_v == v ) return;
+	
+	if (p_eSet)
+		p_eSet->setCacheInvalidated();
+	
+	m_v = v;
+	add_initial_dc();
+}
+
+
+void VoltageSource::add_map()
+{
+	if (!b_status) return;
+	
+	if ( !p_cnode[0]->isGround )
+	{
+		p_A->setUse_b( p_cnode[0]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[0]->n(), Map::et_constant, true );
+	}
+	
+	if ( !p_cnode[1]->isGround )
+	{
+		p_A->setUse_b( p_cnode[1]->n(), p_cbranch[0]->n(), Map::et_constant, true );
+		p_A->setUse_c( p_cbranch[0]->n(), p_cnode[1]->n(), Map::et_constant, true );
+	}
+}
+
+
+void VoltageSource::add_initial_dc()
+{
+	if (!b_status)
+		return;
+	
+	A_b( 0, 0 ) = -1;
+	A_c( 0, 0 ) = -1;
+	A_b( 1, 0 ) = 1;
+	A_c( 0, 1 ) = 1;
+	
+	b_v( 0 ) = m_v;
+}
+
+void VoltageSource::updateCurrents()
+{
+	if (!b_status) return;
+	m_cnodeI[0] = p_cbranch[0]->i;
+	m_cnodeI[1] = -m_cnodeI[0];
+}
+
+
diff --git a/src/electronics/simulation/voltagesource.h b/src/electronics/simulation/voltagesource.h
new file mode 100644
index 0000000..9b27b0d
--- /dev/null
+++ b/src/electronics/simulation/voltagesource.h
@@ -0,0 +1,38 @@
+/***************************************************************************
+ *   Copyright (C) 2003-2004 by David Saxton                               *
+ *   david@bluehaze.org                                                    *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ ***************************************************************************/
+
+#ifndef VOLTAGESOURCE_H
+#define VOLTAGESOURCE_H
+
+#include "element.h"
+
+/**
+CNode n0 is the negative terminal, CNode n1 is the positive terminal
+@short Voltage Source
+*/
+class VoltageSource : public Element
+{
+public:
+	VoltageSource( const double voltage );
+	virtual ~VoltageSource();
+	
+	virtual Type type() const { return Element_VoltageSource; }
+	void setVoltage( const double v );
+	virtual void add_map();
+	
+protected:
+	virtual void updateCurrents();
+	virtual void add_initial_dc();
+	
+private:
+	double m_v; // Voltage
+};
+
+#endif