WIP approximation for zamtubeapproxtube

2 files changed, 70 insertions, 15 deletions

diff --git a/plugins/ZamTube/lut_current.m b/plugins/ZamTube/lut_current.m
index 3d42d6e..6546e61 100644
--- a/plugins/ZamTube/lut_current.m
+++ b/plugins/ZamTube/lut_current.m
@@ -1,10 +1,10 @@
 1;
 
 function I = ip (vgk, vpk)
-  mu=103.2
-  kx=1.26
+  mu=100.0
+  kx=1.4
   kg1=446.0
-  kp=3.4
+  kp=600.0
   kvb=300.0
   e1=vpk .* log(1. + exp(kp .* (1.0/mu + vgk ./ (kvb + vpk .* vpk) .^ 0.5))) / kp
   if e1 < 0
@@ -14,12 +14,33 @@ function I = ip (vgk, vpk)
   I = double((e1 .^ kx)/kg1)
 endfunction
 
-gridsizeg = 2001
-gridsizep = 2001
-vgk = linspace(-20, 20, gridsizeg)
-vpk = linspace(0, 500, gridsizep)
+function f1 (x,y)
+f1 = 1./446 .* (y./sqrt(300).-3.*exp(1).*x.*sqrt(300)./100).^1.4
+endfunction
+
+function f2 (x,y)
+f2 = 1./446.*(y.*log(1.+6.*exp(1))./sqrt(300)).^1.4
+endfunction
+
+function f4 (x,y)
+if (-x./y < 17)
+  f4 = f1(x,y)
+  return
+endif
+if ((x > -0.1) && (y > 250))
+  f4 = f2(x,y)
+  return
+endif
+f4 = ip(x,y)
+endfunction
+
+gridsizeg = 101
+gridsizep = 101
+vgk = linspace(-10, 0, gridsizeg)
+vpk = linspace(0, 300, gridsizep)
 [vg, vp] = ndgrid(double(vgk), double(vpk))
-cur = 1e+6*ip(double(vg), double(vp))
+cur = 1000000*ip(double(vg), double(vp))
+approxcur = 1000000*f4(vg, vp)
 surf(vg,vp,cur)
 
 more off
diff --git a/plugins/ZamTube/triode.cpp b/plugins/ZamTube/triode.cpp
index 9796063..ec42960 100644
--- a/plugins/ZamTube/triode.cpp
+++ b/plugins/ZamTube/triode.cpp
@@ -31,13 +31,47 @@ T Triode::getIa(T Vgk, T Vpk) const {
 		Vgk = 0.0;
 	}
 
-	/* exact solution (expensive) */
-	T ee1 = Vpk*log1p(exp(kp*(1./mu+Vgk/sqrt(kvb+Vpk*Vpk))))/kp;
-	if (ee1 < 0) {
-		return 0.;
-	}
-	//printf("Vpk=%f ans=%f e1=%f exact_e1=%f\n", Vpk, ans, e1, ee1);
-	return 1e+6*pow(ee1, kx) / kg1;
+	/* f(x,y) = (y*log(1+kp/mu*exp(1+mu*x/sqrt(kvb+y*y))))^kx
+	 *        = (y*log(1+6*exp(1+100*x/sqrt(300+y*y))))^1.4
+	 * 
+	 * Let y' = y * sqrt(300)
+	 * Let x' = x * 100 / sqrt(300)
+	 *
+	 * f(x',y') = (y'*log(1+6*exp(1+x'/sqrt(1+y'*y'))))^1.4
+	 *
+	 * When -x'/y' < 17
+	 * 	f(x',y') ~= (y' - 3*exp(1)*x')^1.4
+	 * 	f(x,y) ~= (y/sqrt(300) - 3*exp(1)*x*sqrt(300)/100)^1.4
+	 *
+	 * When |x'| is small && y' is large
+	 * 	f(x',y') ~= (y'*log(1 + 6*exp(1)))^1.4
+	 * 	f(x,y) ~=(y/sqrt(300)*log(1 + 6*exp(1)))^1.4
+	 *
+	 * Otherwise, use exact solution
+	 */
+	float f = 0.f;
+	float sc = 1e+6 / kg1;
+/*
+#define G_NEG_SMALL	-1.f
+#define P_LARGE		150.f
+	float e1 = expf(1.f);
+	float r300 = sqrtf(300.f);
+
+	if (-Vgk/Vpk < 17.) {
+		f = sc * powf(Vpk / r300 - 3.f * e1 * Vgk * r300 / 100., 1.4f);
+	} else if ((Vgk > G_NEG_SMALL) && (Vpk > P_LARGE)) {
+		f = sc * powf(Vpk * log1pf(6.f*e1) / r300, 1.4f);
+	} else {
+*/
+		/* exact solution (expensive) */
+		float ee1 = Vpk*log1pf(expf(kp*(1./mu+Vgk/sqrtf(kvb+Vpk*Vpk))))/kp;
+		if (ee1 < 0) {
+			return 0.;
+		}
+		f = sc * powf(ee1, kx);
+//	}
+	//printf("Vpk=%f e1=%f exact_e1=%f\n", Vpk, ans, e1, ee1);
+	return f;
 }
 
 T Triode::iterateNewtonRaphson(T x, T dx, T Vgk, T a, T R) const {