function result = calibrate()

% Solves for the dc voltages and currents in the preamp model
% and iterates the transistor Vbe and Ie operating points to
% give the correct small-signal ac behavior after linearization.

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

  global Q0Vbe0 Q0Ie0 B0 Q0;
  global Q1Vbe0 Q1Ie0 B1 Q1;
  global Q2Vbe0 Q2Ie0 B2 Q2;
  global Q3Vbe0 Q3Ie0 B3 Q3;
  global Q4Vbe0 Q4Ie0 B4 Q4;
  global Q5Vbe0 Q5Ie0 B5 Q5;
  global Q6Vbe0 Q6Ie0 B6 Q6;
  global Q7Vbe0 Q7Ie0 B7 Q7;
  global Q8Vbe0 Q8Ie0 B8 Q8;
  %if isempty(Q0Vbe0)
    Q0Vbe0 = 0.66;
    Q0Ie0 = 2.1e-3;
    Q1Vbe0 = 0.74;
    Q1Ie0 = 2.1e-3;
    Q2Vbe0 = 0.64;
    Q2Ie0 = 1.0e-3;
    Q3Vbe0 = 0.64;
    Q3Ie0 = 1.0e-3;
    Q4Vbe0 = 0.64;
    Q4Ie0 = 1.0e-3;
    Q8Vbe0 = 0.73;
    Q8Ie0 = 1.5e-3;
    Q5Vbe0 = 0.73;
    Q5Ie0 = 1.8e-3;
    Q6Vbe0 = 0.63;
    Q6Ie0 = 1.0e-3;
    Q7Vbe0 = 0.68;
    Q7Ie0 = 3.6e-3;
    B0 = 100;
    B1 = 100;
    B2 = 100;
    B3 = 100;
    B4 = 100;
    B8 = 100;
    B5 = 100;
    B6 = 100;
    B7 = 100;
    Q0 = Q0Ie0/VT;
    Q1 = Q1Ie0/VT;
    Q2 = Q2Ie0/VT;
    Q3 = Q3Ie0/VT;
    Q4 = Q4Ie0/VT;
    Q8 = Q8Ie0/VT;
    Q5 = Q5Ie0/VT;
    Q6 = Q6Ie0/VT;
    Q7 = Q7Ie0/VT;
  %end

% NPN transistor: BFR92P (Infineon) ------------------
  BFR92P.IS = 121.3e-18;
  BFR92P.BF = 94.733; %30;
  BFR92P.NF = 1.0947;
  BFR92P.VAF = 30;
  BFR92P.IKF = 0.46227;
  BFR92P.ISE = 129.55e-15;
  BFR92P.NE = 1.9052;
  BFR92P.BR = 10.729;
  BFR92P.NR = 0.8983;
  BFR92P.VAR = 14.599;
  BFR92P.IKR = 0.01;
  BFR92P.ISC = 0.75557e-15;
  BFR92P.NC = 1.371;
  BFR92P.RB = 14.998;
  BFR92P.IRB = 0.01652e-3;
  BFR92P.RBM = 7.8145;
  BFR92P.RE = 0.29088;
  BFR92P.RC = 0.13793;
  BFR92P.XTB = 0;
  BFR92P.EG = 1.11;
  BFR92P.XTI = 3;
  BFR92P.CJE = 10.416e-15;
  BFR92P.VJE = 0.70618;
  BFR92P.MJE = 0.34686;
  BFR92P.TF = 26.796-12;
  BFR92P.XTF = 0.3817;
  BFR92P.VTF = 0.32861;
  BFR92P.ITF = 4.4601e-3;
  BFR92P.PTF = 0;
  BFR92P.CJC = 946.47e-15;
  BFR92P.VJC = 0.84079;
  BFR92P.MJC = 0.4085;
  BFR92P.XCJC = 0.13464;
  BFR92P.TR = 1.2744e-9;
  BFR92P.CJS = 0e-15;
  BFR92P.VJS = 0.75;
  BFR92P.MJS = 0;
  BFR92P.FC = 0.99545;

% PNP transistor: BFT92W (NXP) -----------------------
  BFT92W.IS = 437.5e-18;
  BFT92W.BF = 33.58; %20;
  BFT92W.NF = 1.009;
  BFT92W.VAF = 23.39;
  BFT92W.IKF = 99.53e-3;
  BFT92W.ISE = 87.05e-15;
  BFT92W.NE = 1.943;
  BFT92W.BR = 4.947;
  BFT92W.NR = 1.002;
  BFT92W.VAR = 3.903;
  BFT92W.IKR = 5.281e-3;
  BFT92W.ISC = 35.88e-15;
  BFT92W.NC = 1.393;
  BFT92W.RB = 5.000;
  BFT92W.IRB = 1.000;
  BFT92W.RBM = 5.000;
  BFT92W.RE = 1.000;
  BFT92W.RC = 10.00;
  BFT92W.XTB = 0.000;
  BFT92W.EG = 1.110;
  BFT92W.XTI = 3.000;
  BFT92W.CJE = 746.6e-15;
  BFT92W.VJE = 600.0e-3;
  BFT92W.MJE = 0.357;
  BFT92W.TF = 17.49e-12;
  BFT92W.XTF = 1.354;
  BFT92W.VTF = 155.6e-3;
  BFT92W.ITF = 1.000e-3;
  BFT92W.PTF = 45.00;
  BFT92W.CJC = 937.1e-3;
  BFT92W.VJC = 396.4e-3;
  BFT92W.MJC = 0.200;
  BFT92W.XCJC = 0.106;
  BFT92W.TR = 8.422e-9;
  BFT92W.CJS = 0.000e-15;
  BFT92W.VJS = 750.0e-3;
  BFT92W.MJS = 0.000;
  BFT92W.FC = 0.768;

% solve for the dc operating point
  omegaDC = 1; % make it "close enough" to dc
  Q0Vbec.Base = equations('VQ0b','constant',omegaDC);
  Q0Vbec.Emitter = equations('VQ0e','constant',omegaDC);
  Q0Vbec.Collector = equations('VQ0c','constant',omegaDC);
  Q1Vbec.Base = equations('VQ1b','constant',omegaDC);
  Q1Vbec.Emitter = equations('VQ1e','constant',omegaDC);
  Q1Vbec.Collector = equations('VQ1c','constant',omegaDC);
  Q2Vbec.Base = equations('VQ2b','constant',omegaDC);
  Q2Vbec.Emitter = equations('VQ2e','constant',omegaDC);
  Q2Vbec.Collector = equations('VQ2c','constant',omegaDC);
  Q3Vbec.Base = equations('VQ3b','constant',omegaDC);
  Q3Vbec.Emitter = equations('VQ3e','constant',omegaDC);
  Q3Vbec.Collector = equations('VQ3c','constant',omegaDC);
  Q4Vbec.Base = equations('VQ4b','constant',omegaDC);
  Q4Vbec.Emitter = equations('VQ4e','constant',omegaDC);
  Q4Vbec.Collector = equations('VQ4c','constant',omegaDC);
  Q5Vbec.Base = equations('VQ5b','constant',omegaDC);
  Q5Vbec.Emitter = equations('VQ5e','constant',omegaDC);
  Q5Vbec.Collector = equations('VQ5c','constant',omegaDC);
  Q6Vbec.Base = equations('VQ6b','constant',omegaDC);
  Q6Vbec.Emitter = equations('VQ6e','constant',omegaDC);
  Q6Vbec.Collector = equations('VQ6c','constant',omegaDC);
  Q7Vbec.Base = equations('VQ7b','constant',omegaDC);
  Q7Vbec.Emitter = equations('VQ7e','constant',omegaDC);
  Q7Vbec.Collector = equations('VQ7c','constant',omegaDC);
  Q8Vbec.Base = equations('VQ8b','constant',omegaDC);
  Q8Vbec.Emitter = equations('VQ8e','constant',omegaDC);
  Q8Vbec.Collector = equations('VQ8c','constant',omegaDC);
  newQ0Vbe0 = Q0Vbec.Emitter - Q0Vbec.Base;
  newQ1Vbe0 = Q1Vbec.Base - Q1Vbec.Emitter;
  newQ2Vbe0 = Q2Vbec.Emitter - Q2Vbec.Base;
  newQ3Vbe0 = Q3Vbec.Emitter - Q2Vbec.Base;
  newQ4Vbe0 = Q4Vbec.Emitter - Q4Vbec.Base;
  newQ5Vbe0 = Q5Vbec.Base - Q5Vbec.Emitter;
  newQ6Vbe0 = Q6Vbec.Emitter - Q6Vbec.Base;
  newQ7Vbe0 = Q7Vbec.Emitter - Q7Vbec.Base;
  newQ8Vbe0 = Q8Vbec.Base - Q8Vbec.Emitter;
  newQ0Ie0 = equations('Je0','constant',omegaDC);
  newQ1Ie0 = equations('Je1','constant',omegaDC);
  newQ2Ie0 = equations('Je2','constant',omegaDC);
  newQ3Ie0 = equations('Je3','constant',omegaDC);
  newQ4Ie0 = equations('Je4','constant',omegaDC);
  newQ5Ie0 = equations('Je5','constant',omegaDC);
  newQ6Ie0 = equations('Je6','constant',omegaDC);
  newQ7Ie0 = equations('Je7','constant',omegaDC);
  newQ8Ie0 = equations('Je8','constant',omegaDC);

  resolution = 1e-6;
  while abs(newQ0Vbe0-Q0Vbe0)>resolution || ...
        abs(newQ1Vbe0-Q1Vbe0)>resolution || ...
        abs(newQ2Vbe0-Q2Vbe0)>resolution || ...
        abs(newQ3Vbe0-Q3Vbe0)>resolution || ...
        abs(newQ4Vbe0-Q4Vbe0)>resolution || ...
        abs(newQ5Vbe0-Q5Vbe0)>resolution || ...
        abs(newQ6Vbe0-Q6Vbe0)>resolution || ...
        abs(newQ7Vbe0-Q7Vbe0)>resolution || ...
        abs(newQ8Vbe0-Q8Vbe0)>resolution || ...
        abs(newQ0Ie0-Q0Ie0)>resolution || ...
        abs(newQ1Ie0-Q1Ie0)>resolution || ...
        abs(newQ2Ie0-Q2Ie0)>resolution || ...
        abs(newQ3Ie0-Q3Ie0)>resolution || ...
        abs(newQ4Ie0-Q4Ie0)>resolution || ...
        abs(newQ5Ie0-Q5Ie0)>resolution || ...
        abs(newQ6Ie0-Q6Ie0)>resolution || ...
        abs(newQ7Ie0-Q7Ie0)>resolution || ...
        abs(newQ8Ie0-Q8Ie0)>resolution
    Q0Vbe0 = (Q0Vbe0+newQ0Vbe0)/2;
    currents = transistor_model(Q0Vbec,BFT92W);
    Q0Ie0 = currents.Emitter;
    B0 = currents.Collector/(currents.Base+1e-30);
    Q0 = Q0Ie0/VT;
    Q1Vbe0 = (Q1Vbe0+newQ1Vbe0)/2;
    currents = transistor_model(Q1Vbec,BFR92P);
    Q1Ie0 = currents.Emitter;
    B1 = currents.Collector/(currents.Base+1e-30);
    Q1 = Q1Ie0/VT;
    Q2Vbe0 = (Q2Vbe0+newQ2Vbe0)/2;
    currents = transistor_model(Q2Vbec,BFT92W);
    Q2Ie0 = currents.Emitter;
    B2 = currents.Collector/(currents.Base+1e-30);
    Q2 = Q2Ie0/VT;
    Q3Vbe0 = (Q3Vbe0+newQ3Vbe0)/2;
    currents = transistor_model(Q3Vbec,BFT92W);
    Q3Ie0 = currents.Emitter;
    B3 = currents.Collector/(currents.Base+1e-30);
    Q3 = Q3Ie0/VT;
    Q4Vbe0 = (Q4Vbe0+newQ4Vbe0)/2;
    currents = transistor_model(Q4Vbec,BFT92W);
    Q4Ie0 = currents.Emitter;
    B4 = currents.Collector/(currents.Base+1e-30);
    Q4 = Q4Ie0/VT;
    Q8Vbe0 = (Q8Vbe0+newQ8Vbe0)/2;
    currents = transistor_model(Q8Vbec,BFR92P);
    Q8Ie0 = currents.Emitter;
    B8 = currents.Collector/(currents.Base+1e-30);
    Q8 = Q8Ie0/VT;
    Q5Vbe0 = (Q5Vbe0+newQ5Vbe0)/2;
    currents = transistor_model(Q5Vbec,BFR92P);
    Q5Ie0 = currents.Emitter;
    B5 = currents.Collector/(currents.Base+1e-30);
    Q5 = Q5Ie0/VT;
    Q6Vbe0 = (Q6Vbe0+newQ6Vbe0)/2;
    currents = transistor_model(Q6Vbec,BFT92W);
    Q6Ie0 = currents.Emitter;
    B6 = currents.Collector/(currents.Base+1e-30);
    Q6 = Q6Ie0/VT;
    Q7Vbe0 = (Q7Vbe0+newQ7Vbe0)/2;
    currents = transistor_model(Q7Vbec,BFT92W);
    Q7Ie0 = currents.Emitter;
    B7 = currents.Collector/(currents.Base+1e-30);
    Q7 = Q7Ie0/VT;

    Q0Vbec.Base = equations('VQ0b','constant',omegaDC);
    Q0Vbec.Emitter = equations('VQ0e','constant',omegaDC);
    Q0Vbec.Collector = equations('VQ0c','constant',omegaDC);
    Q1Vbec.Base = equations('VQ1b','constant',omegaDC);
    Q1Vbec.Emitter = equations('VQ1e','constant',omegaDC);
    Q1Vbec.Collector = equations('VQ1c','constant',omegaDC);
    Q2Vbec.Base = equations('VQ2b','constant',omegaDC);
    Q2Vbec.Emitter = equations('VQ2e','constant',omegaDC);
    Q2Vbec.Collector = equations('VQ2c','constant',omegaDC);
    Q3Vbec.Base = equations('VQ3b','constant',omegaDC);
    Q3Vbec.Emitter = equations('VQ3e','constant',omegaDC);
    Q3Vbec.Collector = equations('VQ3c','constant',omegaDC);
    Q4Vbec.Base = equations('VQ4b','constant',omegaDC);
    Q4Vbec.Emitter = equations('VQ4e','constant',omegaDC);
    Q4Vbec.Collector = equations('VQ4c','constant',omegaDC);
    Q5Vbec.Base = equations('VQ5b','constant',omegaDC);
    Q5Vbec.Emitter = equations('VQ5e','constant',omegaDC);
    Q5Vbec.Collector = equations('VQ5c','constant',omegaDC);
    Q6Vbec.Base = equations('VQ6b','constant',omegaDC);
    Q6Vbec.Emitter = equations('VQ6e','constant',omegaDC);
    Q6Vbec.Collector = equations('VQ6c','constant',omegaDC);
    Q7Vbec.Base = equations('VQ7b','constant',omegaDC);
    Q7Vbec.Emitter = equations('VQ7e','constant',omegaDC);
    Q7Vbec.Collector = equations('VQ7c','constant',omegaDC);
    Q8Vbec.Base = equations('VQ8b','constant',omegaDC);
    Q8Vbec.Emitter = equations('VQ8e','constant',omegaDC);
    Q8Vbec.Collector = equations('VQ8c','constant',omegaDC);
    newQ0Vbe0 = Q0Vbec.Emitter - Q0Vbec.Base;
    newQ1Vbe0 = Q1Vbec.Base - Q1Vbec.Emitter;
    newQ2Vbe0 = Q2Vbec.Emitter - Q2Vbec.Base;
    newQ3Vbe0 = Q3Vbec.Emitter - Q2Vbec.Base;
    newQ4Vbe0 = Q4Vbec.Emitter - Q4Vbec.Base;
    newQ5Vbe0 = Q5Vbec.Base - Q5Vbec.Emitter;
    newQ6Vbe0 = Q6Vbec.Emitter - Q6Vbec.Base;
    newQ7Vbe0 = Q7Vbec.Emitter - Q7Vbec.Base;
    newQ8Vbe0 = Q8Vbec.Base - Q8Vbec.Emitter;
    newQ0Ie0 = equations('Je0','constant',omegaDC);
    newQ1Ie0 = equations('Je1','constant',omegaDC);
    newQ2Ie0 = equations('Je2','constant',omegaDC);
    newQ3Ie0 = equations('Je3','constant',omegaDC);
    newQ4Ie0 = equations('Je4','constant',omegaDC);
    newQ5Ie0 = equations('Je5','constant',omegaDC);
    newQ6Ie0 = equations('Je6','constant',omegaDC);
    newQ7Ie0 = equations('Je7','constant',omegaDC);
    newQ8Ie0 = equations('Je8','constant',omegaDC);
  end
