diff options
author | Damien Zammit <damien@zamaudio.com> | 2019-01-03 17:18:26 +1100 |
---|---|---|
committer | Damien Zammit <damien@zamaudio.com> | 2019-01-06 21:49:16 +1100 |
commit | 63f36b6bdd937d61a24a4e8cb7030fb53a7b0189 (patch) | |
tree | d2feda0f3ef79103b20512886b23967922c8069a /plugins | |
parent | 96ec0c7dbc9c034bc8716309b92fda84bb27cf92 (diff) |
Simpler circuit, amplification is good
Diffstat (limited to 'plugins')
-rw-r--r-- | plugins/ZamTube/ZamTubePlugin.cpp | 43 | ||||
-rw-r--r-- | plugins/ZamTube/ZamTubePlugin.hpp | 5 | ||||
-rw-r--r-- | plugins/ZamTube/triode.cpp | 657 | ||||
-rw-r--r-- | plugins/ZamTube/triode.h | 33 | ||||
-rw-r--r-- | plugins/ZamTube/wdfcircuits.h | 167 |
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 |