From 54bc96b04622c6c85bebd1e6af1fcd2f7eb748c9 Mon Sep 17 00:00:00 2001
From: Damien Zammit <damien@zamaudio.com>
Date: Mon, 14 Jan 2019 23:56:10 +1100
Subject: Introduce taylor expansion to reduce load by > 300x

---
 plugins/ZamTube/ZamTubePlugin.cpp | 13 +++++++------
 plugins/ZamTube/triode.cpp        | 33 ++++++++++++++++++++++-----------
 plugins/ZamTube/triode.h          |  8 --------
 plugins/ZamTube/wdfcircuits.h     |  3 +--
 4 files changed, 30 insertions(+), 27 deletions(-)

diff --git a/plugins/ZamTube/ZamTubePlugin.cpp b/plugins/ZamTube/ZamTubePlugin.cpp
index 1eb1328..719a2ae 100644
--- a/plugins/ZamTube/ZamTubePlugin.cpp
+++ b/plugins/ZamTube/ZamTubePlugin.cpp
@@ -251,7 +251,7 @@ void ZamTubePlugin::activate()
 	co[0] = 10e-9;
 	ro[0] = 1e+6;
 
-	/* Matt's preamp
+	/* Matt's preamp */
 	ci[1] = 100e-9;
 	rg[1] = 1.;
 	rk[1] = 820.;
@@ -260,9 +260,8 @@ void ZamTubePlugin::activate()
 	er[1] = 120e+3;
 	co[1] = 4.7e-9;
 	ro[1] = 470e+3;
-	*/
 
-	/* CLA's preamp (not good with all tonestacks)
+	/* CLA's preamp
 	ci[1] = 1.0e-7;
 	rg[1] = 5.6e+3;
 	rk[1] = 1.5e+3;
@@ -340,11 +339,13 @@ void ZamTubePlugin::run(const float** inputs, float** outputs, uint32_t frames)
 
 	float tubeout = 0.f;
 	
-	if (insaneold != insane) {
-		ckt.t.insane = insane;
+	if (insane != insaneold) {
+		ckt.reset_tubes();
+		ckt.warmup_tubes();
+		fRec0[3] = fRec0[2] = fRec0[1] = fRec0[0] = 0.f;
 		insaneold = insane;
 	}
-	int pre = 0;
+	int pre = insane < 0.5 ? 0 : 1;
 	float volumepot = tubedrive / 11. * 1e+6;
 	ckt.updateRValues(ci[pre], ck[pre], co[pre], e[pre], er[pre], rg[pre], volumepot, rk[pre], 136e+3, ro[pre], getSampleRate(), v);
 
diff --git a/plugins/ZamTube/triode.cpp b/plugins/ZamTube/triode.cpp
index 35e0b35..dfa2505 100644
--- a/plugins/ZamTube/triode.cpp
+++ b/plugins/ZamTube/triode.cpp
@@ -24,6 +24,18 @@ T Triode::compute(T a, T R, T Vg, T Vk) {
 }
 
 T Triode::getIa(T Vgk, T Vpk) {
+	static bool prepared = false;
+	static double coeff[3];
+
+	if (!prepared) {
+		const double L2 = log(2.0);
+		const double scale = 2e+9*pow(L2, kx-2.0)/(8.0*pow(kp, kx));
+		coeff[0] = 8.0*L2*L2*scale;
+		coeff[1] = kx*kp*L2*4.0*scale;
+		coeff[2] = (kp*kp*kx*kx + L2*kp*kp*kx - kp*kp*kx) * scale;
+		prepared = true;
+	}
+
 	if (Vpk < 0.0) {
 		//printf("Less than zero!\n");
 		Vpk = 0.0;
@@ -32,14 +44,16 @@ T Triode::getIa(T Vgk, T Vpk) {
 		Vgk = 0.0;
 	}
 
-	T e1 = Vpk*log1p(exp(kp*(1./mu+Vgk/pow(kvb+Vpk*Vpk, 0.5))))/kp;
-	if (e1 < 0) {
-		return 0.;
-	}
-	T test = 1e+6*pow(e1, kx) / kg1;
-	//printf("XXX Ip=%f lut=%f diff=%f  vgk=%f vpk=%f\n", test, test2, test2-test, Vgk, Vpk);
-	return test;
+	double A = 1./mu + Vgk / sqrt(kvb + Vpk*Vpk);
+	return (coeff[0] + coeff[1]*A + coeff[2]*A*A) / kg1;
 
+	/* exact solution (takes > 300x longer)
+		e1 = Vpk*log1p(exp(kp*(1./mu+Vgk/sqrt(kvb+Vpk*Vpk))))/kp;
+		if (e1 < 0) {
+			return 0.;
+		}
+		return 1e+6*pow(e1, kx) / kg1;
+	*/
 }
 
 T Triode::evaluateImplicitEquation(T Vak, T Vgk, T a, T R){
@@ -55,13 +69,10 @@ T Triode::iterateNewtonRaphson(T x, T dx, T Vgk, T a, T R){
 
 Triode::Triode()
 {
-	insane = true;
-
-	//12AX7 EHX-1
+	//12AX7 RSD-1
 	kvb = 300.;
 	mu = 103.2;
 	kx = 1.26;
 	kg1 = 446.0;
 	kp = 3.4;
 }
-
diff --git a/plugins/ZamTube/triode.h b/plugins/ZamTube/triode.h
index 5a4e3d3..ae552fe 100644
--- a/plugins/ZamTube/triode.h
+++ b/plugins/ZamTube/triode.h
@@ -30,13 +30,6 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 typedef double T;
 
-static inline float
-sanitize_denormal(float v) {
-        if(!std::isnormal(v) || !std::isfinite(v))
-                return 0.f;
-        return v;
-}
-
 class Triode {
 public:
 	T mu, kp, kvb, kg1, kx;
@@ -46,7 +39,6 @@ public:
 	T getIa(T Vgk, T Vak);
 	T evaluateImplicitEquation(T Vak, T Vgk, T a, T R);
 	T iterateNewtonRaphson(T x, T dx, T Vgk, T a, T R);
-	const float *lut_ip;
 	bool insane;
 };
 
diff --git a/plugins/ZamTube/wdfcircuits.h b/plugins/ZamTube/wdfcircuits.h
index 02dfee8..37203a1 100644
--- a/plugins/ZamTube/wdfcircuits.h
+++ b/plugins/ZamTube/wdfcircuits.h
@@ -12,6 +12,7 @@ public:
 	TubeStageCircuit() {
 		on = false;
 		reset_tubes();
+		warmup_tubes();
 	}
 
 	void warmup_tubes(void) {
@@ -29,10 +30,8 @@ public:
 		Coa = 0.0;
 		Vg = 0.0;
 		Vk = 0.0;
-		warmup_tubes();
 	}
 
-
 	void updateRValues(Real C_Ci, Real C_Ck, Real C_Co, Real E_E500, Real R_E500, Real R_Rg, Real R_Ri, Real R_Rk, Real R_Vi, Real R_Ro, Real sampleRate, Triode& tube) {
 		t = tube;
 		Real ViR = R_Vi;
-- 
cgit v1.2.3