Simpler circuit, amplification is good

5 files changed, 175 insertions, 730 deletions

diff --git a/plugins/ZamTube/ZamTubePlugin.cpp b/plugins/ZamTube/ZamTubePlugin.cpp
index ca149f2..1fa17ef 100644
--- a/plugins/ZamTube/ZamTubePlugin.cpp
+++ b/plugins/ZamTube/ZamTubePlugin.cpp
@@ -203,15 +203,16 @@ void ZamTubePlugin::activate()
 	T Fs = getSampleRate();
 	
 	// Passive components
-	ci = 0.0000001	;	//100nF
-	ck = 0.00001;		//10uF
-	co = 0.00000001;	//10nF
-	ro = 1000000.0;		//1Mohm
-	rp = 100000.0;		//100kohm
-	rg = 20000.0;		//20kohm
-	ri = 1000000.0;		//1Mohm
-	rk = 1000.0;		//1kohm
-	e = 250.0;		//250V
+	ci = 1.0e-7;
+	ri = 1.0e+6;
+	rg = 2.0e+4;
+	rk = 1.0e+3;
+	ck = 1.0e-5;
+	e = 500.0;
+	rp = 1.0e+5;
+	er = rp;
+	co = 1.0e-8;
+	ro = 1.0e+6;
 
 	// 12AX7 triode model RSD-1
 	v.g = 2.242e-3;
@@ -233,7 +234,7 @@ void ZamTubePlugin::activate()
 	v.cg2 = 11.99;
 	v.ig02 = 3.917e-8;
 
-	ckt.updateRValues(ci, ck, co, e, rp, rg, ri, rk, ro, 1000.0, Fs, v);
+	ckt.updateRValues(ci, ck, co, e, er, rg, ri, rk, 1., ro, Fs, v);
 	ckt.warmup_tubes();
 
         fSamplingFreq = Fs;
@@ -949,15 +950,18 @@ void ZamTubePlugin::run(const float** inputs, float** outputs, uint32_t frames)
 
 		// protect against overflowing circuit
 		in = fabs(in) < DANGER ? in : 0.f; 
-
-		double ViE = in*from_dB(tubedrive - 10.);
-		tubeout = e * ckt.advanc(ViE) * from_dB(40. - tubedrive)*from_dB(-3.);
-		if (!ckt.on) {
-			tubeout = 0.0;
+	
+		double pregain = from_dB(tubedrive*3.6364); // 0 to +40dB
+		double ViE = in*pregain;
+		double postgain = from_dB(tubedrive/2.) * from_dB(mastergain);
+		tubeout = ckt.advanc(ViE) * postgain / (e * 0.5);
+		/*
+		if (insane) {
+			tubeout = in*from_dB(-50.)*postgain;
 		} else {
-			tubeout = sanitize_denormal(tubeout);
+			tubeout = ckt.advanc(ViE) * postgain - in*from_dB(-50.)*postgain;
 		}
-
+		*/
 		//Tone Stack sim
 		fRec0[0] = sanitize_denormal((float)tubeout - (fSlow16 * (((fSlow14 * fRec0[2]) + (fSlow13 * fRec0[1])) + (fSlow11 * fRec0[3]))));
 		fRec1[0] = sanitize_denormal((float)tubeout - (fSlow47 * (((fSlow45 * fRec1[2]) + (fSlow44 * fRec1[1])) + (fSlow42 * fRec1[3]))));
@@ -984,8 +988,9 @@ void ZamTubePlugin::run(const float** inputs, float** outputs, uint32_t frames)
 		fRec22[0] = sanitize_denormal((float)tubeout - (fSlow602 * (((fSlow600 * fRec22[2]) + (fSlow599 * fRec22[1])) + (fSlow597 * fRec22[3]))));
 		fRec23[0] = sanitize_denormal((float)tubeout - (fSlow626 * (((fSlow624 * fRec23[2]) + (fSlow623 * fRec23[1])) + (fSlow621 * fRec23[3]))));
 		fRec24[0] = sanitize_denormal((float)tubeout - (fSlow652 * (((fSlow650 * fRec24[2]) + (fSlow649 * fRec24[1])) + (fSlow647 * fRec24[3]))));
-		outputs[0][i] = (FAUSTFLOAT)((fSlow664 * ((fSlow663 * fRec24[0]) + ((fSlow662 * fRec24[1]) + ((fSlow660 * fRec24[3]) + (fSlow658 * fRec24[2]))))) + ((fSlow638 * ((fSlow637 * fRec23[0]) + ((fSlow636 * fRec23[1]) + ((fSlow634 * fRec23[3]) + (fSlow632 * fRec23[2]))))) + ((fSlow614 * ((fSlow613 * fRec22[0]) + ((fSlow612 * fRec22[1]) + ((fSlow610 * fRec22[3]) + (fSlow608 * fRec22[2]))))) + ((fSlow588 * ((fSlow587 * fRec21[0]) + ((fSlow586 * fRec21[1]) + ((fSlow584 * fRec21[3]) + (fSlow582 * fRec21[2]))))) + ((fSlow561 * ((fSlow560 * fRec20[0]) + ((fSlow559 * fRec20[1]) + ((fSlow558 * fRec20[3]) + (fSlow556 * fRec20[2]))))) + ((fSlow540 * ((fSlow539 * fRec19[0]) + ((fSlow538 * fRec19[1]) + ((fSlow537 * fRec19[3]) + (fSlow535 * fRec19[2]))))) + ((fSlow519 * ((fSlow518 * fRec18[0]) + ((fSlow517 * fRec18[1]) + ((fSlow515 * fRec18[3]) + (fSlow513 * fRec18[2]))))) + ((fSlow493 * ((fSlow492 * fRec17[0]) + ((fSlow491 * fRec17[1]) + ((fSlow489 * fRec17[3]) + (fSlow487 * fRec17[2]))))) + ((fSlow463 * ((fSlow462 * fRec16[0]) + ((fSlow461 * fRec16[1]) + ((fSlow459 * fRec16[3]) + (fSlow457 * fRec16[2]))))) + ((fSlow435 * ((fSlow434 * fRec15[0]) + ((fSlow433 * fRec15[1]) + ((fSlow431 * fRec15[3]) + (fSlow429 * fRec15[2]))))) + ((fSlow409 * ((fSlow408 * fRec14[0]) + ((fSlow407 * fRec14[1]) + ((fSlow405 * fRec14[3]) + (fSlow403 * fRec14[2]))))) + ((fSlow385 * ((fSlow384 * fRec13[0]) + ((fSlow383 * fRec13[1]) + ((fSlow381 * fRec13[3]) + (fSlow379 * fRec13[2]))))) + ((fSlow358 * ((fSlow357 * fRec12[0]) + ((fSlow356 * fRec12[1]) + ((fSlow354 * fRec12[3]) + (fSlow352 * fRec12[2]))))) + ((fSlow332 * ((fSlow331 * fRec11[0]) + ((fSlow330 * fRec11[1]) + ((fSlow328 * fRec11[3]) + (fSlow326 * fRec11[2]))))) + ((fSlow305 * ((fSlow304 * fRec10[0]) + ((fSlow303 * fRec10[1]) + ((fSlow301 * fRec10[3]) + (fSlow299 * fRec10[2]))))) + ((fSlow277 * ((fSlow276 * fRec9[0]) + ((fSlow275 * fRec9[1]) + ((fSlow273 * fRec9[3]) + (fSlow271 * fRec9[2]))))) + ((fSlow250 * ((fSlow249 * fRec8[0]) + ((fSlow248 * fRec8[1]) + ((fSlow246 * fRec8[3]) + (fSlow244 * fRec8[2]))))) + ((fSlow224 * ((fSlow223 * fRec7[0]) + ((fSlow222 * fRec7[1]) + ((fSlow220 * fRec7[3]) + (fSlow218 * fRec7[2]))))) + ((fSlow198 * ((fSlow197 * fRec6[0]) + ((fSlow196 * fRec6[1]) + ((fSlow194 * fRec6[3]) + (fSlow192 * fRec6[2]))))) + ((fSlow171 * ((fSlow170 * fRec5[0]) + ((fSlow169 * fRec5[1]) + ((fSlow167 * fRec5[3]) + (fSlow165 * fRec5[2]))))) + ((fSlow140 * ((fSlow139 * fRec4[0]) + ((fSlow138 * fRec4[1]) + ((fSlow136 * fRec4[3]) + (fSlow134 * fRec4[2]))))) + ((fSlow114 * ((fSlow113 * fRec3[0]) + ((fSlow112 * fRec3[1]) + ((fSlow110 * fRec3[3]) + (fSlow108 * fRec3[2]))))) + ((fSlow88 * ((fSlow87 * fRec2[0]) + ((fSlow86 * fRec2[1]) + ((fSlow84 * fRec2[3]) + (fSlow82 * fRec2[2]))))) + ((fSlow60 * ((fSlow59 * fRec1[0]) + ((fSlow58 * fRec1[1]) + ((fSlow56 * fRec1[3]) + (fSlow54 * fRec1[2]))))) + (fSlow32 * ((fSlow30 * fRec0[0]) + ((fSlow29 * fRec0[1]) + ((fSlow27 * fRec0[3]) + (fSlow25 * fRec0[2])))))))))))))))))))))))))))))* from_dB(mastergain);
-		outputs[0][i] = sanitize_denormal(outputs[0][i]);
+		double tonestackout = (FAUSTFLOAT)((fSlow664 * ((fSlow663 * fRec24[0]) + ((fSlow662 * fRec24[1]) + ((fSlow660 * fRec24[3]) + (fSlow658 * fRec24[2]))))) + ((fSlow638 * ((fSlow637 * fRec23[0]) + ((fSlow636 * fRec23[1]) + ((fSlow634 * fRec23[3]) + (fSlow632 * fRec23[2]))))) + ((fSlow614 * ((fSlow613 * fRec22[0]) + ((fSlow612 * fRec22[1]) + ((fSlow610 * fRec22[3]) + (fSlow608 * fRec22[2]))))) + ((fSlow588 * ((fSlow587 * fRec21[0]) + ((fSlow586 * fRec21[1]) + ((fSlow584 * fRec21[3]) + (fSlow582 * fRec21[2]))))) + ((fSlow561 * ((fSlow560 * fRec20[0]) + ((fSlow559 * fRec20[1]) + ((fSlow558 * fRec20[3]) + (fSlow556 * fRec20[2]))))) + ((fSlow540 * ((fSlow539 * fRec19[0]) + ((fSlow538 * fRec19[1]) + ((fSlow537 * fRec19[3]) + (fSlow535 * fRec19[2]))))) + ((fSlow519 * ((fSlow518 * fRec18[0]) + ((fSlow517 * fRec18[1]) + ((fSlow515 * fRec18[3]) + (fSlow513 * fRec18[2]))))) + ((fSlow493 * ((fSlow492 * fRec17[0]) + ((fSlow491 * fRec17[1]) + ((fSlow489 * fRec17[3]) + (fSlow487 * fRec17[2]))))) + ((fSlow463 * ((fSlow462 * fRec16[0]) + ((fSlow461 * fRec16[1]) + ((fSlow459 * fRec16[3]) + (fSlow457 * fRec16[2]))))) + ((fSlow435 * ((fSlow434 * fRec15[0]) + ((fSlow433 * fRec15[1]) + ((fSlow431 * fRec15[3]) + (fSlow429 * fRec15[2]))))) + ((fSlow409 * ((fSlow408 * fRec14[0]) + ((fSlow407 * fRec14[1]) + ((fSlow405 * fRec14[3]) + (fSlow403 * fRec14[2]))))) + ((fSlow385 * ((fSlow384 * fRec13[0]) + ((fSlow383 * fRec13[1]) + ((fSlow381 * fRec13[3]) + (fSlow379 * fRec13[2]))))) + ((fSlow358 * ((fSlow357 * fRec12[0]) + ((fSlow356 * fRec12[1]) + ((fSlow354 * fRec12[3]) + (fSlow352 * fRec12[2]))))) + ((fSlow332 * ((fSlow331 * fRec11[0]) + ((fSlow330 * fRec11[1]) + ((fSlow328 * fRec11[3]) + (fSlow326 * fRec11[2]))))) + ((fSlow305 * ((fSlow304 * fRec10[0]) + ((fSlow303 * fRec10[1]) + ((fSlow301 * fRec10[3]) + (fSlow299 * fRec10[2]))))) + ((fSlow277 * ((fSlow276 * fRec9[0]) + ((fSlow275 * fRec9[1]) + ((fSlow273 * fRec9[3]) + (fSlow271 * fRec9[2]))))) + ((fSlow250 * ((fSlow249 * fRec8[0]) + ((fSlow248 * fRec8[1]) + ((fSlow246 * fRec8[3]) + (fSlow244 * fRec8[2]))))) + ((fSlow224 * ((fSlow223 * fRec7[0]) + ((fSlow222 * fRec7[1]) + ((fSlow220 * fRec7[3]) + (fSlow218 * fRec7[2]))))) + ((fSlow198 * ((fSlow197 * fRec6[0]) + ((fSlow196 * fRec6[1]) + ((fSlow194 * fRec6[3]) + (fSlow192 * fRec6[2]))))) + ((fSlow171 * ((fSlow170 * fRec5[0]) + ((fSlow169 * fRec5[1]) + ((fSlow167 * fRec5[3]) + (fSlow165 * fRec5[2]))))) + ((fSlow140 * ((fSlow139 * fRec4[0]) + ((fSlow138 * fRec4[1]) + ((fSlow136 * fRec4[3]) + (fSlow134 * fRec4[2]))))) + ((fSlow114 * ((fSlow113 * fRec3[0]) + ((fSlow112 * fRec3[1]) + ((fSlow110 * fRec3[3]) + (fSlow108 * fRec3[2]))))) + ((fSlow88 * ((fSlow87 * fRec2[0]) + ((fSlow86 * fRec2[1]) + ((fSlow84 * fRec2[3]) + (fSlow82 * fRec2[2]))))) + ((fSlow60 * ((fSlow59 * fRec1[0]) + ((fSlow58 * fRec1[1]) + ((fSlow56 * fRec1[3]) + (fSlow54 * fRec1[2]))))) + (fSlow32 * ((fSlow30 * fRec0[0]) + ((fSlow29 * fRec0[1]) + ((fSlow27 * fRec0[3]) + (fSlow25 * fRec0[2])))))))))))))))))))))))))))));
+		
+		outputs[0][i] = sanitize_denormal(tonestackout);
 
 		// post processing
 		for (int i=3; i>0; i--) fRec24[i] = sanitize_denormal(fRec24[i-1]);
diff --git a/plugins/ZamTube/ZamTubePlugin.hpp b/plugins/ZamTube/ZamTubePlugin.hpp
index b6956c4..f583da4 100644
--- a/plugins/ZamTube/ZamTubePlugin.hpp
+++ b/plugins/ZamTube/ZamTubePlugin.hpp
@@ -41,13 +41,14 @@ class ZamTubePlugin : public Plugin
 public:
 	Triode v;
 	TubeStageCircuit ckt;
-	T e;
 	T ci;
 	T ck;
 	T co;
-	T ro;
+	T e;
+	T er;
 	T rp;
 	T rg;
+	T ro;
 	T ri;
 	T rk;
 
diff --git a/plugins/ZamTube/triode.cpp b/plugins/ZamTube/triode.cpp
index 045be35..3b70782 100644
--- a/plugins/ZamTube/triode.cpp
+++ b/plugins/ZamTube/triode.cpp
@@ -3,205 +3,64 @@
 #include <cmath>
 #include "triode.h"
 using std::abs;
-
-T Triode::getC(void)
-{
-	return Kb;
-}
-
-T Triode::getP(void)
-{
-	return Pb;
-}
-
-T Triode::getG(void)
-{
-	return Gb;
-}
-
-void Triode::compute(T Pbb, T Gbb, T Kbb)
-{
-//	T Kb_o = Kb;
-//	T Gb_o = Gb;
-//	T Pb_o = Pb;
-
-//	Kb = (2.0*vk-Kbb);
-//	Gb = (2.0*vg-Gbb);
-//	Pb = (2.0*vp-Pbb);
-
-	Kb = Kbb;
-	Gb = Gbb;
-	Pb = Pbb;
-
-	//Step 3: compute wave reflections inside the triode
-	T vg0, vg1, vp0, vp1;
-
-	vg0 = -10.0;
-	vg1 = 10.0;
-	vg = sanitize_denormal(zeroffg(vg0, vg1, TOLERANCE));
-	//v.vg = v.secantfg(&vg0,&vg1);
-
-	vp0 = e;
-	vp1 = 0.0;
-	if (insane) {
-		vp = sanitize_denormal(zeroffp_insane(vp0,vp1,TOLERANCE));
-	} else {
-		vp = sanitize_denormal(zeroffp(vp0,vp1,TOLERANCE));
+#define DUMP(x) x
+
+T Triode::compute(T a, T R, T Vgate, T Vk) {
+	T VakGuess = 100.;
+	
+	T Vgk = Vgate - Vk;
+
+	T Vak = VakGuess; // initial guess
+	int iteration = 0;
+	T err = 1e6;
+	T Iak = 0.0;
+	while (fabs(err)/fabs(Vak) > 1e-9){
+		VakGuess = iterateNewtonRaphson(Vak, 1e-6, Vgk, a, R);
+		err = Vak - VakGuess;
+		Vak = VakGuess;
+
+		if (iteration > 100){
+			printf("Convergence failure!");
+			break;
+		}
+		++iteration;
 	}
-	//v.vp = v.secantfp(&vp0,&vp1);
+	T b = Vak - R*Iak;
 
-	vk = sanitize_denormal(ffk());
-
-	Kb = (2.0*vk-Kb);
-	Gb = (2.0*vg-Gb);
-	Pb = (2.0*vp-Pb);
-
-	//printf("vg = %f  vp = %f  vk = %f  Gb = %f  Kb = %f  Pb = %f\n", vg, vp, vk, Gb, Kb, Pb);
+	//printf("Vgate=%f Vk=%f  Vgk=%f b=%f\n", Vgate, Vk, Vgk, b);
+	return b;
 }
 
-T Triode::ffg(T VG) {
-        return (Gb-Gr*(gg*pow(log(1.0+exp(cg*VG))/cg,e)+ig0)-VG);
-}
-
-T Triode::fgdash(T VG) {
-        T a1 = exp(cg*VG);
-        T b1 = -e*pow(log(a1+1.0)/cg,e-1.0);
-        T c1 = a1/(a1+1.0)*gg*Gr;
-        return (b1*c1);
-}
-
-T Triode::ffp(T VP) {
-	static bool prepared = false;
-	if(!prepared) {
-		//go go series expansion
-		const double L2 = log(2.0);
-
-		const double scale = pow(L2,gamma-2)/(8.0*pow(c,gamma));
-		ffp_raw[0] = 8.0*L2*L2*scale;
-		ffp_raw[1] = gamma*c*L2*4*scale;
-		ffp_raw[2] = (c*c*gamma*scale)*(gamma+L2-1.);
-		prepared = true;
+T Triode::getIa(T Vgk, T Vak) {
+	if (Vak < 0.0) {
+		printf("Less than zero!\n");
+		Vak = 0.0;
 	}
-
-	return ffp_coeff[0]+ffp_coeff[1]*VP+ffp_coeff[2]*VP*VP;
-}
-
-void Triode::prepare(void)
-{
-    const double c0 = ffp_raw[0],
-                 c1 = ffp_raw[1],
-                 c2 = ffp_raw[2];
-    ffp_coeff[0] = Pb + Pr*(Gb/Gr + c2*vg*vg - vg/Gr + c1*vg + c0);
-    ffp_coeff[1] = (2.*c2*vg + c1)*(Pr/mu) - 1.;
-    ffp_coeff[2] = c2*Pr/(mu*mu);
-}
-
-T Triode::ffp_insane(T VP) {
-	return (Pb+Pr*((g*pow(log(1.0+exp(c*(VP/mu+vg)))/c,gamma))+(Gb-vg)/Gr)-VP);
-}
-
-T Triode::fpdash(T VP) {
-        T a1 = exp(c*(vg+VP/mu));
-        T b1 = a1/(mu*(a1+1.0));
-        T c1 = g*gamma*Pr*pow(log(a1+1.0)/c,gamma-1.0);
-        return (c1*b1);
-}
-
-T Triode::ffk() {
-        return (Kb - Kr*(g*pow(log(1.0+exp(c*(vp/mu+vg)))/c,gamma)));
-}
-/*
-T Triode::secantfg(T *i1, T *i2) {
-        T vgn = 0.0;
-        T init = *i1;
-        for (int i = 0; i<ITER; ++i) {
-                vgn = *i1 - fg(*i1)*(*i1-*i2)/(fg(*i1)-fg(*i2));
-                *i2 = *i1;
-                *i1 = vgn;
-                if ((fabs(fg(vgn))) < EPSILON) break;
-        }
-        if ((fabs(fg(vgn)) >= EPSILON)) {
-                DUMP(fprintf(stderr,"Vg did not converge\n"));
-		return init;
+	if (Vgk > 0.0) {
+		Vgk = 0.0;
 	}
-        return vgn; 
-}               
-        
-T Triode::newtonfg(T *i1) {
-        T init = *i1;
-        if (fabs(fg(*i1)) < EPSILON || fgdash(*i1)==0.0) return *i1;
-        T vgn = 0.0;
-        for (int i = 0; i<ITER; ++i) {
-                vgn = *i1 - fg(*i1)/fgdash(*i1);
-                *i1 = vgn;
-                if (fabs(fg(vgn)) < EPSILON) break;
-        } 
-        if ((fabs(fg(vgn)) >= EPSILON)) {
-//                vgn = init;
-                DUMP(fprintf(stderr,"Vg did not converge\n"));
-        }       
-        return vgn;
-}
-
-T Triode::newtonfp(T *i1) {
-        T init = *i1;
-        if (fabs(fp(*i1)) < EPSILON || fpdash(*i1)==0.0) return *i1;
-        T vpn = 0.0;
-        for (int i = 0; i<ITER; ++i) {
-                vpn = *i1 - fp(*i1)/fpdash(*i1);
-                *i1 = vpn;
-                if (fabs(fp(vpn)) < EPSILON) break;
-        }
-        if ((fabs(fp(vpn)) >= EPSILON)) {
-//                vpn = init;
-                DUMP(fprintf(stderr,"Vp did not converge\n"));
-        }
-        return vpn;
+	T mu = 100.;
+	T kx = 1.4;
+	T kg1 = 3.981e-8;
+	T kp = 600.;
+	T kvb = 300.;
+
+	T e1 = Vak/kp*log(1.+exp(kp*(1./mu+Vgk/pow(kvb+Vak*Vak, 0.5))));
+	if (e1 < 0) {
+		return 0.;
+	}
+	return (pow(e1, kx) / kg1);
 }
 
-T Triode::secantfp(T *i1, T *i2) {
-        T vpn = 0.0;
-        for (int i = 0; i<ITER; ++i) {
-                vpn = *i1 - fp(*i1)*(*i1-*i2)/(fp(*i1)-fp(*i2));
-                *i2 = *i1;
-                *i1 = vpn;
-                if ((fabs(fp(vpn))) < EPSILON) break;
-        }
-
-        if ((fabs(fp(vpn)) >= EPSILON))
-                DUMP(fprintf(stderr,"Vp did not converge\n"));
-        return vpn;
+T Triode::evaluateImplicitEquation(T Vak, T Vgk, T a, T R){
+	T Iak = getIa(Vgk, Vak);
+	return Vak + R*Iak - a;
 }
-*/
 
-//****************************************************************************80
-//	Purpose:
-//
-//		Brent's method root finder.
-//
-//	Licensing:
-//
-//		This code below is distributed under the GNU LGPL license.
-//
-//	Author:
-//
-//		Original FORTRAN77 version by Richard Brent.
-//		C++ version by John Burkardt.
-//		Adapted for zamvalve by Damien Zammit.
-
-T Triode::r8_abs ( T x )
-{
-	T value;
-
-	if ( 0.0 <= x )
-	{
-		value = x;
-	}
-	else
-	{
-		value = - x;
-	}
-	return value;
+T Triode::iterateNewtonRaphson(T x, T dx, T Vgk, T a, T R){
+	T F = evaluateImplicitEquation(x, Vgk, a, R);
+	T xNew = x - dx*F/(evaluateImplicitEquation(x + dx, Vgk, a, R) - F);
+	return xNew;
 }
 
 Triode::Triode()
@@ -210,424 +69,4 @@ Triode::Triode()
 	vk = 0.0;
 	vp = 0.0;
 	insane = false;
-
-	T r = 1.0;
-
-	while ( 1.0 < ( T ) ( 1.0 + r )	)
-	{
-		r = r / 2.0;
-	}
-
-	r *= 2.0;
-	r8_epsilon = r;
-}
-
-T Triode::r8_max ( T x, T y )
-{
-	T value;
-
-	if ( y < x )
-	{
-		value = x;
-	}
-	else
-	{
-		value = y;
-	}
-	return value;
-}
-
-T Triode::r8_sign ( T x )
-{
-	T value;
-
-	if ( x < 0.0 )
-	{
-		value = -1.0;
-	}
-	else
-	{
-		value = 1.0;
-	}
-	return value;
-}
-
-T Triode::zeroffp_insane ( T a, T b, T t )
-{
-	T c;
-	T d;
-	T e;
-	T fa;
-	T fb;
-	T fc;
-	T m;
-	T macheps;
-	T p;
-	T q;
-	T r;
-	T s;
-	T sa;
-	T sb;
-	T tol;
-//
-//	Make local copies of A and B.
-//
-	sa = a;
-	sb = b;
-	fa = ffp_insane ( sa );
-	fb = ffp_insane ( sb );
-
-	c = sa;
-	fc = fa;
-	e = sb - sa;
-	d = e;
-
-	macheps = r8_epsilon;
-
-	for ( ; ; )
-	{
-		if ( abs ( fc ) < abs ( fb ) )
-		{
-			sa = sb;
-			sb = c;
-			c = sa;
-			fa = fb;
-			fb = fc;
-			fc = fa;
-		}
-
-		tol = 2.0 * macheps * abs ( sb ) + t;
-		m = 0.5 * ( c - sb );
-
-		if ( abs ( m ) <= tol || fb == 0.0 )
-		{
-			break;
-		}
-
-		if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )
-		{
-			e = m;
-			d = e;
-		}
-		else
-		{
-			s = fb / fa;
-
-			if ( sa == c )
-			{
-				p = 2.0 * m * s;
-				q = 1.0 - s;
-			}
-			else
-			{
-				q = fa / fc;
-				r = fb / fc;
-				p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
-				q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
-			}
-
-			if ( 0.0 < p )
-			{
-				q = - q;
-			}
-			else
-			{
-				p = - p;
-			}
-
-			s = e;
-			e = d;
-
-			if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) &&
-				p < abs ( 0.5 * s * q ) )
-			{
-				d = p / q;
-			}
-			else
-			{
-				e = m;
-				d = e;
-			}
-		}
-		sa = sb;
-		fa = fb;
-
-		if ( tol < abs ( d ) )
-		{
-			sb = sb + d;
-		}
-		else if ( 0.0 < m )
-		{
-			sb = sb + tol;
-		}
-		else
-		{
-			sb = sb - tol;
-		}
-
-		fb = ffp_insane ( sb );
-
-		if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
-		{
-			c = sa;
-			fc = fa;
-			e = sb - sa;
-			d = e;
-		}
-	}
-	return sb;
-}
-
-T Triode::zeroffp ( T a, T b, T t )
-{
-	T c;
-	T d;
-	T e;
-	T fa;
-	T fb;
-	T fc;
-	T m;
-	T macheps;
-	T p;
-	T q;
-	T r;
-	T s;
-	T sa;
-	T sb;
-	T tol;
-
-	prepare();
-//
-//	Make local copies of A and B.
-//
-	sa = a;
-	sb = b;
-	fa = ffp ( sa );
-	fb = ffp ( sb );
-
-	c = sa;
-	fc = fa;
-	e = sb - sa;
-	d = e;
-
-	macheps = r8_epsilon;
-
-	for ( ; ; )
-	{
-		if ( abs ( fc ) < abs ( fb ) )
-		{
-			sa = sb;
-			sb = c;
-			c = sa;
-			fa = fb;
-			fb = fc;
-			fc = fa;
-		}
-
-		tol = 2.0 * macheps * abs ( sb ) + t;
-		m = 0.5 * ( c - sb );
-
-		if ( abs ( m ) <= tol || fb == 0.0 )
-		{
-			break;
-		}
-
-		if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )
-		{
-			e = m;
-			d = e;
-		}
-		else
-		{
-			s = fb / fa;
-
-			if ( sa == c )
-			{
-				p = 2.0 * m * s;
-				q = 1.0 - s;
-			}
-			else
-			{
-				q = fa / fc;
-				r = fb / fc;
-				p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
-				q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
-			}
-
-			if ( 0.0 < p )
-			{
-				q = - q;
-			}
-			else
-			{
-				p = - p;
-			}
-
-			s = e;
-			e = d;
-
-			if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) &&
-				p < abs ( 0.5 * s * q ) )
-			{
-				d = p / q;
-			}
-			else
-			{
-				e = m;
-				d = e;
-			}
-		}
-		sa = sb;
-		fa = fb;
-
-		if ( tol < abs ( d ) )
-		{
-			sb = sb + d;
-		}
-		else if ( 0.0 < m )
-		{
-			sb = sb + tol;
-		}
-		else
-		{
-			sb = sb - tol;
-		}
-
-		fb = ffp ( sb );
-
-		if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
-		{
-			c = sa;
-			fc = fa;
-			e = sb - sa;
-			d = e;
-		}
-	}
-	return sb;
-}
-
-T Triode::zeroffg ( T a, T b, T t )
-{
-	T c;
-	T d;
-	T e;
-	T fa;
-	T fb;
-	T fc;
-	T m;
-	T macheps;
-	T p;
-	T q;
-	T r;
-	T s;
-	T sa;
-	T sb;
-	T tol;
-//
-//	Make local copies of A and B.
-//
-	sa = a;
-	sb = b;
-	fa = ffg ( sa );
-	fb = ffg ( sb );
-
-	c = sa;
-	fc = fa;
-	e = sb - sa;
-	d = e;
-
-	macheps = r8_epsilon;
-
-	for ( ; ; )
-	{
-		if ( abs ( fc ) < abs ( fb ) )
-		{
-			sa = sb;
-			sb = c;
-			c = sa;
-			fa = fb;
-			fb = fc;
-			fc = fa;
-		}
-
-		tol = 2.0 * macheps * abs ( sb ) + t;
-		m = 0.5 * ( c - sb );
-
-		if ( abs ( m ) <= tol || fb == 0.0 )
-		{
-			break;
-		}
-
-		if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )
-		{
-			e = m;
-			d = e;
-		}
-		else
-		{
-			s = fb / fa;
-
-			if ( sa == c )
-			{
-				p = 2.0 * m * s;
-				q = 1.0 - s;
-			}
-			else
-			{
-				q = fa / fc;
-				r = fb / fc;
-				p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );
-				q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );
-			}
-
-			if ( 0.0 < p )
-			{
-				q = - q;
-			}
-			else
-			{
-				p = - p;
-			}
-
-			s = e;
-			e = d;
-
-			if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) &&
-				p < abs ( 0.5 * s * q ) )
-			{
-				d = p / q;
-			}
-			else
-			{
-				e = m;
-				d = e;
-			}
-		}
-		sa = sb;
-		fa = fb;
-
-		if ( tol < abs ( d ) )
-		{
-			sb = sb + d;
-		}
-		else if ( 0.0 < m )
-		{
-			sb = sb + tol;
-		}
-		else
-		{
-			sb = sb - tol;
-		}
-
-		fb = ffg ( sb );
-
-		if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )
-		{
-			c = sa;
-			fc = fa;
-			e = sb - sa;
-			d = e;
-		}
-	}
-	return sb;
 }
diff --git a/plugins/ZamTube/triode.h b/plugins/ZamTube/triode.h
index a3b50ba..36e1036 100644
--- a/plugins/ZamTube/triode.h
+++ b/plugins/ZamTube/triode.h
@@ -53,38 +53,13 @@ public:
 	T g1, mu1, gamma1, c1, gg1, e1, cg1, ig01;
 	T g2, mu2, gamma2, c2, gg2, e2, cg2, ig02;
 
-	T ffg(T VG);
-	T fgdash(T VG);
-	T ffp(T VP);
-	T ffp_insane(T VP);
-	T fpdash(T VP);
-	T ffk();
-	T secantfg(T *i1, T *i2);
-	T newtonfg(T *i1);
-	T secantfp(T *i1, T *i2);
-	T newtonfp(T *i1);
 	bool insane;
 	
 	Triode();
-	void compute(T Kbb, T Gbb, T Pbb);
-	void prepare(void);
-	T getC(void);
-	T getG(void);
-	T getP(void);
-
-	//Brent's method
-	T r8_abs ( T x );
-	T r8_epsilon;
-	T r8_max ( T x, T y );
-	T r8_sign ( T x );
-	T zeroffp ( T a, T b, T t );
-	T zeroffp_insane ( T a, T b, T t );
-	T zeroffg ( T a, T b, T t );
-
-private:
-	//Taylor series coefficients for fast calculations
-	double ffp_raw[3];
-	double ffp_coeff[3];
+	T compute(T Kbb, T Gbb, T Pbb, T R);
+	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);
 };
 
 #endif
diff --git a/plugins/ZamTube/wdfcircuits.h b/plugins/ZamTube/wdfcircuits.h
index f8d8a49..2f774f5 100644
--- a/plugins/ZamTube/wdfcircuits.h
+++ b/plugins/ZamTube/wdfcircuits.h
@@ -13,28 +13,34 @@ public:
 		Cia = 0.0;
 		Cka = 0.0;
 		Coa = 0.0;
+		Vk = 0.0;
+		Vg = 0.0;
 		on = false;
 	}
-	TubeStageCircuit(Real C_Ci, Real C_Ck, Real C_Co, Real E_E250, Real R_E250, Real R_Rg, Real R_Ri, Real R_Rk, Real R_Ro, Real R_Vi, Real sampleRate, Triode& tube) {
-		Cia = 0.0;
-		Cka = 0.0;
-		Coa = 0.0;
+
+	void warmup_tubes(void) {
+		int i;
 		on = false;
-		updateRValues(C_Ci, C_Ck, C_Co, E_E250, R_E250, R_Rg, R_Ri, R_Rk, R_Ro, R_Vi, sampleRate, tube);
+		for (i = 0; i < 8000; i++) {
+			advanc(0.0);
+		}
+		on = true;
 	}
 
-	void updateRValues(Real C_Ci, Real C_Ck, Real C_Co, Real E_E250, Real R_E250, Real R_Rg, Real R_Ri, Real R_Rk, Real R_Ro, Real R_Vi, Real sampleRate, Triode& tube) {
+
+
+	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;
-		CiR = 1.0 / (2.0*C_Ci*sampleRate);
-		CkR = 1.0 / (2.0*C_Ck*sampleRate);
-		CoR = 1.0 / (2.0*C_Co*sampleRate);
-		Real RoR = R_Ro;
-		Real RgR = R_Rg;
+		Real CiR = 1.0 / (2.0*C_Ci*sampleRate);
 		Real RiR = R_Ri;
+		Real RgR = R_Rg;
+		Real RoR = R_Ro;
 		Real RkR = R_Rk;
-		Real E250R = R_E250;
-		E250E = E_E250;
+		Real CkR = 1.0 / (2.0*C_Ck*sampleRate);
+		Real E500R = R_E500;
+		E500E = E_E500;
+		Real CoR = 1.0 / (2.0*C_Co*sampleRate);
 		Real S0_3R = (CiR + ViR);
 		S0_3Gamma1 = CiR/(CiR + ViR);
 		Assert(S0_3Gamma1 >= 0.0 && S0_3Gamma1 <= 1.0);
@@ -45,82 +51,101 @@ public:
 		Assert(P0_3Gamma1 >= 0.0 && P0_3Gamma1 <= 1.0);
 		S1_3Gamma1 = RgR/(RgR + P0_3R);
 		Assert(S1_3Gamma1 >= 0.0 && S1_3Gamma1 <= 1.0);
-		Real I3_1R = CkR;
-		Real I3_2R = RkR;
-		I3_3Gamma1 = 1.0 / I3_1R/(1.0 / I3_1R + 1.0 / I3_2R);
-		Assert(I3_3Gamma1 >= 0.0 && I3_3Gamma1 <= 1.0);
-		Real S2_3R = (CoR + RoR);
-		S2_3Gamma1 = CoR/(CoR + RoR);
-		Assert(S2_3Gamma1 >= 0.0 && S2_3Gamma1 <= 1.0);
-		Real P2_1R = S2_3R;
-		Real P2_2R = E250R;
+		Real P1_1R = CkR;
+		Real P1_2R = RkR;
+		Real P1_3R = 1.0 /(1.0 / P1_1R + 1.0 / P1_2R);
+		P1_3Gamma1 = 1.0 / P1_1R/(1.0 / P1_1R + 1.0 / P1_2R);
+		Assert(P1_3Gamma1 >= 0.0 && P1_3Gamma1 <= 1.0);
+		Real S3_3R = (CoR + RoR);
+		S3_3Gamma1 = CoR/(CoR + RoR);
+		Assert(S3_3Gamma1 >= 0.0 && S3_3Gamma1 <= 1.0);
+		Real P2_1R = S3_3R;
+		Real P2_2R = E500R;
+		Real P2_3R = 1.0 /(1.0 / P2_1R + 1.0 / P2_2R);
 		P2_3Gamma1 = 1.0 / P2_1R/(1.0 / P2_1R + 1.0 / P2_2R);
 		Assert(P2_3Gamma1 >= 0.0 && P2_3Gamma1 <= 1.0);
-		t.Kr = sanitize_denormal(I3_3Gamma1);
-		t.Pr = sanitize_denormal(P2_3Gamma1);
-		t.Gr = sanitize_denormal(S1_3Gamma1);
-		//printf("Kr = %f  Pr = %f  Gr = %f\n", t.Kr, t.Pr, t.Gr);
-	}
-
-	void warmup_tubes(void) {
-		int i;
-		on = false;
-		for (i = 0; i < 8000; i++) {
-			advanc(0.0);
-		}
-		on = true;
-	}
+		S2_3Gamma1 = P1_3R/(P1_3R + P2_3R);
+		Assert(S2_3Gamma1 >= 0.0 && S2_3Gamma1 <= 1.0);
+	}	
 
-	Real advanc(Real VE){
-		ViE = VE;
+	Real advanc(Real ViE) {
+		//Get Bs
+		//P1_3GetB
 		Real Ckb = Cka;
-		Real I3_3b3 = I3_3Gamma1 * Ckb / CkR;
-		Real Cib = Cia;
-		Real S0_3b3 = (ViE + S0_3Gamma1 * Cib / CiR);
-		Real P0_3b3 = -P0_3Gamma1 * S0_3b3;
-		Real S1_3b3 = (P0_3b3 + S1_3Gamma1*(P0_3b3));
+		//P1_1SetA
+		//RkGetB
+		//P1_2SetA
+		Real P1_3b3 = -P1_3Gamma1*(-Ckb);
 		Real Cob = Coa;
-		Real S2_3b3 = Cob;
-		Real P2_3b3 = (E250E - P2_3Gamma1 * S2_3b3);
-		//Tube:    K       G      P
-		t.compute(I3_3b3,-S1_3b3,P2_3b3);
-		Real b1 = t.getC();
-		Real b2 = t.getG();
-		Real b3 = t.getP();
+		//S3_1SetA
+		//RoGetB
+		//S3_2SetA
+		Real S3_3b3 = -(Cob);
+		//P2_1SetA
+		//E500GetB
+		//P2_2SetA
+		Real P2_3b3 = E500E - P2_3Gamma1*(E500E - S3_3b3);
+		//S2_2SetA
+		Real S2_3b3 = -(P1_3b3 + P2_3b3);
+		//S1_3GetB
+		//RgGetB
+		//S1_1SetA
+		//P0_3GetB
+		//S0_3GetB
+		Real Cib = Cia;
+		//S0_1SetA
+		//ViGetB
+		//S0_2SetA
+		Real S0_3b3 = -(Cib + ViE);
+		//P0_1SetA
+		//RiGetB
+		//P0_2SetA
+		Real P0_3b3 = -P0_3Gamma1*(-S0_3b3);
+		//S1_2SetA
+		Real S1_3b3 = -(P0_3b3);
+		//Call tube model
+		Vg = -S1_3b3;
+		Real b = t.compute(-S2_3b3, S2_3Gamma1, Vg, Vk);
 		//Set As
-		Real I3_3b1 = (b1 - I3_3Gamma1*(b1 + Ckb));
-		Cka = I3_3b1;
-		Real S1_3b2 = (P0_3b3 - b2 - S1_3Gamma1*(P0_3b3 - b2));
-		Real P0_3b1 = (S1_3b2 + (-S0_3b3) - P0_3Gamma1*(S1_3b2 - S0_3b3));
-		Real S0_3b1 = (Cib - P0_3b1 - S0_3Gamma1*(Cib - P0_3b1));
-		Cia = S0_3b1;
-		Real P2_3b1 = (b3 - S2_3b3 - P2_3Gamma1*(b3 - S2_3b3));
-		Real S2_3b2 = (Cob - S2_3Gamma1*(Cob - P2_3b1));
-		Coa = S2_3b2;
-		Real S2_3b1 = (Cob - P2_3b1 - S2_3Gamma1*(Cob - P2_3b1));
-		Real Roa = S2_3b1;
-		//printf("K=%f G=%f P=%f out=%f\n", I3_3b3,-S1_3b3,P2_3b3, Roa);
-		return Roa;
+		//S2_3SetA
+		Real S2_3b1 = P1_3b3 - S2_3Gamma1*(P1_3b3 + P2_3b3 + b);
+		//P1_3SetA
+		Real P1_3b1 = S2_3b1  - Ckb - P1_3Gamma1*(-Ckb);
+		Cka = P1_3b1;
+		//RkSetA
+		Real S2_3b2 = P1_3b3 + b - S2_3Gamma1*(P1_3b3 + P2_3b3 + b);
+		//P2_3SetA
+		Real P2_3b1 = S2_3b2 + E500E - S3_3b3 - P2_3Gamma1*(E500E - S3_3b3);
+		//S3_3SetA
+		Real S3_3b1 = Cob - S3_3Gamma1*(Cob + P2_3b1);
+		Coa = S3_3b1;
+		Real S3_3b2 = Cob + P2_3b1 - S3_3Gamma1*(Cob + P2_3b1);
+		//RoSetA
+		Real Roa = S3_3b2;
+		Vk = -P1_3b3;
+		//printf("K=%f G=%f P=%f  b=%f  in=%f out=%f\n", P1_3b3,S1_3b3,P2_3b3,b, ViE,-Roa);
+		return -(Roa);
 	}
 
 private:
 	//State variables
+	Real Coa;
 	Real Cia;
 	Real Cka;
-	Real Coa;
+	Real Vk;
+	Real Vg;
 
 	//R values
-	Real CiR;
-	Real CkR;
-	Real CoR;
-	Real I3_3Gamma1;
-	Real E250E;
-	Real ViE;
+	Real P3_3Gamma1;
 	Real P0_3Gamma1;
+	Real S1_3Gamma1;
+	Real S3_3Gamma1;
+	Real P1_3Gamma1;
+	Real ViE;
 	Real S0_3Gamma1;
 	Real S2_3Gamma1;
-	Real S1_3Gamma1;
 	Real P2_3Gamma1;
+	Real E500E;
 };
 
 #endif