function currents = transistor_model(Vbec, model_parms)
% Model of a bipolar transistor, returns currents object with 3 members:
%   currents.Base = base current (A)
%   currents.Emitter = emitter current (A)
%   currents.Collector = collector current (A)
% in terms of one independent variable object with 3 members:
%   Vbec.Base = base potential (V)
%   Vbec.Emitter = emitter potential (V)
%   Vbec.Collector = collector potential (V)

% The model parameters are those of the SPICE BJT model:
%   model_parms.IS = Transport saturation current (A)
%   model_parms.BF = Ideal maximum forward beta
%   model_parms.NF = Forward current emission coefficient
%   model_parms.VAF = Forward early voltage (V)
%   model_parms.IKF = Corner for forward beta high current roll-off (A)
%   model_parms.ISE = B-E leakage saturation current (A)
%   model_parms.NE = B-E leakage emission coefficient
%   model_parms.BR = Ideal maximum reverse beta
%   model_parms.NR = Reverse current emission coefficient
%   model_parms.VAR = Reverse early voltage (V)
%   model_parms.IKR = Corner for reverse beta high current roll-off (A)
%   model_parms.ISC = B-C leakage saturation current (A)
%   model_parms.NC = B-C leakage emission coefficient
%   model_parms.RB = Base ohmic resistance (Ohms)
%   model_parms.IRB = Current where base resistance falls halfway to its min value (A)
%   model_parms.RBM = Minimum base resistance at high currents (Ohms)
%   model_parms.RE = Emitter resistance (Ohms)
%   model_parms.RC = Collector resistance (Ohms)
%   model_parms.CJE = B-E zero-bias depletion capacitance (F)
%   model_parms.VJE = B-E built-in potential (V)
%   model_parms.MJE = B-E junction exponential factor
%   model_parms.TF = Ideal forward transit time (s)
%   model_parms.XTF = Coefficient for bias dependence on TF
%   model_parms.VTF = Voltage describing VBC dependence on TF (V)
%   model_parms.ITF = High-current parameter for effect on TF (A)
%   model_parms.PTF = Excess phase at freq = 1/(2 pi TF) (degrees)
%   model_parms.CJC = B-C zero-bias depletion capacitance (F)
%   model_parms.VJC = B-C built-in potential (V)
%   model_parms.MJC = B-C junction exponential factor
%   model_parms.XCJC = Fraction of B-C depletion capacitance connected to internal base node
%   model_parms.TR = Ideal reverse transit time (s)
%   model_parms.CJS = Zero-bias collector-substrate capacitance (F)
%   model_parms.VJS = Substrate junction built-in potential
%   model_parms.MJS = Substrate junction exponential factor
%   model_parms.XTB = Forward and reverse beta temperature exponent
%   model_parms.EG = Energy gap for temperature effect on IS (eV)
%   model_parms.XTI = Temperature exponent for effect on 1S
%   model_parms.KF = Flicker-noise coefficient
%   model_parms.AF = Flicker-noise exponent
%   model_parms.FC = Coefficient for forward-bias depletion
%
% Most of these parameters are not currently used, but they are included
% in the model specification so that future enhancements may use them.
%
% The following model is a simplified version of the SPICE model
% for the behavior of bipolar junction transistors.  The basic model
% consists of standard diode currents flowing on the base-emitter and
% base-collector junctions, and current sources which transmit a
% fraction alpha of that current through the base to the opposite 
% side, either emitter or collector.  The current on the base-emitter
% junction is called "forward" because that junction is normally
% forward-biased, whereas the current on the base-collector junction
% is called "reverse".  The model takes advantage of constraints in
% BJT transistor construction to designate a single saturation current
% for both junctions, whether in reverse or forward.
%
% The SPICE formulas are taken from the reference "Fundamentals of
% Industrial Electronics", Robert W. Erickson, Dragan Maksimovic,
% CRC Press, 2011 (Industrial Electronics Handbook, Second Edition).
% http://www.eng.auburn.edu/~wilambm/pap/2011/K10147_C009.pdf
% from which I took page 9-14 and kept a copy with this source.
%
% The model covers both PNP and NPN transistors, chosing the PNP model
% if Ve > Vb and NPN if Ve <= Vb.
%

% physical constant kT/e
  VT = 25.6e-3;

  Vbe = Vbec.Base - Vbec.Emitter;
  Vbc = Vbec.Base - Vbec.Collector;

  if (Vbe < 0)     % PNP model
    Vbe = -Vbe;
    Vbc = -Vbc;
  end

  Ib = (model_parms.IS/model_parms.BF)*(exp(Vbe/(model_parms.NF*VT)) - 1) + ...
       (model_parms.IS/model_parms.BR)*(exp(Vbc/(model_parms.NR*VT)) - 1) + ...
       model_parms.ISE*(exp(Vbe/(model_parms.NE*VT)) - 1) + ...
       model_parms.ISC*(exp(Vbc/(model_parms.NC*VT)) - 1);

  Qx = (model_parms.IS/model_parms.IKF)*(exp(Vbe/(model_parms.NF*VT)) - 1) + ...
       (model_parms.IS/model_parms.IKR)*(exp(Vbc/(model_parms.NR*VT)) - 1);
  Qb = ((1 + sqrt(1 + 4*Qx))/2) / ...
       (1 - Vbe/model_parms.VAR - Vbc/model_parms.VAF);
  Ic = (model_parms.IS/Qb)*(exp(Vbe/(model_parms.NF*VT)) - ...
                            exp(Vbc/(model_parms.NR*VT))) - ...
       (model_parms.IS/model_parms.BR)*(exp(Vbc/(model_parms.NR*VT)) - 1) - ...
       model_parms.ISC*(exp(Vbc/(model_parms.NC*VT)) - 1);
 
  currents.Base = Ib;
  currents.Collector = Ic;
  currents.Emitter = Ib + Ic;
