function response = stiffener(var,inp,f,Iin)
% Expresses the variables and constraints of a linear circuit model.
% Input frequency parameter f is in Hz.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file was automatically generated by mpp.pl
% based on stiffener.mpp so do not make changes here!
% Edit stiffener.mpp to make changes, and regenerate
% this file using the mpp.pl tool.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  w = 2*pi*f;

% resistor values (Ohms)
% effective resistance of the circuit is very roughly -R1*R3/(R1+R2)
  R1 = 2e3;
  R2 = 10e3;
  R3 = 72e3;
  R4 = 10e3;

% dc bias levels (V)
  VDAC = 80 + 0;
  VRef = 0.696 + 0;
  %Iin = 50e-6 + 0;

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

  % set the initial operating points for transistors
  % currents are in Amperes
  % junction potentials are in V
  % junction conductances 1/Ohms
  global sQ0Vbe0 sQ0Ie0 sB0 sQ0;
  global sQ1Vbe0 sQ1Ie0 sB1 sQ1;
  if isempty(sQ0Vbe0)
    sQ0Vbe0 = 0.66;
    sQ0Ie0 = 2.1e-3;
    sQ1Vbe0 = 0.66;
    sQ1Ie0 = 2.1e-3;
    sB0 = 100;
    sB1 = 100;
    sQ0 = Q0Ie0/VT;
    sQ1 = Q1Ie0/VT;
    warning('calling calib_stiff');
    calib_stiff;
  end

% equations
% Internal voltages are named according to the following rules:
%  1) if connected to the base/emitter/collector of transitor Qn,
%	call it VQnb/VQne/VQnc.
%  2) if located on a junction between resistors Ra and Rb,
%	call it VRaRb.
%  3) if located on a junction between resistor Ra and capacitor Cb,
%	call it VRaCb.
%  4) if located on a junction between two capacitors Ca and Cb,
%	call it VCaCb.
%  5) if part of the power distribution, give it a unique name,
%	for example, VPa, VPs, etc.
% Internal currents are named according to the following rules:
%  1) if flowing on the base/emitter/collector of transistor Qn,
%	call it Jbn/Jen/Jcn.
%  2) if flowing on a resistor or capacitor Ra/Cb,
%	call it Ia/Hb.
%  3) sign convention is + for current flowing down or right on the
%	Altium circuit schematic for the preamp.

%resistors

%transistors

%current junctions

% Zener equation
  Nvariables = 24;
% these are the constraints
  Nconstraints = 20;
  constraint = zeros(Nconstraints,Nvariables);
%   0 = VDAC-VQ0e-I1*R1
  constraint(1,1) = -1; % VQ0e
  constraint(1,15) = -R1; % I1
  constraint(1,21) = +VDAC; % constant
  constraint(1,22) = +1; % VDAC
%   0 = VDAC-VQ1e-I2*R2
  constraint(2,2) = -1; % VQ1e
  constraint(2,17) = -R2; % I2
  constraint(2,21) = +VDAC; % constant
  constraint(2,22) = +1; % VDAC
%   0 = VQ0c-VZ0r-I3*R3
  constraint(3,5) = +1; % VQ0c
  constraint(3,11) = -R3; % I3
  constraint(3,19) = -1; % VZ0r
%   0 = VZ0r-Vout-I4*R4
  constraint(4,13) = -R4; % I4
  constraint(4,19) = +1; % VZ0r
  constraint(4,20) = -1; % Vout
%   0 = Jb0+Jc0-Je0
  constraint(5,7) = +1; % Jc0
  constraint(5,8) = +1; % Jb0
  constraint(5,16) = -1; % Je0
%   0 = Jc0-sB0*Jb0
  constraint(6,7) = +1; % Jc0
  constraint(6,8) = -sB0; % Jb0
%   0 = VQ0e*sQ0-VQ0b*sQ0-sQ0Vbe0*sQ0+sQ0Ie0-Je0
  constraint(7,1) = +sQ0; % VQ0e
  constraint(7,6) = -sQ0; % VQ0b
  constraint(7,16) = -1; % Je0
  constraint(7,21) = -sQ0Vbe0*sQ0+sQ0Ie0; % constant
%   0 = Jb1-Jc1+Je1
  constraint(8,9) = +1; % Jb1
  constraint(8,10) = -1; % Jc1
  constraint(8,18) = +1; % Je1
%   0 = Jc1+sB1*Jb1
  constraint(9,9) = +sB1; % Jb1
  constraint(9,10) = +1; % Jc1
%   0 = VQ1e*sQ1-VQ1b*sQ1-sQ1Vbe0*sQ1+sQ1Ie0-Je1
  constraint(10,2) = +sQ1; % VQ1e
  constraint(10,3) = -sQ1; % VQ1b
  constraint(10,18) = -1; % Je1
  constraint(10,21) = -sQ1Vbe0*sQ1+sQ1Ie0; % constant
%   0 = VQ0b-VQ1b
  constraint(11,3) = -1; % VQ1b
  constraint(11,6) = +1; % VQ0b
%   0 = VQ1c-VZ0r
  constraint(12,4) = +1; % VQ1c
  constraint(12,19) = -1; % VZ0r
%   0 = VQ0c-VQ0b
  constraint(13,5) = +1; % VQ0c
  constraint(13,6) = -1; % VQ0b
%   0 = Iin-I1-I2
  constraint(14,15) = -1; % I1
  constraint(14,17) = -1; % I2
  constraint(14,21) = +Iin; % constant
  constraint(14,24) = +1; % Iin
%   0 = Jc0+Jb0-Jb1-I3-IZ0
  constraint(15,7) = +1; % Jc0
  constraint(15,8) = +1; % Jb0
  constraint(15,9) = -1; % Jb1
  constraint(15,11) = -1; % I3
  constraint(15,12) = -1; % IZ0
%   0 = Jc1+I3-I4
  constraint(16,10) = +1; % Jc1
  constraint(16,11) = +1; % I3
  constraint(16,13) = -1; % I4
%   0 = IZ0+I4-Iout
  constraint(17,12) = +1; % IZ0
  constraint(17,13) = +1; % I4
  constraint(17,14) = -1; % Iout
%   0 = I1-Je0
  constraint(18,15) = +1; % I1
  constraint(18,16) = -1; % Je0
%   0 = I2-Je1
  constraint(19,17) = +1; % I2
  constraint(19,18) = -1; % Je1
%   0 = VZ0r-Vout-VRef
  constraint(20,19) = +1; % VZ0r
  constraint(20,20) = -1; % Vout
  constraint(20,21) = -VRef; % constant
  constraint(20,23) = -1; % VRef

% now solve for the desired variable
% in response to the designated input signal
  ivar = index_lookup(var);
  iinp = index_lookup(inp);
  solution = constraint(:,1:Nconstraints)\constraint(:,iinp);
  if ivar == 0
    response = -solution;
  else
    response = -solution(ivar);
  end

  function ivar = index_lookup(var)
    if isnumeric(var)
      ivar = var;
    elseif strcmp(var,'VQ0e')
      ivar = 1;
    elseif strcmp(var,'VQ1e')
      ivar = 2;
    elseif strcmp(var,'VQ1b')
      ivar = 3;
    elseif strcmp(var,'VQ1c')
      ivar = 4;
    elseif strcmp(var,'VQ0c')
      ivar = 5;
    elseif strcmp(var,'VQ0b')
      ivar = 6;
    elseif strcmp(var,'Jc0')
      ivar = 7;
    elseif strcmp(var,'Jb0')
      ivar = 8;
    elseif strcmp(var,'Jb1')
      ivar = 9;
    elseif strcmp(var,'Jc1')
      ivar = 10;
    elseif strcmp(var,'I3')
      ivar = 11;
    elseif strcmp(var,'IZ0')
      ivar = 12;
    elseif strcmp(var,'I4')
      ivar = 13;
    elseif strcmp(var,'Iout')
      ivar = 14;
    elseif strcmp(var,'I1')
      ivar = 15;
    elseif strcmp(var,'Je0')
      ivar = 16;
    elseif strcmp(var,'I2')
      ivar = 17;
    elseif strcmp(var,'Je1')
      ivar = 18;
    elseif strcmp(var,'VZ0r')
      ivar = 19;
    elseif strcmp(var,'Vout')
      ivar = 20;
    elseif strcmp(var,'constant')
      ivar = 21;
    elseif strcmp(var,'VDAC')
      ivar = 22;
    elseif strcmp(var,'VRef')
      ivar = 23;
    elseif strcmp(var,'Iin')
      ivar = 24;
    end
  end
end
