## Copyright (C) 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_quadric (@var{fs}, @var{x}, @var{kp}, @var{kp2}, @var{kpg}, @var{E}, @var{Rp}, @var{Rk}, @var{Ck}, @var{Rg})
##
## Simulate the common-cathode triode stage circuit described in
## @quotation
## R. Giampiccolo, S. D'Angelo, A. Bernardini, and A. Sarti, ``A Quadric Surface
## Model of Vacuum Tubes for Virtual Analog Applications,''
## @emph{submitted paper}.
## @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{kp}, @var{kp2}, and @var{kpg} are scalars representing the tube
## parameters.
##
## @var{E}, @var{Rp}, @var{Rk}, @var{Ck}, and @var{Rg} are scalars representing
## circuit component values.
##
## @end deftypefn

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

function y = triode_stage_quadric(fs, x, kp, kp2, kpg, E, Rp, Rk, Ck, Rg)

  #### Input arguments

  # Checks

  if (nargin != 10)
    print_usage();
  elseif (!isscalar(fs) || !isnumeric(fs) || !isreal(fs) || fs <= 0)
    error("triode_stage_quadric: FS must be a positive real");
  elseif (isscalar(x) || !isvector(x) || !isnumeric(x) || !isreal(x))
    error("triode_stage_quadric: X must be a vector of real values");
  elseif (!isscalar(kp) || !isnumeric(kp) || !isreal(kp))
    error("triode_stage_quadric: KP must be a real scalar");
  elseif (!isscalar(kp2) || !isnumeric(kp2) || !isreal(kp2))
    error("triode_stage_quadric: KP2 must be a real scalar");
  elseif (!isscalar(kpg) || !isnumeric(kpg) || !isreal(kpg))
    error("triode_stage_quadric: KPG must be a real scalar");
  elseif (!isscalar(E) || !isnumeric(E) || !isreal(E))
    error("triode_stage_quadric: E must be a real scalar");
  elseif (!isscalar(Rp) || !isnumeric(Rp) || !isreal(Rp))
    error("triode_stage_quadric: RP must be a real scalar");
  elseif (!isscalar(Rk) || !isnumeric(Rk) || !isreal(Rk))
    error("triode_stage_quadric: RK must be a real scalar");
  elseif (!isscalar(Ck) || !isnumeric(Ck) || !isreal(Ck))
    error("triode_stage_quadric: CK must be a real scalar");
  elseif (!isscalar(Rg) || !isnumeric(Rg) || !isreal(Rg))
    error("triode_stage_quadric: 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
  
  k1 = kpg / (2 * kp2) + Rp / Rk + 1;
  k2 = k1 * (kp / kp2 + 2 * E) * kp2;
  k3 = Rk * k2 + 1;
  Vk0 = (k3 - sign(k1) * sqrt(2 * k3 - 1)) / (2 * Rk * k1 * k1 * kp2);
  Vp0 = E - Rp / Rk * Vk0;

  # 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,
                                   kp, kp2, kpg);
    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, kp, kp2, kpg)
  # These could be precomputed
  bk_bp = R0k / R0p;
  k_eta = 1 / (bk_bp * (0.5 * kpg + kp2) + kp2); # = 1 / (2 * kp2 * gamma)
  k_delta = kp2 * k_eta * k_eta / (R0p + R0p);
  k_bp_s = k_eta * sqrt((kp2 + kp2) / R0p);
  bp_k = 1 / (R0p + R0k);
  bp_ap_0 = bp_k * (R0k - R0p);
  bp_ak_0 = bp_k * (R0p + R0p);
  
  v1 = 0.5 * ap;
  v2 = ak + v1 * bk_bp;
  alpha = kpg * (ag - v2) + kp;
  beta = kp2 * (v1 - v2);
  eta = k_eta * (beta + beta + alpha);
  v3 = eta + k_delta;
  delta = ap + v3;
  
  if (delta >= 0)
    bp = k_bp_s * sqrt(delta) - v3 - k_delta;
    d = bk_bp * (ap - bp);
    bk = ak + d;
    Vpk2 = ap + bp - ak - bk;
    
    if (kpg * (ag - ak - 0.5 * d) + kp2 * Vpk2 + kp < 0)
      bp = ap;
      bk = ak;
      Vpk = ap - ak;
    else
      Vpk = 0.5 * Vpk2;
    endif
  else
    bp = ap;
    bk = ak;
    Vpk = ap - ak;
  endif
  
  if (Vpk < 0)
    bp = bp_ap_0 * ap + bp_ak_0 * ak;
  endif
  
  bg = ag;
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
