## Copyright (C) 2015, 2023 Stefano D'Angelo
##
## Permission is hereby granted, free of charge, to any person obtaining a copy
## of this software and associated documentation files (the "Software"), to deal
## in the Software without restriction, including without limitation the rights
## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
## copies of the Software, and to permit persons to whom the Software is
## furnished to do so, subject to the following conditions:
##
## The above copyright notice and this permission notice shall be included in
## all copies or substantial portions of the Software.
##
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
## OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
## THE SOFTWARE.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{y} =} triode_stage (@var{fs}, @var{x}, @var{G}, @var{u}, @var{h}, @var{D}, @var{K}, @var{Voff}, @var{E}, @var{Rp}, @var{Rk}, @var{Ck}, @var{Rg}, @var{gc})
##
## Simulate the common-cathode triode stage circuit described in
## @quotation
## S. D'Angelo, J. Pakarinen, and V. Valimaki, ``New Family of Wave-Digital
## Triode Models,'' @emph{IEEE Trans. Audio, Speech, and Lang. Process.},
## vol. 21, no. 2, pp. 313-321, February 2013.
## @end quotation
##
## It produces the output vector @var{y} from the audio input vector @var{x}.
## @var{fs} is the sample rate (in Hz).
##
## @var{G}, @var{u}, and @var{h} are vectors of 4 real values representing
## polynomial coefficients for the respective tube parameter expressions.
##
## @var{D}, @var{K}, and @var{Voff} are scalars representing the remaining
## tube parameters.
##
## @var{E}, @var{Rp}, @var{Rk}, @var{Ck}, and @var{Rg} are scalars representing
## circuit component values.
##
## If @var{gc} is true, the grid current effect is simulated, otherwise it is
## not.
##
## @end deftypefn

## Author: Stefano D'Angelo <zanga.mail@gmail.com>
## Maintainer: Stefano D'Angelo <zanga.mail@gmail.com>
## Version: 1.1.0
## Keywords: vacuum tube triode guitar amplifier

function y = triode_stage(fs, x, G, u, h, D, K, Voff, E, Rp, Rk, Ck, Rg, gc)

  #### Input arguments

  # Checks

  if (nargin != 14)
    print_usage();
  elseif (!isscalar(fs) || !isnumeric(fs) || !isreal(fs) || fs <= 0)
    error("triode_stage: FS must be a positive real");
  elseif (isscalar(x) || !isvector(x) || !isnumeric(x) || !isreal(x))
    error("triode_stage: X must be a vector of real values");
  elseif (!isvector(G) || !isnumeric(G) || !isreal(G) || length(G) != 4)
    error("triode_stage: G must be a vector of 4 real values");
  elseif (!isvector(u) || !isnumeric(u) || !isreal(u) || length(u) != 4)
    error("triode_stage: U must be a vector of 4 real values");
  elseif (!isvector(h) || !isnumeric(h) || !isreal(h) || length(h) != 4)
    error("triode_stage: H must be a vector of 4 real values");
  elseif (!isscalar(D) || !isnumeric(D) || !isreal(D))
    error("triode_stage: D must be a real scalar");
  elseif (!isscalar(K) || !isnumeric(K) || !isreal(K))
    error("triode_stage: K must be a real scalar");
  elseif (!isscalar(Voff) || !isnumeric(Voff) || !isreal(Voff))
    error("triode_stage: VOFF must be a real scalar");
  elseif (!isscalar(E) || !isnumeric(E) || !isreal(E))
    error("triode_stage: E must be a real scalar");
  elseif (!isscalar(Rp) || !isnumeric(Rp) || !isreal(Rp))
    error("triode_stage: RP must be a real scalar");
  elseif (!isscalar(Rk) || !isnumeric(Rk) || !isreal(Rk))
    error("triode_stage: RK must be a real scalar");
  elseif (!isscalar(Ck) || !isnumeric(Ck) || !isreal(Ck))
    error("triode_stage: CK must be a real scalar");
  elseif (!isscalar(Rg) || !isnumeric(Rg) || !isreal(Rg))
    error("triode_stage: RG must be a real scalar");
  endif

  # Enforce bounds

  E  = limit(E,  1e-6,  inf, "E");
  Rp = limit(Rp, 1e-6,  inf, "Rp");
  Rk = limit(Rk, 1e-6,  inf, "Rk");
  Ck = limit(Ck, 1e-18, inf, "Ck");

  #### Constants

  # Circuit components

  Ri = 1e6;
  Ci = 100e-9;
  Ro = 1e6;
  Co = 10e-9;
  
  # Initial voltages
  
  Vk0 = 0.5 * Rk / (Rp + Rk);
  Vp0 = 0.5 * (E + Rk / (Rp + Rk));
  function v = f_init_cond(x)
    Gmin = 1e-9;
    umin = 1e-9;
    
    Vp = x(1);
    Vk = x(2);
    Vpk = Vp - Vk;
    Vgk = -Vk;
    Vgk2 = Vgk * Vgk;
    Vgk3 = Vgk2 * Vgk;
    
    xG = G(1) + G(2) * Vgk + G(3) * Vgk2 + G(4) * Vgk3;
    xu = u(1) + u(2) * Vgk + u(3) * Vgk2 + u(4) * Vgk3;
    xh = h(1) + h(2) * Vgk + h(3) * Vgk2 + h(4) * Vgk3;

    if (xG < Gmin)
      xG = Gmin;
    endif
    if (xu < umin)
      xu = umin;
    endif
    
    Ik = xG * (Vgk + Vpk / xu + xh) ^ 1.5;
    if (Vgk <= Voff || !gc)
      Ig = 0;
    else
      Ig = Ik / (1 + D * (Vpk / (Vgk - Voff)) ^ K);
    endif
    Ip = Ik + Ig;
    
    v = [Ip - (E - Vp) / Rp, Ik - Vk / Rk];
  endfunction
  v = fsolve(@f_init_cond, [Vp0, Vk0], optimset("AutoScaling", "on"));
  Vp0 = v(1);
  Vk0 = v(2);

  # WDF element values

  wVi_R  = 1e-6;
  wCi_R  = 1 / (2 * fs * Ci);
  wCk_R  = 1 / (2 * fs * Ck);
  wCo_R  = 1 / (2 * fs * Co);
  wsi_kl = wCi_R / (wCi_R + wVi_R);
  wsi_R  = wCi_R + wVi_R;
  wpg_kt = wsi_R / (wsi_R + Ri);
  wpg_R  = (wsi_R * Ri) / (wsi_R + Ri);
  wsg_kl = Rg / (Rg + wpg_R);
  wsg_R  = Rg + wpg_R;
  wpk_kt = wCk_R / (Rk + wCk_R);
  wpk_R  = (Rk * wCk_R) / (Rk + wCk_R);
  wsp_kl = wCo_R / (wCo_R + Ro);
  wsp_R  = wCo_R + Ro;
  wpp_kt = wsp_R / (wsp_R + Rp);
  wpp_R  = (wsp_R * Rp) / (wsp_R + Rp);

  # Filter coefficients

  kTxCi  = 1 - wpg_kt;
  kTCk   = 1 - wpk_kt;
  kTCo   = 1 - wpp_kt;
  kT0    = wpp_kt * E;
  kyT    =  0.5 * (1 - wsp_kl);
  kyCo   = -0.5 * (1 - wsp_kl) * (1 + wpp_kt);
  ky0    =  0.5 * (1 - wsp_kl) * wpp_kt * E;
  kCiT   = wsi_kl * (1 - wsg_kl);
  kCixCi = wsi_kl * ((1 - wpg_kt) * (wsg_kl + 1) - 2);
  kCoCo  = 1 - wsp_kl * (1 + wpp_kt);
  kCo0   = wsp_kl * wpp_kt * E;

  #### Filter

  y     = zeros(1, length(x));
  wCi_s = 0;
  wCk_s = Vk0;
  wCo_s = Vp0;

  for i = 1:length(x)

    ### Grid subcircuit (upwards)

    ## Vi, Ci, Rg, Ri outs
    # wVi_b = x(i)
    # wCi_b = wCi_s
    # wRg_b = 0
    # wRi_b = 0

    ## Series junction i 
    # wsi_al = wCi_b = wCi_s
    # wsi_ar = wVi_b = x(i)
    # wsi_bt = -(wsi_al + wsi_ar) = -(x(i) + wCi_s)

    ## Polarity inverter i
    # wii_ab = wsi_bt  = -(x(i) + wCi_s)
    # wii_bt = -wii_ab =   x(i) + wCi_s

    ## Parallel junction g
    # wpg_al = wii_bt = x(i) + wCi_s
    # wpg_ar = wRi_b  = 0
    # wpg_bt = wpg_al + wpg_kt * (wpg_ar - wpg_al) ...
    #        = (1 - wpg_kt) * (x(i) + wCi_s)

    ## Series junction g
    # wsg_al = wRg_b  = 0
    # wsg_ar = wpg_bt = (1 - wpg_kt) * (x(i) + wCi_s)
    # wsg_bt = -(wsg_al + wsg_ar) = (wpg_kt - 1) * (x(i) + wCi_s)

    ## Polarity inverter g
    # wig_ab =  wsg_bt = (wpg_kt - 1) * (x(i) + wCi_s)
    # wig_bt = -wig_ab = (1 - wpg_kt) * (x(i) + wCi_s)

    ### Cathode subcircuit (upwards)

    ## Rk, Ck outs
    # wRk_b = 0
    # wCk_b = wCk_s

    ## Parallel junction k
    # wpk_al = wCk_b = wCk_s
    # wpk_ar = wRk_b = 0
    # wpk_bt = wpk_al + wpk_kt * (wpk_ar - wpk_al) = (1 - wpk_kt) * wCk_s

    ### Plate subcircuit (upwards)

    ## ERp, Ro, Co outs
    # wERp_b = E
    # wRo_b  = 0
    # wCo_b  = wCo_s

    ## Series junction p
    # wsp_al = wCo_b = wCo_s
    # wsp_ar = wRo_b = 0
    # wsp_bt = -(wsp_al + wsp_ar) = -wCo_s

    ## Polarity inverter p
    # wip_ab =  wsp_bt = -wCo_s
    # wip_bt = -wip_ab =  wCo_s

    ## Parallel junction p
    # wpp_al = wip_bt  = wCo_s
    # wpp_ar = wERp0_b = E
    # wpp_bt = wpp_al + wpp_kt * (wpp_ar - wpp_al) ...
    #        = wCo_s + wpp_kt * (E - wCo_s)

    ### Triode

    # wT_ag = wig_bt = (1 - wpg_kt) * (x(i) + wCi_s)
    # wT_ak = wpk_bt = (1 - wpk_kt) * wCk_s
    # wT_ap = wpp_bt = wCo_s + wpp_kt * (E - wCo_s)
    # [wT_bg, wT_bk, wT_bp] = triode(wT_ag, wT_ak, wT_ap, wsg_R, wpk_R, wpp_R,
    #                                G, u, h, D, K, Voff, gc)

    ### Plate subcircuit (downwards)

    ## Parallel junction p
    # wpp_at = wT_bp
    # wpp_bx = wpp_at + wpp_bt = wT_bp + wCo_s + wpp_kt * (E - wCo_s)
    # wpp_bl = wpp_bx - wpp_al = wT_bp + wpp_kt * (E - wCo_s)
    # wpp_br = wpp_bx - wpp_ar = wT_bp - (1 - wpp_kt) * (E - wCo_s)

    ## Polarity inverter p
    # wip_at =  wpp_bl = wT_bp + wpp_kt * (E - wCo_s)
    # wip_bb = -wip_at = wpp_kt * (wCo_s - E) - wT_bp

    ## Series junction p
    # wsp_at = wip_bb = wpp_kt * (wCo_s - E) - wT_bp
    # wsp_bl = wsp_al + wsp_kl * (wsp_bt - wsp_at) ...
    #        = wCo_s + wsp_kl * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)
    # wsp_br = -(wsp_bl + wsp_at) ...
    #        = (1 - wsp_kl) * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)

    ## ERp, Ro, Co ins
    # wERp_a = wpp_br = wT_bp - (1 - wpp_kt) * (E - wCo_s)
    # wRo_a  = wsp_br = (1 - wsp_kl) * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)
    # wCo_a  = wsp_bl = wCo_s + wsp_kl * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)

    ### Cathode subcircuit (downwards)

    ## Parallel junction k
    # wpk_at = wT_bk
    # wpk_bx = wpk_at + wpk_bt = wT_bk + (1 - wpk_kt) * wCk_s
    # wpk_bl = wpk_bx - wpk_al = wT_bk - wpk_kt * wCk_s
    # wpk_br = wpk_bx - wpk_ar = wT_bk + (1 - wpk_kt) * wCk_s

    ## Rk, Ck ins
    # wRk_a = wpk_br = wT_bk + (1 - wpk_kt) * wCk_s
    # wCk_a = wpk_bl = wT_bk - wpk_kt * wCk_s

    ### Grid subcircuit (downwards)

    ## Polarity inverter g
    # wig_at = wT_bg
    # wig_bb = -wig_at = -wT_bg

    ## Series junction g
    # wsg_at = wig_bb = -wT_bg
    # wsg_bl = wsg_al + wsg_kl * (wsg_bt - wsg_at) ...
    #        = wsg_kl * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))
    # wsg_br = -(wsg_bl + wsg_at) ...
    #        = wT_bg - wsg_kl * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))

    ## Parallel junction g
    # wpg_at = wsg_br = wT_bg - wsg_kl * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))
    # wpg_bx = wpg_at + wpg_bt ...
    #        = 2 * wT_bg ...
    #          - (1 + wsg_kl) * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))
    # wpg_bl = wpg_bx - wpg_al ...
    #        = (1 - wsg_kl) * wT_bg ...
    #          + ((1 - wpg_kt) * wsg_kl - wpg_kt) * (x(i) + wCi_s)
    # wpg_br = wpg_bx - wpg_ar ...
    #        = 2 * wT_bg ...
    #          - (1 + wsg_kl) * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))

    ## Polarity inverter i
    # wii_at =  wpg_bl ...
    #        = (1 - wsg_kl) * wT_bg ...
    #          + ((1 - wpg_kt) * wsg_kl - wpg_kt) * (x(i) + wCi_s)
    # wii_bb = -wii_at ...
    #        = (wsg_kl - 1) * wT_bg ...
    #          + (wpg_kt - (1 - wpg_kt) * wsg_kl) * (x(i) + wCi_s)

    ## Series junction i
    # wsi_at = wii_bb ...
    #        = (wsg_kl - 1) * wT_bg ...
    #          + (wpg_kt - (1 - wpg_kt) * wsg_kl) * (x(i) + wCi_s)
    # wsi_bl = wsi_al + wsi_kl * (wsi_bt - wsi_at) ...
    #        = wsi_kl * ((1 - wsg_kl) * wT_bg ...
    #                    + ((1 - wpg_kt) * (1 + wsg_kl) - 2) ...
    #                    * (x(i) + wCi_s)) + wCi_s
    # wsi_br unused

    ## Vi, Ci, Rg, Ri ins
    # wVi_a unused
    # wCi_a = wsi_bl ...
    #       = wsi_kl * ((1 - wsg_kl) * wT_bg ...
    #                   + ((1 - wpg_kt) * (wsg_kl + 1) - 2) ...
    #                   * (wCi_s + x(i))) + wCi_s
    # wRg_a = wsg_bl ...
    #        = wsg_kl * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))
    # wRi_a = wpg_br ...
    #       = 2 * wT_bg - (1 + wsg_kl) * (wT_bg - (1 - wpg_kt) * (x(i) + wCi_s))

    ### Global output

    # Vo = 0.5 * (wRo_a + wRo_b) ...
    #    = 0.5 * (1 - wsp_kl) * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)

    ### State updates

    # wCi_s = wCi_a ...
    #       = wsi_kl * ((1 - wsg_kl) * wT_bg ...
    #                   + ((1 - wpg_kt) * (wsg_kl + 1) - 2) ...
    #                   * (wCi_s + x(i))) + wCi_s
    # wCk_s = wCk_a = wpk_bl = wT_bk - wpk_kt * wCk_s
    # wCo_s = wCo_a = wCo_s + wsp_kl * (wT_bp + wpp_kt * (E - wCo_s) - wCo_s)

    ### Summing up

    xCi   = x(i) + wCi_s;
    wT_ag = kTxCi * xCi;
    wT_ak = kTCk * wCk_s;
    wT_ap = kTCo * wCo_s + kT0;
    [wT_bg, wT_bk, wT_bp] = triode(wT_ag, wT_ak, wT_ap, wsg_R, wpk_R, wpp_R,
                                   G, u, h, D, K, Voff, gc);
    y(i)  = kyT * wT_bp + kyCo * wCo_s + ky0;
    wCi_s = kCiT * wT_bg + kCixCi * xCi + wCi_s;
    wCk_s = wT_bk - wpk_kt * wCk_s;
    wCo_s = wsp_kl * wT_bp + kCoCo * wCo_s + kCo0;

  endfor

endfunction

function [bg, bk, bp] = triode(ag, ak, ap, R0g, R0k, R0p, G, u, h, D, K, Voff,
                               gc)
  max_err_k = 1e-5;
  max_err_g = 1e-5;
  dfVkmin   = 1e-9;
  dfVgmin   = 1e-9;

  ## Vk with no grid current

  Vg      = ag;
  Vkprev  = ak;
  fVkprev = triode_fVk(Vg, Vkprev, ag, ak, ap, R0g, R0k, R0p, G, u, h);
  Vk      = ak + abs(fVkprev);
  do
    fVk  = triode_fVk(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, G, u, h);
    dfVk = fVk - fVkprev;
    if (abs(dfVk) <= dfVkmin)
      break;
    endif

    Vknext = Vk - fVk * (Vk - Vkprev) / dfVk;
    if (Vknext < ak)
      # The original algorithm would limit the value and break. This solution
      # gives higher precision (but higher computational load too).
      Vknext = ak;
    endif

    Vkprev  = Vk;
    Vk      = Vknext;
    fVkprev = fVk;
    err_k   = Vk - Vkprev;
  until (abs(err_k) <= max_err_k)

  ## Grid current

  if (gc && (Vg - Vk > Voff))
    max_iter = 5;
    err_k = ak - Vk;
    do

      # Compute Vg

      Vgmin   = ag + R0g / R0k * (ak - Vk);
      Vgstart = Vg;
      Vgprev  = Vg;
      fVgprev = triode_fVg(Vgprev, Vk, ag, ak, ap, R0g, R0k, R0p, D, K, Voff);
      # First, see Errata. Then, again, the following is another accuracy
      # improvement similar to what was done before (limiting vs breaking) ...
      if (Vg == ag)
        Vg = Vg - abs(fVgprev);
      elseif (Vg == Vgmin)
        Vg = Vg + abs(fVgprev);
      else
        Vg = Vg + fVgprev;
        if (Vg > ag)
          Vg = ag;
        elseif (Vg < Vgmin)
          Vg = Vgmin;
        endif
      endif
      err_g   = Vg - Vgprev;
      do
        fVg = triode_fVg(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, D, K, Voff);
        dfVg = fVg - fVgprev;
        if (abs(dfVg) <= dfVgmin)
          break;
        endif

        Vgnext = Vg - fVg * (Vg - Vgprev) / dfVg;
        # ... and again...
        if (Vgnext > ag)
          Vgnext = ag;
        elseif (Vgnext < Vgmin)
          Vgnext = Vgmin;
        endif

        Vgprev  = Vg;
        Vg      = Vgnext;
        fVgprev = fVg;
        err_g   = Vg - Vgprev;
      until (abs(err_g) <= max_err_g)
      err_g = Vgstart - Vg;  # See Errata

      if (abs(err_k) <= max_err_k && abs(err_g) <= max_err_g)
        break;
      endif

      # Compute Vk

      Vkmin   = ak + R0k / R0g * (ag - Vg);
      Vkstart = Vk;
      Vkprev  = Vk;
      fVkprev = triode_fVk(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, G, u, h);
      # ... and again...
      if (Vk == ak || Vk == Vkmin)
        Vk = Vk + abs(fVkprev);
      else
        Vk = Vk + fVkprev;
        if (Vk < ak)
          Vk = ak;
        endif
        if (Vk < Vkmin)
          Vk = Vkmin;
        endif
      endif
      err_k   = Vk - Vkprev;
      do
        fVk = triode_fVk(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, G, u, h);
        dfVk = fVk - fVkprev;
        if (abs(dfVk) <= dfVkmin)
          break;
        endif

        Vknext = Vk - fVk * (Vk - Vkprev) / dfVk;
        # ... and again.
        if (Vknext < ak)
          Vknext = ak;
        endif
        if (Vknext < Vkmin)
          Vknext = Vkmin;
        endif

        Vkprev  = Vk;
        Vk      = Vknext;
        fVkprev = fVk;
        err_k   = Vk - Vkprev;
      until (abs(err_k) <= max_err_k)
      err_k = Vkstart - Vk;  # See Errata
      max_iter--;
    until (abs(err_k) <= max_err_k && abs(err_g) <= max_err_g || max_iter == 0)
  endif

  ## Compute outputs

  Vkak = Vk - ak;
  Vgag = Vg - ag;
  bg   = Vg + Vgag;
  bk   = Vk + Vkak;
  bpx  = R0p * (Vgag / R0g + Vkak / R0k);
  bp   = ap - bpx - bpx;
endfunction

function fVk = triode_fVk(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, G, u, h)
  Gmin = 1e-9;
  umin = 1e-9;

  agVg = ag - Vg;

  Vgk  = Vg - Vk;
  Vgk2 = Vgk * Vgk;
  Vgk3 = Vgk2 * Vgk;

  G = G(1) + G(2) * Vgk + G(3) * Vgk2 + G(4) * Vgk3;
  u = u(1) + u(2) * Vgk + u(3) * Vgk2 + u(4) * Vgk3;
  h = h(1) + h(2) * Vgk + h(3) * Vgk2 + h(4) * Vgk3;

  if (G < Gmin)
    G = Gmin;
  endif
  if (u < umin)
    u = umin;
  endif

  alpha = (Vk - ak) / (G * R0k);
  if (alpha < 0)
    alpha = 0;
  else
    alpha = alpha ^ (2 / 3);
  endif

  # See Errata
  fVk = (R0g * (R0p * ak + R0k * (ap + u * (Vg + h - alpha))) ...
         + R0k * R0p * agVg) / (R0g * (R0p + (u - 1.0) * R0k)) - Vk;
endfunction

function fVg = triode_fVg(Vg, Vk, ag, ak, ap, R0g, R0k, R0p, D, K, Voff)
  agVg = ag - Vg;
  akVk = ak - Vk;

  if (agVg == 0)
    fVg = Voff + Vk - Vg;
    return;
  endif

  t = -(R0g / R0k * akVk / agVg + 1.0) / D;
  if (t < 0)
    fVg = (R0k * (R0g * (ap - Vk) + R0p * ag) + R0g * R0p * akVk) ...
          / (R0k * R0p) - Vg;
  else
    t = t ^ (1.0 / K);
    if (isinf(t))
      fVg = Voff + Vk - Vg;
    else
      fVg = (R0k * (R0g * ((Voff + Vk) * t - Vk + ap) + R0p * ag) ...
             + R0g * R0p * akVk) / (R0k * (R0g * t + R0p)) - Vg;
    endif
  endif
endfunction

function y = limit(x, low, up, name)
  y = x;
  if (min(y) < low || max(y) > up)
    warning(["triode_stage: limiting " name " to [" num2str(low) ", " ...
             num2str(up) "]"]);
    y(y < low) = low;
    y(y > up)  = up;
  endif
endfunction
