1 files changed, 97 insertions, 0 deletions
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;
+	}
+}
+