function [inputFilterModules,outputFilterModules,osemModel,coilDriverModel,acBandLimFilterModules,dcBandLimFilterModules] = make_OSEM_filter_model(freq,susType,isoStage,varargin)
% [inputFilterModules,outputFilterModules,osemModel,coilDriverModel,acBandLimFilterModules] = make_OSEM_filter_model(freq,susType,isoStage,varargin)
%
% Model of the OSEM Flag/LED Sensor and OSEM Coil/Magnet Actuator
%
% INPUT:
% [[vv mandatory vv]]
% freq - desired frequency vector
% susType - four character string indicating type of aLIGO suspension, i.e.
%           'QUAD','BSFM','HLTS','HSTS','OMCS','TMTS','HTTS','HAUX',
%           'HXDS','HSDS','HPDS','HTDS','HDDS','OFIS','OPOS'
% isoType - 2 or 5 character string indicating the isolation stage of an aLIGO
%           suspension, i.e.
%           'M0','R0','L1','L2','L3','M1','M2','M3',''
%           ('' is reserved for single suspensions such as the HTTS and HAUX).
%           Note, for drivers that have been permanently modified, prefix
%           with 'MOD' for correct response; possibilities include
%           'MODM2' for HSTS, HLTS, and BSFM (TACQ drivers modified as per L1200226)
%           'MODM3' for HSTS, and HLTS (TACQ drivers modified as per L1200226)
%           'MODL1' for QUAD (UIM driver modified as per T1400223)
%           'MODM1' for HXDS, HTTS, HAUX, OFIS, or OPOS (HAM-A driver modified as per T2100410)
% [[vv optional vv]]
% coilDriverState - Defining the state for a given coil driver. Options:
%   TOP (triple or quad) - {'acq','lp'}
%   UIM - {'acq','lp1','lp2','lp3'}
%   PUM - {'acqOfflpOff','acqOnlpOff','acqOfflpOn','acqOnLpOn'}
%   TAQ - {'acqOfflpOff','acqOnlpOff','acqOfflpOn','acqOnLpOn'}
%   HMA - {'acq','lp'}
% plotFigs - boolean to trigger plotting of coil driver response
%
% OUTPUT:
% [[Note, currently digital filters are hard-coded to be discretized at
% 16384 Hz]] 
%
% inputFilterModules = (SatAmp compensation filters)
% 1x10 struct array with fields:
%     name        : Name of Filter Module name for EPICs
%     zeros_Hz    : list of zeros in the filter [Hz]
%     poles_Hz    : list of poles in the filter [Hz]
%     zeros_radps   : list of zeros in the filter [rad] (-2*pi*zeros_Hz)
%     poles_radps   : list of poles in the filter [rad] (-2*pi*zeros_Hz)
%     norm        : gain which normalizes the frequency response 
%     c           : continuous zpk object of the filter 
%     d           : discrete zpk object of the filter (to be used for FOTON)
%     f           : frequency repsonse of filter, squeezeed over freq
%     used        : boolean array indicating which filter modules are used
%
% outputFilterModules = (CoilDriver compensation filters)
% 1x10 struct array with fields: (same as inputFilterModules)
%
% osemModel = 
%     ctsperm: (sensor calibration in counts per meter [assume satamp freq resp is compensated])
%               c: continuous zpk object of the filter 
%               f: frequency repsonse of filter, squeezeed over freq
%     ctspernm: (sensor calibration in counts per nanometer [no satamp freq resp])
%               c: continuous zpk object of the filter 
%               f: frequency repsonse of filter, squeezeed over freq
%        NperA: Force Coefficient for coil magnet pair [N/A]
%       Nperct: (actuator calibration in Newtons per count [assume coil driver response is compensated])
%               c: continuous zpk object of the filter 
%               f: frequency repsonse of filter, squeezeed over freq
%
% coilDriverModel = 
%         type  : name of coil driver type
%        modes  : list of possible driver states
%           Rc  : coil resistance (from coil driver to OSEM coil and back)
%       dcGain  : gain of driver in [A/V]
%     zeros_Hz  : list of zeros in the filter [Hz]
%     poles_Hz  : list of poles in the filter [Hz]
%    zeros_radps  : list of zeros in the filter [rad] (-2*pi*zeros_Hz)
%    poles_radps  : list of poles in the filter [rad] (-2*pi*zeros_Hz)
%         norm  : gain which normalizes the frequency response 
%            c  : [1x1 zpk]
%            d  : [1x1 zpk]
%            f  : [500x1 double]
%       fmName: {' '  'acq'}


%% Generic Parameters
sampleFreq = 16384;
global sampleTime 
global radFreq 
global plotFigs
sampleTime = 1/sampleFreq;
radFreq = 2*pi*freq;

global digitizationType
digitizationType = 'tustin';

%% Variable inputs
%Coildriver type and plot flags
if length(varargin) < 2
    plotFigs = false;
else
    plotFigs = varargin{2};
end

if length(varargin) < 1
    coilDriverState = 'None';
else
    coilDriverState = varargin{1};
end


%% %% FILTER MATLAB STRUCTURE %% %%
global emptyFilter;
emptyFilter.name = '';
emptyFilter.zeros_Hz = [];
emptyFilter.poles_Hz = [];
emptyFilter.zeros_radps = [];
emptyFilter.poles_radps = [];
emptyFilter.norm = [];
emptyFilter.c = [];
emptyFilter.d = [];
emptyFilter.f = [];
emptyFilter.used = false;


%% Generic Electronics %%
% General Standards, 32 Channel, 16-bit ADC
adc.gain = (2^16) / 40; % [cts / V]

% General Standards, 8 Channel, 18-bit DAC
dac.gain = 20 / (2^18); % [V / ct];

% aLIGO AA Filter D070081
% Frequency Response of 
% 3rd order Butterwoth Filter at 10kHz, w/ notch at 65535 Hz
% is ignored
antialiasing.gain = +1; % [V/V]

% aLIGO AI Filter D070081
antiimaging.gain = +1; % [V/V]

%% SENSOR SIDE %%%
osemType = [susType isoStage];
osem.sensitivity = getosemsensitivity(susType,isoStage); %[A/m]

satAmp = getsatamp(susType,isoStage);


% cts / meter
%osemModel.ctsperm.c = osem.sensitivity *...   % [A/m]
%                      satAmp.dcGain * ...     % [V/A]
%                      zpk([],[],1) * ...      % []
%                      antialiasing.gain * ... % [V/V]
%                      adc.gain;               % [cts/V]
                  
                  
% cts / meter

%See LLO alog 4291  %https://alog.ligo-la.caltech.edu/aLOG/index.php?callRep=4291
%Ideal osem counts per meter
osemModel.ctsperm.c = 30000 / 0.0007 * zpk([],[],1); %Ideal 30,000 counts per 0.7 mm linear range (0.0007 m)

% Convert to cts / nanometer
osemModel.ctspernm.c = 1e-9 * osemModel.ctsperm.c;  % [cts/nm]      
          
%
osemModel.ctsperm.f  = squeeze(freqresp(osemModel.ctsperm.c,2*pi*freq));   % [cts/m]
osemModel.ctspernm.f = squeeze(freqresp(osemModel.ctspernm.c,2*pi*freq));  % [cts/nm]


%% ACTUATOR SIDE %%

% Measured by Stuart Aston,
% <https://alog.ligo-la.caltech.edu/aLOG/index.php?callRep=3340>
Rc.bosem = 42.7; % [Ohms] 
Rc.aosem = 19.8; % [Ohms] 
Rc.cosem = 31.2; % [Ohms] LHO aLOG 62013
Lc.cosem = 4.93e-3; % [H] LHO aLOG 62013
Lc.bosem = 11.9e-3; % [H] T1200264 (but given recent experience with AOSEMs and COSEMs, probably more)
Lc.aosem = 4.47e-3; % [H] LHO aLOG 61726

osem.forceCoeff = getosemforcecoeff(susType,isoStage,coilDriverState); % [N/A]

[outputFilterModules, coilDriverModel] = getcoildriver(freq,susType,isoStage,coilDriverState,Rc,Lc); 

coilDriverModel = getdacbitcount(susType,isoStage,coilDriverModel);

if strcmp(coilDriverState, 'None') ~= 1
  osemModel.NperA = osem.forceCoeff;  
  
  switch osemType
      case {'HDDSM2','HCDSM2'} % The only AC coupled coil system we have thus far
          % Not really a continuous zpk, since I use 
          %      coilDriverModel.transConductance 
          % rather than
          %     zpk([],[],1)     * ...       % []
          %      coilDriverModel.dcTransConductance * ... % [A/V]
          % but I'm just trying to make the code run for now...
          osemModel.Nperct.c = dac.gain * ...             % [V/ct]
                               antiimaging.gain * ...       % [V/V]
                               coilDriverModel.transConductance * ... % [A/V]
                               osem.forceCoeff;             % [N/A]
                           
          
          osemModel.Nperct.f = osemModel.Nperct.c;
      otherwise
          % Newtons / ct
          osemModel.Nperct.c = dac.gain * ...             % [V/ct]
                               antiimaging.gain * ...       % [V/V]
                               zpk([],[],1)     * ...       % []
                               coilDriverModel.dcTransConductance * ... % [A/V]
                               osem.forceCoeff;             % [N/A]
          
         osemModel.Nperct.f = squeeze(freqresp(osemModel.Nperct.c,2*pi*freq));  
  end
                
end

switch osemType
    case {'QUADL3','BSFMM3'} % No local sensors.
        inputFilterModules(5).c = zpk([],[],1);
        inputFilterModules(5).d = zpk(c2d(ss(inputFilterModules(5).c),sampleTime,digitizationType));
        inputFilterModules(5).f = squeeze(freqresp(inputFilterModules(5).d,radFreq));
        inputFilterModules(5).name = '';
        inputFilterModules(5).used = false;
    otherwise 
        %Create input filters
        inputFilterModules = repmat(emptyFilter,1,10);
        inputFilterModules(1).c = 1/satAmp.c;
        inputFilterModules(1).d = zpk(c2d(ss(inputFilterModules(1).c),sampleTime,digitizationType));
        inputFilterModules(1).f = squeeze(freqresp(inputFilterModules(1).d,radFreq));
        inputFilterModules(1).name = '10:0.4';
        inputFilterModules(1).used = true;
        
        inputFilterModules(4).c = 1/osemModel.ctspernm.c;
        inputFilterModules(4).d = zpk(c2d(ss(inputFilterModules(4).c),sampleTime,digitizationType));
        inputFilterModules(4).f = squeeze(freqresp(inputFilterModules(4).d,radFreq));
        inputFilterModules(4).name = 'to_nm';
        inputFilterModules(4).used = true;
        
        inputFilterModules(5).c = 1/osemModel.ctspernm.c/1000.0;
        inputFilterModules(5).d = zpk(c2d(ss(inputFilterModules(5).c),sampleTime,digitizationType));
        inputFilterModules(5).f = squeeze(freqresp(inputFilterModules(5).d,radFreq));
        inputFilterModules(5).name = 'to_um';
        inputFilterModules(5).used = true;
    
end

%Create Watchdog AC band limited filters

acBandLimFilterModules = repmat(emptyFilter,1,10);
acBandLimFilterModules(1).zeros = -2*pi*[0];
acBandLimFilterModules(1).poles = -2*pi*[0.1;10;10];

acBandLimFilterModules(1).norm = 3987.4;


acBandLimFilterModules(1).c = zpk(acBandLimFilterModules(1).zeros,acBandLimFilterModules(1).poles,acBandLimFilterModules(1).norm);
acBandLimFilterModules(1).d = zpk(c2d(ss(acBandLimFilterModules(1).c),sampleTime,digitizationType));
acBandLimFilterModules(1).f = squeeze(freqresp(acBandLimFilterModules(1).d,radFreq));
acBandLimFilterModules(1).name = 'acBandLim';
acBandLimFilterModules(1).used = true;

%Create Watchdog DC band limited filters

dcBandLimFilterModules = repmat(emptyFilter,1,10);
dcBandLimFilterModules(1).zeros = -2*pi*[];
dcBandLimFilterModules(1).poles = -2*pi*[10];

dcBandLimFilterModules(1).norm = 62.8319;


dcBandLimFilterModules(1).c = zpk(dcBandLimFilterModules(1).zeros,dcBandLimFilterModules(1).poles,dcBandLimFilterModules(1).norm);
dcBandLimFilterModules(1).d = zpk(c2d(ss(dcBandLimFilterModules(1).c),sampleTime,digitizationType));
dcBandLimFilterModules(1).f = squeeze(freqresp(dcBandLimFilterModules(1).d,radFreq));
dcBandLimFilterModules(1).name = 'dcBandLim';
dcBandLimFilterModules(1).used = true;

%% Figures
if plotFigs
    figure;
    ll=loglog(freq,abs(osemModel.ctspernm.f));
    set(gca,'FontSize',16)
    title({'Model OSEM Flag/LED/PD Sensitivity';...
           [susType ' ' isoStage ', ' coilDriverState ' Mode']})
    set(ll,'LineWidth',2)
    ylabel('Magnitude (ct/nm)')
    xlabel('Frequency Hz')
    xlim([freq(1) freq(end)])
    grid on
    
    figure
    ll = loglog(freq,abs(osemModel.Nperct.f));
    set(gca,'FontSize',16)
    title({'Model Magnet/Coil Drive Strength';...
           [susType ' ' isoStage ', ' coilDriverState ' Mode']})
    set(ll,'LineWidth',2)
    ylabel('Magnitude (N/ct)')
    xlabel('Frequency (Hz)')
    xlim([freq(1) freq(end)])
    grid on
end



function sensitivity = getosemsensitivity(susType,isoStage)
% sensitivity = getosemsensitivity(susType,isoStage)
% 
% susType options: 'QUAD','BSFM','HLTS','HSTS','OMCS','TMTS','HTTS','HAUX'
% isoStage options: 'M0','R0','L1','L2','L3','M1','M2','M3',''
% ('' is reserved for single suspensions such as the HTTS and HAUX).
% 
% The sensitivity of aLIGO OSEM sensors, once normalized is 95 uA/mm
%
% The sensitivity is returned in [A/m]

sens.A = 95e-6/1e-3; % [A/m] <https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=3056>
sens.B = 95e-6/1e-3; % [A/m] ~20 V/m with a satellite amp gain of 320k Ohms


osemType = [susType isoStage];

switch osemType 
    case {'QUADM0','QUADR0','QUADL1','BSFMM1','BSFMM2','HLTSM1','HSTSM1','TMTSM1','OMCSM1','HTTSM1','HSDSM1','HPDSM1','HTDSM1','HDDSM1','HXDSM1','HCDSM1',...
          'QUADMODL1','HSTSMODM1','HTTSMODM1','HSDSMODM1','HPDSMODM1','HTDSMODM1','HDDSMODM1','HXDSMODM1','BBSSM1','BBSSM2','HRTSM1','HATSM1'}
        sensitivity = sens.B;
    case {'QUADL2','HSTSM2','HSTSM3', 'HLTSM2','HLTSM3','HAUXM1','OFISM1','OPOSM1',...
          'HSTSMODM2','HSTSMODM3','HLTSMODM2','HLTSMODM3','BSFMMODM2','HAUXMODM1','OFISMODM1','OPOSMODM1','BBSSM3','HATSM3','BHSSM1'}
        sensitivity = sens.A;
    case {'QUADL3','BSFMM3','TMTSM2','OMCSM2','HRTSM2','HATSM2','HDDSM2','HCDSM2'}
        disp(['WARNING: There are no sensors on a ' susType ' ' isoStage ' stage. Setting sensitivity to 0'])
        sensitivity = 0;
    otherwise
        error([isoStage ' is not a recognized stage of a ' susType])
end


function forceCoeff = getosemforcecoeff(susType,isoStage,cdState)
% forceCoeff = getosemforcecoeff(susType,isoStage,cdState)
%
% There are 6 types of Coil/Magnet Pairs for aLIGO Suspensions. This
% function returns the force coefficients (in [N/A]) of a given isoStage
% for a given susType.
% susType options: 'QUAD','BSFM','HLTS','HSTS','OMCS','TMTS','HTTS','HAUX'
% isoStage options: 'M0','R0','L1','L2','L3','M1','M2','M3'
% ('M1' should be used for single suspensions such as the HTTS and HAUX).
% 
% Force coefficients (for perfectly aligned, sweet spot positions)
% are modelled in T1000164.
%
% For QUADL3, getosemforcecoeff returns a (linear) force coefficient for
% the ESD, assuming a bias voltage of 100 [V]. 

%      Coeff [N/A]      Diameter [mm]   Thickness [mm]    Coil Type      Uses   
fmax(1) = 1.694;   %       10.0              10.0           BOSEM        QUAD M0, QUAD R0, QUAD L1, TMTS M0, BSFM M1, HLTS M1, OMCS M1, BBSS M1
fmax(2) = 0.963;   %       10.0               5.0           BOSEM        BSFM M2, HSTS M1, HXDS M1, BBSS M2
fmax(12) = 0.39254;%        4.0               4.0           BOSEM        HRTS M1, HATS M1
fmax(4) = 0.021;   %        2.0               3.0           BOSEM        HTTS M1
fmax(7) = 0.0695;  %        3.0               6.0           AOSEM        OFIS M1
fmax(3) = 0.0309;  %        2.0               6.0           AOSEM        QUAD L2, OPOS M1, HDDS M2, BHSS M1
fmax(5) = 0.0158;  %        1.9               3.2           AOSEM        HSTS M2, HLTS M2, HAUX
fmax(6) = 0.00281; %        2.0               0.5           AOSEM        HSTS M3, HLTS M3, HATS M3, BBSS M3

% ESD  Coeff [N/V] (Assumes bias voltage of 400 [V])

% These [N/V^2] values are not for the production patterns, but "good
% enough," as confirmed by M. Evans and J. Miller (in "Force Calculations from
% "Production" ERMs and CPs" email exchange to JSK June 13/14 2011) 
% Also, the cited [N/V^2] coefficient is for all L motion, i.e. four
% quadrants at once, with the same voltage. We want to eventually do P and
% Y, so I divide by four to get the per-quadrant number (with the intent to
% use it assuming that the quadrant is a point actuator, with the acutation
% point at the "center" (1/2*(OD-ID) and at 45, 135, 225, and 315 degs) of
% the quadrant's pattern). 

% To linearize, we assume the maximum possible force by taking
% maxForce_N = [N/V^2] * (dacMax_V * esdDriverGain_VpV * (quadrant - bias))^2 / nQuadrants
%            =         * (10 [V]   * 40 [V/V]          *         2        )^2 / 4
%
% (Note, the example is for the aquisition driver, with esdDriverGain_VpV = 40 [V/V]. 
% T0900567 proposes a low noise driver with gain of 1.5 [V/V].
%
% which we then distribute over 1/2 the dacMax_V, assuming we'll
% operate with the quadrants biased to 1/2 the range so we can get both
% attractive and repulsive force
% linearResponse = maxForce_N / (2*dacMax_V)
%
% (I know it looks like I'm multiplying the dacMax_V by 2, but you'll be
% multiplying this [N/V] by the dacMax_V to get the full range, and you 
% want to *reduce* the [N] you get by 2, so, you *divide* the response by
% 2)  

linearizingFactor.acq  = (10 * 40  * 2)^2 /(4*2*10); % [V] acquisition driver, T1000222
linearizingFactor.ln   = (10 * 1.1 * 2)^2 /(4*2*10); % [V] low noise driver, T0900567, T1200479
alpha.ITM = 7.5e-12; % [N/V^2] (pg 4 of T1000119)
alpha.ETM = 4.2e-10; % [N/V^2] (similar to pg 7 of G0900956)

        % [N/V^2] * [V]
fmax(8) = alpha.ITM * linearizingFactor.acq; % [N/V] 20 [mm] gap, aquisition driver
fmax(9) = alpha.ITM * linearizingFactor.ln;  % [N/V] 20 [mm] gap, low noise driver 
fmax(10) = alpha.ETM * linearizingFactor.acq; % [N/V] 5 [mm] gap, acquisition driver
fmax(11)= alpha.ETM * linearizingFactor.ln;  % [N/V] 5 [mm] gap, low noise driver 

osemType = [susType isoStage];

switch osemType 
    case {'QUADM0', 'QUADR0', 'QUADL1','BSFMM1', 'HLTSM1',...
            'TMTSM1','OMCSM1','QUADMODL1','BBSSM1'}
        forceCoeff = fmax(1);
    case {'BSFMM2', 'HSTSM1','HSTSMODM1','BSFMMODM2',...
          'HXDSM1','HSDSM1','HPDSM1','HTDSM1','HDDSM1','HCDSM1',...
          'HXDSMODM1','HSDSMODM1','HPDSMODM1','HTDSMODM1',...
          'HDDSMODM1', 'BBSSM2'}
        forceCoeff = fmax(2);
    case {'QUADL2','OPOSM1','OPOSMODM1','HDDSM2','HCDSM2','BHSSM1'}
        forceCoeff = fmax(3);
    case {'HTTSM1','HTTSMODM1'}
        forceCoeff = fmax(4);
    case {'HSTSM2', 'HLTSM2', 'HAUXM1','HSTSMODM2','HLTSMODM2'}
        forceCoeff = fmax(5);
    case {'HSTSM3','HLTSM3','HSTSMODM3','HLTSMODM3','HATSM3','BBSSM3'}
        forceCoeff = fmax(6);
    case {'OFISM1','OFISMODM1'}
        forceCoeff = fmax(7);
    case {'QUADL3'}
        switch cdState
            case 'acq20mm'
                forceCoeff = fmax(8);
            case 'ln20mm'
                forceCoeff = fmax(9);
            case 'acq5mm'
                forceCoeff = fmax(10);
            case 'ln5mm'
                forceCoeff = fmax(11);
        end
    case {'HRTSM1','HATSM1'}
        forceCoeff = fmax(12);
    otherwise
        error([isoStage ' is not a recognized stage of a ' susType])
end

function [filterModules,coilDriver] = getcoildriver(freq,susType,isoStage,coilDriverState,Rc,Lc)

global emptyFilter;
global radFreq;
global sampleTime;
global digitizationType;
osemType = [susType isoStage];

switch osemType
    case {'QUADM0', 'QUADR0','BSFMM1','HLTSM1','HSTSM1','TMTSM1','OMCSM1','HSTSMODM1','BBSSM1'}
        
       % TOP DRIVER
       % Circuit Schematic    D0902747
       % Design Description   T0900232
       % User's Guide         T0900251
       % Board Test Plan      T0900231
       % Chassis Test Plan    T0900230 
       
       coilDriver.modes      = {'acq';... % = Zsw open;
                                'lp'};    % = Zsw closed;
       
       %                                                                               {gnd-|Inf|-}+--|R29|--+
       %                                                                                           |         |
       %                                                                                           +-|N      |
       % o --|R7|--+-----------------------|P           dac                                             |    >--+--- o --+--|R90|---|R91|--+----+
       % |         |                       |    >--+--|R8 |--+------|Zrc2|--+------------------------|P           |   |                 |    |
       % |         |                     +-|N      |         |              |     gnd                             |   +--|R5 |-+-|R1 |--+    |
       % |        _|_                    |         |        _|_            _|_    _|_                             |            |            _|_
       % |       |   |         +--|R10|--+--|R3 |--+       |   |          |   |  |   |                            |            |           |   | |
       % Vin     |Zsw|    gnd--|                           |Zrc|          |Zop|  |Zop|                           Vout        C200          |Zc | | Ic
       % |       |___|         +--|R15|--+--|R20|--+       |_1_|          |___|  |___|                            |            |           |___| V
       % |         |                     |         |         |              |      |                              |            |             |
       % |         |                     +-|P      |         |             gnd     |                              |   +--|R5 |-+-|R1 |--+    |
       % |         |                       |    >--+--|R26|--+------|Zrc2|---------+-----------------|P           |   |                 |    |
       % o --|R25|-+-----------------------|N                                                        |    >--+--- o --+--|R90|---|R91|--+----+
       %                                                                                           +-|N      |
       %                                                                                           |         |
       %                                                                               {gnd-|Inf|-}+--|R30|--+
       %
       %  |--no DC gain--|------------ Gb1 -----------|-no DC gain-|----no DC gain----|--------- Gout ---------| |-------- Zout --------|
       %  |----------------------------------------------- Gtot -----------------------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R3  = 2.2e3;                         % [V/A] ( = R20)
       R10 = 10e3;                          % [V/A] ( = R15)
       Gb1 = (1 + R3/R10);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.22
       
       R29  = 100;                          % [V/A] ( = R30)
       Gout = (1 + R29/inf);                % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       coilDriver.dcVoltageGain = Gb1 * Gout;        % [Vin/Vout] 
       
       % Calculate DC output impedance, Zout
       % Yeah, that's right. I'm not fucking kidding. This is only difference between the boards.
       switch osemType
           case {'QUADM0', 'QUADR0','TMTSM1'} % QUAD TOP Driver
               coilDriver.type       = 'QUAD TOP';
               R5  = 22;      % [V/A]
               R1  = 22;      % [V/A]
               R90 = 220;     % [V/A]
               R91 = 220;     % [V/A]
           case {'BSFMM1','HLTSM1','HSTSM1','HSTSMODM1','BBSSM1'} % Triple TOP Driver
               coilDriver.type       = 'Triple TOP';
               R5  = 10;      % [V/A] 
               R1  = 22;      % [V/A]
               R90 = 220;     % [V/A]
               R91 = 220;     % [V/A]
           case {'OMCSM1'}
               coilDriver.type       = 'OMC TOP';
               R5  = 22;      % [V/A]
               R1  = 270;     % [V/A]
               R90 = 3.3e3;   % [V/A]
               R91 = 3.3e3;   % [V/A]
       end
       
       Rs1 = R5+R1;   % [V/A] Equivalent resistence along upper path
       Rs2 = R90+R91; % [V/A] Equivalent resistence along lower parallel path

       coilDriver.dcZout = Rs2*Rs1 / (Rs2 + Rs1); % [V/A] Equivalent output resistance of two parallel paths
                                                  %       C200 is ignored it only add a pole at many MHz, with 
                                                  %       no DC gain
       
       % Calculate Total DC Current Gain (Ic / Vin), in [A/V] 
       coilDriver.dcZc            = Rc.bosem;      % [V/A] At DC, we can ignore 
                                                   % the inductance, so Zc = Rc 
       
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
       
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6 
       zeros = 0.9;
       poles = 31;
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 
       zeros = 10;
       poles = 1;
       filterModules(2) = fillFilter(zeros,poles,['SimLP' isoStage]);

       %FM module 7
       zeros = 1;
       poles = 10;
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);

       % Sets total coil
       switch coilDriverState
           case coilDriver.modes(1)
               % ACQUIRE Mode (Relays NOT Energized, B I/O Switched to 0
               % (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(1)];
           case coilDriver.modes(2)
               % RUN Mode (Relays Energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);              
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(2)];
           case 'None'
               
           otherwise
               error(['Unrecognized coilDriverState "' coilDriverState '" for osemType "' osemType '"'])
       end

    case {'QUADL1'}
        
       % UIM Driver
       % Circuit Schematic D070481
       % Design Description T0900233
       % User's Guide T0900249
       % Board Test Plan T0900237
       % Chassis Test Plan  T0900236
       
       coilDriver.type  = 'UIM';
       coilDriver.modes = {'acq';'lp1';'lp2';'lp3'};
       
       % Calculate total DC Voltage Gain, Gtot in [V/V] 
       % (Differential Volts across just before output impedance /
       % Differential Volts across input)          
       %                                                     {gnd-|Inf|-}+--|R102|--+
       %                                                                 |          |
       %                                                                 +-|P       |                                                       +--|C12|--|R4 |--+
       % o --|R8 |--+-----------------|P                                   |    >---+--|R6 |--+--------------------|P                       |                |
       % |          |                 |    >---+--|R7 |--+-----------------|N                 |                    |    >--|R104|--+--- o --+------|R5 |-----+----+
       % |          |               +-|N       |         |                                    |                  +-|N              |    |                         |
       % |         _|_              |          |        _|_                                  _|_                 |                 |    |                        _|_
       % |        |   |             +--|R100|--+       |   |                                |   |      +--|R10|--+--|R3 |----------+    |                       |   | |
       % Vin      |Zsw| {gnd-|Inf|--|}                 |Zsw|                                |Zsw|  gnd-|                               Vout                     |Zc | |Ic
       % |        |_1_|             +--|R101|--+       |_2_|                                |_3_|      +--|R15|--+--|R20|----------+    |                       |___| V
       % |          |               |          |         |                                    |                  |                 |    |                         |
       % |          |               +-|P       |         |                                    |                  +-|N              |    |                         |
       % |          |                 |    >---+--|R25|--+-----------------|P                 |                    |    >--|R105|--+--- o --+------|R23|-----+----+
       % o --|R26|--+-----------------|N                                   |    >---+--|R24|--+--------------------|P                       |                |
       %                                                                 +-|N       |                                                       +--|C26|--|R27|--+                  
       %                                                                 |          |               
       %                                                     {gnd-|Inf|-}+--|R103|--+
       %
       %  |-no DC gain-|--------- Gb1 ----------|-no DC gain-|--------- Gb2 ----------|-no DC gain-|--------------- Gout ---------------|   |------ Zout ------|
       %  |-------------------------------------------------------- Gtot ---------------------------------------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot = Vin / Vout
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       R100 = 510;              % [V/A] ( = R101)
       Gb1  = (1 + R100/inf);   % [V/V] Unity Gain "Follower" Buffer Amp. Nominally 1.00
       
       R102 = 510;              % [V/A] ( = R103)
       Gb2  = (1 + R102/inf);   % [V/V] Unity Gain "Follower" Buffer Amp. Nominally 1.00
       
       R3  = 3.3e3;             % [V/A] ( = R20)
       R10 = 10e3;              % [V/A] ( = R15)
       Gout = (1 + R3/R10);     % [V/V] Non-inverting Amp. Gain. Nominally 1.33
       
       coilDriver.dcVoltageGain = Gb1 * Gb2 * Gout; % [Vin/Vout] 
       
       % Calculate DC output impedance, Zout
       R104 = 10;
       % R5 = 3.9e3;
       R5 = 4.3e3; % 10% increase, to better match transconductance data (see LLO aLOG 4495)
       
       
       coilDriver.dcZout = R104 + R5;           % [V/A] Equivalent output resistance of two parallel paths
                                                %       Note that R4 + C12 path is not included, because its 
                                                %       a low pass filter with freq.resp [z:p] = [300:50] 
                                                %       and unity gain.
       
       % Calculate total Current Gain (Ic / Vin), in [A/V]
       coilDriver.dcZc       = Rc.bosem;        % [V/A] At DC, we can ignore 
                                                % the inductance, so Zc = Rc 
                                                
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
        
        
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = 300;
       poles = 50;
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 
       zeros = 10.5;
       poles = 1;
       filterModules(2) = fillFilter(zeros,poles,['SimLP1' isoStage]);

       %FM module 7
       zeros = 1;
       poles = 10.5;
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP1' isoStage]);
        
       
       %Fill in FM module 3 
       zeros = 10.5;
       poles = 1;
       filterModules(3) = fillFilter(zeros,poles,['SimLP2' isoStage]);

       %FM module 8
       zeros = 1;
       poles = 10.5;
       filterModules(8)  = fillFilter(zeros,poles,['AntiLP2' isoStage]);
        
       %Fill in FM module 4 
       zeros = 10.5;
       poles = 1;
       filterModules(4) = fillFilter(zeros,poles,['SimLP2' isoStage]);

       %FM module 9
       zeros = 1;
       poles = 10.5;
       filterModules(9)  = fillFilter(zeros,poles,['AntiLP2' isoStage]);
       
        switch coilDriverState
            case 'acq' % DEFAULT Mode (NO relays energized, all B I/O Switched to 0 (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];

            case 'lp1' % LP1 ON (LP1 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
          
               
            case 'lp2' % LP2 ON (LP1 and LP2 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
               
                
            case 'lp3' % LP3 ON (LP1, LP2, and LP3 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7 8 9];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
               
            case 'None'
               
        end
        
    case {'QUADMODL1'}
        
       % Modified UIM Driver
       % ECR E1400164
       % Circuit Schematic D070481 + T1400223
       % Design Description T1400223
       % User's Guide ????
       % Board Test Plan ????
       % Chassis Test Plan  ????
       
       coilDriver.type  = 'UIM';
       coilDriver.modes = {'acq';'lp1';'lp2';'lp3'};
       
       % Calculate total DC Voltage Gain, Gtot in [V/V] 
       % (Differential Volts across just before output impedance /
       % Differential Volts across input)          
       %                                                     {gnd-|Inf|-}+--|R102|--+                         ECR: Change R10 & R15 to 2.2k, Change R5&R23 to 2k 
       %                                                                 |          |
       %                                                                 +-|P       |                                                       +--|C12|--|R4 |--+
       % o --|R8 |--+-----------------|P                                   |    >---+--|R6 |--+--------------------|P                       |                |
       % |          |                 |    >---+--|R7 |--+-----------------|N                 |                    |    >--|R104|--+--- o --+------|R5 |-----+----+
       % |          |               +-|N       |         |                                    |                  +-|N              |    |                         |
       % |         _|_              |          |        _|_                                  _|_                 |                 |    |                        _|_
       % |        |   |             +--|R100|--+       |   |                                |   |      +--|R10|--+--|R3 |----------+    |                       |   | |
       % Vin      |Zsw| {gnd-|Inf|--|}                 |Zsw|                                |Zsw|  gnd-|                               Vout                     |Zc | |Ic
       % |        |_1_|             +--|R101|--+       |_2_|                                |_3_|      +--|R15|--+--|R20|----------+    |                       |___| V
       % |          |               |          |         |                                    |                  |                 |    |                         |
       % |          |               +-|P       |         |                                    |                  +-|N              |    |                         |
       % |          |                 |    >---+--|R25|--+-----------------|P                 |                    |    >--|R105|--+--- o --+------|R23|-----+----+
       % o --|R26|--+-----------------|N                                   |    >---+--|R24|--+--------------------|P                       |                |
       %                                                                 +-|N       |                                                       +--|C26|--|R27|--+                  
       %                                                                 |          |               
       %                                                     {gnd-|Inf|-}+--|R103|--+
       %
       %  |-no DC gain-|--------- Gb1 ----------|-no DC gain-|--------- Gb2 ----------|-no DC gain-|--------------- Gout ---------------|   |------ Zout ------|
       %  |-------------------------------------------------------- Gtot ---------------------------------------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot = Vin / Vout
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       R100 = 510;              % [V/A] ( = R101)
       Gb1  = (1 + R100/inf);   % [V/V] Unity Gain "Follower" Buffer Amp. Nominally 1.00
       
       R102 = 510;              % [V/A] ( = R103)
       Gb2  = (1 + R102/inf);   % [V/V] Unity Gain "Follower" Buffer Amp. Nominally 1.00
       
       R3  = 3.3e3;             % [V/A] ( = R20)
       R10 = 2.2e3;             % [V/A] ( = R15)
       Gout = (1 + R3/R10);     % [V/V] Non-inverting Amp. Gain. Nominally 1.33
       
       coilDriver.dcVoltageGain = Gb1 * Gb2 * Gout; % [Vin/Vout] 
       
       % Calculate DC output impedance, Zout
       R104 = 10;
       % R5 = 3.9e3;
       % R5 = 4.3e3; % 10% increase, to better match transconductance data (see LLO aLOG 4495)
       R5 = 2e3; % See T1400223, as per ECR E1400164
       
       
       coilDriver.dcZout = R104 + R5;           % [V/A] Equivalent output resistance of two parallel paths
                                                %       Note that R4 + C12 path is not included, because its 
                                                %       a low pass filter with freq.resp [z:p] = [300:50] 
                                                %       and unity gain.
       
       % Calculate total Current Gain (Ic / Vin), in [A/V]
       coilDriver.dcZc       = Rc.bosem;        % [V/A] At DC, we can ignore 
                                                % the inductance, so Zc = Rc 
                                                
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
        
        
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = 300;
       poles = 85;   % from 50, see T1400223
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 
       zeros = 10.5;
       poles = 1;
       filterModules(2) = fillFilter(zeros,poles,['SimLP1' isoStage]);

       %FM module 7
       zeros = 1;
       poles = 10.5;
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP1' isoStage]);
        
       
       %Fill in FM module 3 
       zeros = 10.5;
       poles = 1;
       filterModules(3) = fillFilter(zeros,poles,['SimLP2' isoStage]);

       %FM module 8
       zeros = 1;
       poles = 10.5;
       filterModules(8)  = fillFilter(zeros,poles,['AntiLP2' isoStage]);
        
       %Fill in FM module 4 
       zeros = 10.5;
       poles = 1;
       filterModules(4) = fillFilter(zeros,poles,['SimLP2' isoStage]);

       %FM module 9
       zeros = 1;
       poles = 10.5;
       filterModules(9)  = fillFilter(zeros,poles,['AntiLP2' isoStage]);
       
        switch coilDriverState
            case 'acq' % DEFAULT Mode (NO relays energized, all B I/O Switched to 0 (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];

            case 'lp1' % LP1 ON (LP1 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
          
               
            case 'lp2' % LP2 ON (LP1 and LP2 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
               
                
            case 'lp3' % LP3 ON (LP1, LP2, and LP3 relay energized, B I/O Switched to 1 (one))
               modulesOn = [6 7 8 9];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
               
            case 'None'
               
        end
        
        
    case {'QUADL2'}
        
       % Circuit Schematic  D070483
       % Design Description T0900277
       % User's Guide       T0900290
       % Board Test Plan    T0900291
       % Chassis Test Plan  T0900292 
       
       coilDriver.type  = 'PUM';
       coilDriver.modes = {'acqOfflpOff','acqOnlpOff','acqOfflpOn','acqOnlpOn'};
       
       
       %                                             {gnd-|Inf|-}+--|Zrc3|--+
       %      gnd--|R1|--+--|Zrc1|--+                            |          |
       %                 |          |                            +-|N       |
       %                 +-|N       |                              |    >---+--- o --|R130|--+---------|R22|---------+-----+  
       %                   |    >---+--|Zrc2|--+-------------------|P            |           |                       |     |
       % o ----------------|P                  |                                 |           +---------|Zdc|---------+     | 
       % |                                     |                                 |           |                       |     |
       % |                                     |                                 |           +---+--|R24|--+         |     |     
       % |                                    _|_                                |              /          |--|Zac|--+    _|_
       % |                                   |   |                               |               +---------+             |   | |
       % Vin                                 |Zsw|                              Vout                                     |Zc | | Ic
       % |                                   |_1_|                               |               +---------+             |___| V
       % |                                     |                                 |              \          |--|Zac|--+     |
       % |                                     |                                 |           +---+--|R24|--+         |     |
       % |                                     |                                 |           |                       |     | 
       % o ----------------|P                  |                                 |           +---------|Zdc|---------+     |
       %                   |    >---+--|Zrc2|--+-------------------|P            |           |                       |     |
       %                 +-|N       |                              |    >---+--- o --|R29 |--+---------|R27|---------+-----+     
       %                 |          |                            +-|N       |
       %      gnd--|R2|--+--|Zrc1|--+                            |          |
       %                                             {gnd-|Inf|-}+--|Zrc3|--+
       %
       %  |------------ Gb1 -----------|-no DC gain-|----------- Gout ---------|   |-------------- Zout -------------|
       %  |------------------------------- Gtot -------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R1   = 10e3;             % [V/A] ( = R1 )
       R12  = 2e3;              % [V/A] ( = R13) At DC, C100 has no effect, so Zrc1 = R12
       Gb1  = (1 + R12/R1);     % [V/V] Non-inverting Amp Gain. Nominally 1.2 
       
       R28  = 100;              % [V/A] ( = R30) At DC, C21 has no effect, so Zrc3 = R28
       Gout = (1 + R28/inf);    % [V/V] Unity Gain "Follower" Buffer Amp. Nominally 1.00
       
       coilDriver.dcVoltageGain = Gb1 * Gout; % [Vin/Vout] 
       
       % Calculate DC output impedance, Zout
       R130 = 15;
       % R22 = 2.2e3;
       R22 = 2.21e3; % increased by 0.4% to better match data, See LLO aLOG 4495.
       
       % [V/A] Equivalent output resistance 
       coilDriver.dcZout = R130 + R22; % R34 (and R24 when switched in) are ignored, matches data
       
       
       % Calculate total Current Gain (Ic / Vin), in [A/V]
       coilDriver.dcZc = Rc.aosem;    % [V/A] At DC, we can ignore 
                                      % the inductance, so Zc = Rc 
                                          
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
       
       
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6 
       zeros = 110;
       poles = 12;
       filterModules(6) = fillFilter(zeros,poles,['AntiAcqOff' isoStage]);
            
       %Fill in FM module 1 
       zeros = 12;
       poles = 110;
       filterModules(1) = fillFilter(zeros,poles,['SimAcqOff' isoStage]);

       %FM module 7
       zeros = 80.5;
       poles = 1.35;
       filterModules(7)  = fillFilter(zeros,poles,['AntiAcqOn' isoStage]);
       
       %FM module 2
       zeros = 1.35;
       poles = 80.5;
       filterModules(2)  = fillFilter(zeros,poles,['SimAcqOn' isoStage]);
       
       %FM module 8
       zeros = [0.5 250];
       poles = [6 20];
       filterModules(8)  = fillFilter(zeros,poles,['AntiLP' isoStage]);
       
       %FM module 3
       zeros = [6 20];
       poles = [0.5 250];
       filterModules(3)  = fillFilter(zeros,poles,['SimLP' isoStage]);
       
       switch coilDriverState
           case 'acqOfflpOff' % DEFAULT Mode (NO relays energized, all B I/O Switched to 0 (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           
           case 'acqOnlpOff' % Acquire Mode (Acquire relays energized, acquire B I/O Switched to 1 (one))
               modulesOn = [7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOfflpOn' %Low pass enabled (LP relay energized, LP B I/O Switched to 1 (one))
               modulesOn = [6 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOnlpOn' % Acquire Mode with Low Pass enabled (Acquire and LP relays energized, Acquire and LP B I/O Switched to 1 (one))
               modulesOn = [7 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];       
       end
       
    case {'QUADL3'}
        % OK OK, it's not coil driver, but there's no sense in making a whole
        % new function for the ESD driver.
        
        coilDriver.type  = 'ESD';
        coilDriver.modes = {'acq20mm','ln20mm','acq5mm','ln5mm'};
        
        switch coilDriverState
            case {'acq20mm','acq5mm'}
                coilDriver.dcVoltageGain = 40;
            case {'ln20mm','ln5mm'}
                coilDriver.dcVoltageGain = 1.1;
        end
        
        
        coilDriver.dcZout = 1;             % Meaningless non-sense just to satisfy for loop
        coilDriver.dcZc   = 1;             % Meaningless non-sense just to satisfy for loop
        coilDriver.dcAdmittance = 1;       % Meaningless non-sense just to satisfy for loop [A_coil / V_out]                                          
        coilDriver.dcTransConductance = 1; % Meaningless non-sense just to satisfy for loop [A_coil / V_in] 
                                          
        
        % Fills in the associated filter modules
        % Create the empty filter module strcutures
        filterModules = repmat(emptyFilter,1,10);
        
        coilDriver.poles_Hz = 2e3;
        coilDriver.zeros_Hz = [];
        coilDriver.zeros_radps = -2*pi*coilDriver.zeros_Hz;
        coilDriver.poles_radps = -2*pi*coilDriver.poles_Hz;

        coilDriver.norm = prod(abs(coilDriver.poles_radps(~(~abs(coilDriver.poles_radps))))) ...
                          ./ prod(abs(coilDriver.zeros_radps(~(~abs(coilDriver.zeros_radps)))));
                      
        coilDriver.c = coilDriver.dcTransConductance*zpk(coilDriver.zeros_radps,... % [[FIX ME]] DC Gain = 1, more meaningless non-sense just to satisfy a for loop.
                                                         coilDriver.poles_radps,...
                                                         coilDriver.norm); 
        coilDriver.d = zpk(c2d(ss(coilDriver.c),...
                         sampleTime,...
                         digitizationType));
    
        coilDriver.f = squeeze(freqresp(coilDriver.d,radFreq));
        
       
    case {'HLTSM2','HLTSM3','HSTSM2','HSTSM3','BSFMM2','BBSSM2','BBSSM3'}
      
       % TRIPLE ACQUISITION DRIVER 
       % Circuit Schematic  D0901047
       % Design Description T1200300
       % User's Guide       T1100192
       % Board Test Plan    T1000133
       % Chassis Test Plan  T??? 
       
       coilDriver.type  = 'Triple Acquisition D0901047-v4';
       coilDriver.modes = {'acqOfflpOff','acqOnlpOff','acqOfflpOn','acqOnLpOn'};
       
       %                                              gnd-|R5|--+--|Zrc2|--+
       %     {gnd-|Inf|-}+--|R22|--+                            |          |
       %                 |         |                            +-|N       |
       %                 +-|N      |                              |    >---+--- o --|R8 |--+----------|R10|---------+-----+
       %                   |    >--+--|Zrc1|--+-------------------|P            |          |                        |     |
       % o ----------------|P                 |                                 |          +----------|Zdc|---------+     |
       % |                                    |                                 |          |                        |     |
       % |                                    |                                 |          +---+--|R14|--+          |     |
       % |                                   _|_                                |             /          |--|Cnet|--+    _|_
       % |                                  |   |                               |              +---------+             |   | |
       % Vin                                |Zsw|                              Vout                                    |Zc | | Ic
       % |                                  |___|                               |              +---------+             |___| V
       % |                                    |                                 |             \          |--|Cnet|--+     |
       % |                                    |                                 |          +---+--|R15|--+          |     |
       % |                                    |                                 |          |                        |     |
       % o ----------------|P                 |                                 |          +----------|Zdc|---------+     |
       %                   |    >--+--|Zrc1|--+-------------------|P            |          |                        |     |
       %                 +-|N      |                              |    >---+--- o --|R9 |--+----------|R11|---------+-----+
       %                 |         |                            +-|N       |
       %     {gnd-|Inf|-}+--|R24|--+                            |          |
       %                                             gnd--|R6|--+--|Zrc2|--+
       %
       %  |------------ Gb1 ---------|-no DC gain-|----------- Gout -----------|   |-------------- Zout -------------|
       %  |------------------------------- Gtot -------------------------------|
       
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R22  = 499;              % [V/A] ( = R24 )
       Gb1  = (1 + R22 / inf);  % [V/V] Non-inverting Amp Gain. Nominally 1.2 
       
       % High current unity gain follower IC3's input is reduced by half because
       % of OPA544F's input range, but then immediately doubled to recover unity
       % gain. Just in case any of the resistances change, we actually compute the
       % gain (HAM-A driver's the only other one like this, the rest [TOP, UIM,
       % PUM] all just have a normal unity gain follower). This detail is not
       % shown in the toy schematic above merely for consistency.
       % Gic3in
       R29 = 1e3;    % [V/A] (= R32)
       R30 = 1e3;    % [V/A] (= R33)
       Gic3in = R30./(R29 + R30); % [V/V] voltage divider. nominally 0.5 at all frequencies.
       
       % Gic3out
       R12 = 1e3;    % [V/A] (= R20)
       R28 = 1e3;    % [V/A] (= R31)
       Gic3out = (1 + R12/R28);   % [V/V] non-inverting op amp. nominally 2.0 at all frequencies.
       
       % Gic3tot
       Gic3tot = Gic3in * Gic3out; % [V/V] total gain, nominally 1.0 at all frequencies.
       
       % Gout
       R4  = 1.21e3; % [V/A] (= R7)
       R5 = 1e3;     % [V/A] (= R6)
       Gout = Gic3tot*(1 + R4./R5); % [V/V] non-inverting op amp. Nominally 2.21 at DC.
       
       coilDriver.dcVoltageGain = Gb1 * Gout; % [Vin/Vout] 
       
       R8   = 56;    % [V/A] (= R9)
       R21  = 332;   % [V/A] (= R27)
       
       switch susType % Stupid stupid stupid stupid
           case {'HSTS','HLTS','HSTSMOD'}
               R10  = 3.32e3;% [V/A] (= R11)
               coilDriver.dcZc = Rc.aosem;    % [V/A] At DC, we can ignore 
                                              % the inductance, so Zc = Rc 
           case {'BSFM','BBSS'}
               R10  = 3.3e3;% [V/A] (= R11)
               coilDriver.dcZc = Rc.bosem;    % [V/A] At DC, we can ignore 
                                              % the inductance, so Zc = Rc 
       end
       
       coilDriver.dcZout = R8 + R10; % [V/A] (All other paths involve capacitors which are open at DC)
       
       % Calculate total Current Gain (Ic / Vin), in [A/V]
                                          
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
       
       
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6 
       zeros = 82;
       poles = 9;
       filterModules(6) = fillFilter(zeros,poles,['AntiAcqOff' isoStage]);
            
       %Fill in FM module 1 
       zeros = 9;
       poles = 82;
       filterModules(1) = fillFilter(zeros,poles,['SimAcqOff' isoStage]);

       %FM module 7
       zeros = 46;
       poles = 1.05;
       filterModules(7)  = fillFilter(zeros,poles,['AntiAcqOn' isoStage]);
       
       %FM module 2
       zeros = 1.05;
       poles = 46;
       filterModules(2)  = fillFilter(zeros,poles,['SimAcqOn' isoStage]);
       
       %FM module 8
       zeros = [1,210];
       poles = [11,21];
       filterModules(8)  = fillFilter(zeros,poles,['AntiLP' isoStage]);
       
       %FM module 3
       zeros = [11,21];
       poles = [1,210];
       filterModules(3)  = fillFilter(zeros,poles,['SimLP' isoStage]);
       
       switch coilDriverState
           case 'acqOfflpOff' % DEFAULT Mode (NO relays energized, all B I/O Switched to 0 (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           
           case 'acqOnlpOff' % Acquire Mode (Acquire relays energized, acquire B I/O Switched to 1 (one))
               modulesOn = [7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOfflpOn' %Low pass enabled (LP relay energized, LP B I/O Switched to 1 (one))
               modulesOn = [6 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOnlpOn' % Acquire Mode with Low Pass enabled (Acquire and LP relays energized, Acquire and LP B I/O Switched to 1 (one))
               modulesOn = [7 8];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];       
       end
       
    case {'HSTSMODM2','HSTSMODM3','HLTSMODM2','BSFMMODM2'}
               
       % TRIPLE ACQUISITION DRIVER (Modified for 10x gain)
       % Circuit Schematic  L1200226
       % Design Description T1200300
       % User's Guide       ????
       % Board Test Plan    T1000133
       % Chassis Test Plan  T??? 
       
       coilDriver.type  = 'Triple Acquisition L1200226 Mods';
       coilDriver.modes = {'acqOfflpOff','acqOnlpOff','acqOfflpOn','acqOnLpOn'};
       
              %                                       gnd-|R5|--+--|Zrc2|--+                      L1200226 Short Here
       %     {gnd-|Inf|-}+--|R22|--+                            |          |                               |
       %                 |         |                            +-|N       |                               V
       %                 +-|N      |                              |    >---+--- o --|R8 |--+------|R10|----+--------+-----+
       %                   |    >--+--|Zrc1|--+-------------------|P            |          |               |        |     |
       % o ----------------|P                 |                                 |          +------|R21|----+-|C28|--+     |
       % |                                    |                                 |          |                        |     |
       % |                                    |                                 |          +---+--|R14|--+          |     |
       % |                                   _|_                                |             /          |--|Cnet|--+    _|_
       % |                                  |   |                               |              +---------+             |   | |
       % Vin                                |Zsw|                              Vout                                    |Zc | | Ic
       % |                                  |___|                               |              +---------+             |___| V
       % |                                    |                                 |             \          |--|Cnet|--+     |
       % |                                    |                                 |          +---+--|R15|--+          |     |
       % |                                    |                                 |          |                        |     |
       % o ----------------|P                 |                                 |          +------|Zdc|----+-|C37|--+     |
       %                   |    >--+--|Zrc1|--+-------------------|P            |          |               |        |     |
       %                 +-|N      |                              |    >---+--- o --|R9 |--+------|R11|----+--------+-----+
       %                 |         |                            +-|N       |
       %     {gnd-|Inf|-}+--|R24|--+                            |          |
       %                                             gnd--|R6|--+--|Zrc2|--+
       %
       %  |------------ Gb1 ---------|-no DC gain-|----------- Gout -----------|   |-------------- Zout -------------|
       %  |------------------------------- Gtot -------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R22  = 499;              % [V/A] ( = R24 )
       Gb1  = (1 + R22 / inf);  % [V/V] Non-inverting Amp Gain. Nominally 1.0 
       
       % High current unity gain follower IC3's input is reduced by half because
       % of OPA544F's input range, but then immediately doubled to recover unity
       % gain. Just in case any of the resistances change, we actually compute the
       % gain (HAM-A driver's the only other one like this, the rest [TOP, UIM,
       % PUM] all just have a normal unity gain follower). This detail is not
       % shown in the toy schematic above merely for consistency.
       % Gic3in
       R29 = 1e3;    % [V/A] (= R32)
       R30 = 1e3;    % [V/A] (= R33)
       Gic3in = R30./(R29 + R30); % [V/V] voltage divider. nominally 0.5 at all frequencies.
       
       % Gic3out
       R12 = 1e3;    % [V/A] (= R20)
       R28 = 1e3;    % [V/A] (= R31)
       Gic3out = (1 + R12/R28);   % [V/V] non-inverting op amp. nominally 2.0 at all frequencies.
       
       % Gic3tot
       Gic3tot = Gic3in * Gic3out; % [V/V] total gain, nominally 1.0 at all frequencies.
       
       % Gout
       R4  = 1.21e3; % [V/A] (= R7)
       R5 = 1e3;     % [V/A] (= R6)
       Gout = Gic3tot*(1 + R4./R5); % [V/V] non-inverting op amp. Nominally 2.21 at DC.
       
       coilDriver.dcVoltageGain = Gb1 * Gout; % [Vin/Vout] 
       
       R8   = 56;    % [V/A] (= R9)
       % R21  = 332;   % [V/A] (= R27)
       R21  = 360;   % [V/A] (= R27) % increased by 10% to match data, see LLO aLOG 4495
       
       switch susType % Stupid stupid stupid stupid
           case {'HSTS','HLTS','HSTSMOD'}
               R10  = 3.32e3;% [V/A] (= R11)
           case 'BSFM'
               R10  = 3.3e3;% [V/A] (= R11)
       end
       
       coilDriver.dcZout = R8 + (1./R10 + 1/R21)^(-1); % [V/A] (All other paths involve capacitors which are open at DC)
       
       % Calculate total Current Gain (Ic / Vin), in [A/V]
       coilDriver.dcZc = Rc.aosem;    % [V/A] At DC, we can ignore 
                                      % the inductance, so Zc = Rc 
                                          
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
       
       
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6 
       zeros = 65;
       poles = 13;
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 1 
       zeros = 13;
       poles = 65;
       filterModules(1) = fillFilter(zeros,poles,['SimAcq' isoStage]);

       %FM module 7
       zeros = [1,210];
       poles = [11,21];
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);
       
       %FM module 2
       zeros = [11,21];
       poles = [1,210];
       filterModules(2)  = fillFilter(zeros,poles,['SimLP' isoStage]);
       
       switch coilDriverState
           case 'acqOfflpOff' % DEFAULT Mode (NO relays energized, all B I/O Switched to 0 (zero))
               modulesOn = [];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           
           case 'acqOnlpOff' % Acquire Mode (Acquire relays energized, acquire B I/O Switched to 1 (one))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOfflpOn' %Low pass enabled (LP relay energized, LP B I/O Switched to 1 (one))
               modulesOn = [7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];
           case 'acqOnlpOn' % Acquire Mode with Low Pass enabled (Acquire and LP relays energized, Acquire and LP B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriverState];       
       end 
       
    case {'HRTSM1','HATSM1','BHSSM1'}
       
       % HAM-A DRIVER (100 Ohm output impedance)
       % Assembly Drawing     D1100687-v1
       % Circuit Schematic    D1100117-v2
       % Design Description   T1200264
       % User's Guide         ????????
       % Board Test Plan      ????????
       % Chassis Test Plan    ???????? 
       
       coilDriver.type       = 'HAM-A D1100117-v2';
       coilDriver.modes      = {'acq';'lp'};
       
       %                                               gnd-|R7 |--+--|R5 |--+
       %                                                          |         |
       %                                                          +-|N      |
       % o ---------------------|P                                  |    >--+--- o ----|R11|-----+
       % |                      |    >--+--|Zrc1|--+----------------|P           |               |
       % |                    +-|N      |          |                             |               |
       % |                    |         |         _|_                            |              _|_
       % |           +-|inf|-}+--|R9 |--+        |   |                           |             |   | |
       % Vin   {gnd--|                           |Zrc|                          Vout           |Zc | | Ic
       % |           +-|inf|-}+--|R21|--+        |_2_|                           |             |___| V
       % |                    |         |          |                             |               |
       % |                    +-|P      |          |                             |               |
       % |                      |    >--+--|Zrc1|--+----------------|P           |               |
       % o ---------------------|N                                  |    >--+--- o ----|R33|-----+
       %                                                          +-|N      |
       %                                                          |         |
       %                                               gnd-|R24|--+--|R23|--+
       %
       %  |-------------- Gb1 -----------|-no DC gain-|---------- Gout ---------| |--- Zout ----|
       %  |--------------------------------- Gtot ------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R9  = 499 ;                          % [V/A] ( = R21)
       Gb1 = (1 + R9/inf);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       R5  = 1.21e3;                        % [V/A] ( = R23)
       R7  = 1.00e3;                        % [V/A] ( = R24)
       Gout = (1 + R5/R7);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       coilDriver.dcVoltageGain = Gb1 * Gout;        % [Vin/Vout] 
       
       R11 = 100;       % [V/A] (= R33) 

       coilDriver.dcZout = R11; % [V/A]
       
       % Calculate Total DC Current Gain (Ic / Vin), in [A/V] 
       switch susType
           case {'HTTS','HSDS','HPDS','HTDS','HDDS','HXDS','HRTS'}
               coilDriver.dcZc            = Rc.bosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
           case {'HAUX','OFIS','OPOS','BHSS'}
               coilDriver.dcZc            = Rc.aosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
           case {'HATS'}
               switch isoStage
                   case {'M1'}
                       coilDriver.dcZc            = Rc.bosem;      % [V/A] At DC, we can ignore
                                                                   % the inductance, so Zc = Rc
                   case {'M3'}
                       coilDriver.dcZc            = Rc.aosem;      % [V/A] At DC, we can ignore
                                                                   % the inductance, so Zc = Rc
               end
       end
                                                   
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
    
                                     
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = [0.031 20.65]; % [Hz]
       poles = [0.032 20.00]; % [Hz]
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 (plus one because of matlab indexing)
       zeros = [11 21];   % [Hz]
       poles = [1 215]; % [Hz]
       filterModules(2) = fillFilter(zeros,poles,['SimLP' isoStage]);

       %FM module 7
       zeros = [1 215];
       poles = [11 21];
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);

       % Sets total coil
       switch coilDriverState
           case coilDriver.modes(1)
               % ACQUIRE Mode (Relays NOT Energized, B I/O Switched to 0
               % (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(1)];
           case coilDriver.modes(2)
               % RUN Mode (Relays Energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);              
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(2)];
           case 'None'
               
           otherwise
               error(['Unrecognized coilDriverState "' coilDriverState '" for osemType "' osemType '"'])
       end                  
       

    case {'HTTSM1','HSDSM1','HPDSM1','HTDSM1','HDDSM1','HXDSM1','HCDSM1','HAUXM1','OFISM1','OPOSM1','HATSM3'}
       
       % HAM-A DRIVER (1.2k Ohm output impedance)
       % Assembly Drawing     D1100687-v2
       % Circuit Schematic    D1100117-v3 or -v4
       % Design Description   T1200264
       % User's Guide         ????????
       % Board Test Plan      ????????
       % Chassis Test Plan    ???????? 
       
       coilDriver.type       = 'HAM-A D1100117-v3 or -v4';
       coilDriver.modes      = {'acq';'lp'};
       
       %                                               gnd-|R7 |--+--|R5 |--+
       %                                                          |         |
       %                                                          +-|N      |
       % o ---------------------|P                                  |    >--+--- o ----|R11|-----+
       % |                      |    >--+--|Zrc1|--+----------------|P           |               |
       % |                    +-|N      |          |                             |               |
       % |                    |         |         _|_                            |              _|_
       % |           +-|inf|-}+--|R9 |--+        |   |                           |             |   | |
       % Vin   {gnd--|                           |Zrc|                          Vout           |Zc | | Ic
       % |           +-|inf|-}+--|R21|--+        |_2_|                           |             |___| V
       % |                    |         |          |                             |               |
       % |                    +-|P      |          |                             |               |
       % |                      |    >--+--|Zrc1|--+----------------|P           |               |
       % o ---------------------|N                                  |    >--+--- o ----|R33|-----+
       %                                                          +-|N      |
       %                                                          |         |
       %                                               gnd-|R24|--+--|R23|--+
       %
       %  |-------------- Gb1 -----------|-no DC gain-|---------- Gout ---------| |--- Zout ----|
       %  |--------------------------------- Gtot ------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R9  = 499 ;                          % [V/A] ( = R21)
       Gb1 = (1 + R9/inf);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       R5  = 1.21e3;                        % [V/A] ( = R23)
       R7  = 1.00e3;                        % [V/A] ( = R24)
       Gout = (1 + R5/R7);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       coilDriver.dcVoltageGain = Gb1 * Gout;        % [Vin/Vout] 
       
       R11 = 1.2e3;       % [V/A] (= R33) 

       coilDriver.dcZout = R11; % [V/A]
       
       % Calculate Total DC Current Gain (Ic / Vin), in [A/V] 
       switch osemType
           case {'HTTSM1','HSDSM1','HPDSM1','HTDSM1','HDDSM1','HXDSM1','HCDSM1','HATSM1'}
               coilDriver.dcZc            = Rc.bosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
           case {'HAUXM1','OFISM1','OPOSM1','HATSM3'}
               coilDriver.dcZc            = Rc.aosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
       end
                                                   
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
    
                                     
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = [0.031 20.65]; % [Hz]
       poles = [0.032 20.00]; % [Hz]
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 (plus one because of matlab indexing)
       zeros = [11 21];   % [Hz]
       poles = [1 215]; % [Hz]
       filterModules(2) = fillFilter(zeros,poles,['SimLP' isoStage]);

       %FM module 7
       zeros = [1 215];
       poles = [11 21];
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);

       % Sets total coil
       switch coilDriverState
           case coilDriver.modes(1)
               % ACQUIRE Mode (Relays NOT Energized, B I/O Switched to 0
               % (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(1)];
           case coilDriver.modes(2)
               % RUN Mode (Relays Energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);              
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(2)];
           case 'None'
               
           otherwise
               error(['Unrecognized coilDriverState "' coilDriverState '" for osemType "' osemType '"'])
       end           
      
              
    case {'HTTSMODM1','HSDSMODM1','HPDSMODM1','HTDSMODM1','HDDSMODM1','HXDSMODM1','HAUXMODM1','OPOSMODM1','OFISMODM1'}
       
       % HAM-A DRIVER (modified low-pass)
       % Circuit Schematic    ????????
       % Design Description   T2100410-v3 or higher
       % User's Guide         ????????
       % Board Test Plan      ????????
       % Chassis Test Plan    ???????? 
       
       coilDriver.type       = 'HAM-A MOD T2100410-v3';
       coilDriver.modes      = {'acq';'lp'};
       
       %                                               gnd-|R7 |--+--|R5 |--+
       %                                                          |         |
       %                                                          +-|N      |
       % o ---------------------|P                                  |    >--+--- o ----|R11|---+
       % |                      |    >--+--|Zrc1|--+----------------|P           |             |
       % |                    +-|N      |          |                             |             |
       % |                    |         |         _|_                            |            _|_
       % |           +-|inf|-}+--|R9 |--+        |   |                           |           |   | |
       % Vin   {gnd--|                           |Zrc|                          Vout         |Zc | | Ic
       % |           +-|inf|-}+--|R21|--+        |_2_|                           |           |___| V
       % |                    |         |          |                             |             |
       % |                    +-|P      |          |                             |             |
       % |                      |    >--+--|Zrc1|--+----------------|P           |             |
       % o ---------------------|N                                  |    >--+--- o ----|R33|---+
       %                                                          +-|N      |
       %                                                          |         |
       %                                               gnd-|R24|--+--|R23|--+
       %
       %  |-------------- Gb1 -----------|-no DC gain-|---------- Gout ---------| |--- Zout ---|
       %  |--------------------------------- Gtot ------------------------------|
       
       % Calculate total DC Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       R9  = 499 ;                          % [V/A] ( = R21)
       Gb1 = (1 + R9/inf);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       R5  = 1.21e3;                        % [V/A] ( = R23)
       R7  = 1.00e3;                        % [V/A] ( = R24)
       Gout = (1 + R5/R7);                  % [V/V] Non-inverting Amp. Gain. Nominally 1.00
       
       coilDriver.dcVoltageGain = Gb1 * Gout;        % [Vin/Vout] 
       
       R11 = 1.2e3;       % [V/A] (= R33) 

       coilDriver.dcZout = R11; % [V/A]
       
       % Calculate Total DC Current Gain (Ic / Vin), in [A/V] 
       switch susType
           case {'HTTS','HSDS','HPDS','HTDS','HDDS','HXDS'}
               coilDriver.dcZc            = Rc.bosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
           case {'HAUX','OFIS','OPOS'}
               coilDriver.dcZc            = Rc.aosem;      % [V/A] At DC, we can ignore
                                                           % the inductance, so Zc = Rc
       end
                                                   
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]                                          
       coilDriver.dcTransConductance   = coilDriver.dcVoltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.dcAdmittance; 
    
                                     
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = [0.031   50.77]; % [Hz]
       poles = [0.03305 54.10]; % [Hz]
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 (plus one because of matlab indexing)
       zeros = [50.77 52.32];   % [Hz]
       poles = [0.5 3174.41]; % [Hz]
       filterModules(2) = fillFilter(zeros,poles,['SimLP' isoStage]);

       %FM module 7
       zeros = [0.5 3174.41];
       poles = [50.77 52.32];
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);

       % Sets total coil
       switch coilDriverState
           case coilDriver.modes(1)
               % ACQUIRE Mode (Relays NOT Energized, B I/O Switched to 0
               % (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(1)];
           case coilDriver.modes(2)
               % RUN Mode (Relays Energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);              
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(2)];
           case 'None'
               
           otherwise
               error(['Unrecognized coilDriverState "' coilDriverState '" for osemType "' osemType '"'])
       end                              
                                     
    case {'HDDSM2','HCDSM2'}
       
       % HAM-A DITHER DRIVER
       % Assembly Drawing     D1100687-v3
       % Circuit Schematic    D2100489
       % Design Description   T1200264
       % User's Guide         ????????
       % Board Test Plan      ????????
       % Chassis Test Plan    ???????? 
       
       coilDriver.type       = 'HAM-A_Dither';
       coilDriver.modes      = {'acq'};
       
       %                                               gnd-|R7 |--+--|R5 |--+
       %                                                          |         |
       %                                                          +-|N      |
       % o ---------------------|P                                  |    >--+--- o ----|R11|--|C26|---+
       % |                      |    >--+--|Zrc1|--+----------------|P           |                    |
       % |                    +-|N      |          |                             |                    |
       % |                    |         |         _|_                            |                   _|_
       % |           +-|inf|-}+--|R9 |--+        |   |                           |                  |   | |
       % Vin   {gnd--|                           |Zrc|                          Vout                |Zc | | Ic
       % |           +-|inf|-}+--|R21|--+        |_2_|                           |                  |___| V
       % |                    |         |          |                             |                    |
       % |                    +-|P      |          |                             |                    |
       % |                      |    >--+--|Zrc1|--+----------------|P           |                    |
       % o ---------------------|N                                  |    >--+--- o ----|R33|--|C27|---+
       %                                                          +-|N      |
       %                                                          |         |
       %                                               gnd-|R24|--+--|R23|--+
       %
       %  |-------------- Gb1 -----------|--- GSwitch ---|---------- Gout ---------| |------ Zout -------|
       %  |--------------------------------- Gtot ------------------------------|
       
       % Calculate total Voltage Gain, Gtot in [Vin/Vout]
       %     Analyzed by following positive leg, with anything crossing the
       %     legs (measured voltages and impedances) halved and connected to
       %     gnd, to be consistent. 
       %     Works because [1/2 Vout / 1/2 Vin] = [Vin/Vout].
       
       w = 2*pi*freq;
       
       % Gswitch
       C11 = 10.0e-9; % Farad
       R18 = 10.0e3; % Ohm
       
       % Gout
       R5 = 1.21e3;
       R7 = 1.0e3;
       
       % Former switchable low pass is now a high pass filter,
       % using a simple voltage divider configuration
       Gswitch = (R18) ./ ( (R18) + (1./(1i*w*C11)) );
       
       % The usual overly complicated output opamp
       % system can just be treated as a single
       % non-inverting opamp
       Gout = (1.0 + R5 / R7);
       
       coilDriver.voltageGain = Gswitch * Gout;        % [Vin/Vout] 
       
       % HAM-A Output Impedance
       Rout = 100; % Ohm
       Cout = 2.2e-6; % Farad
       
       % Former RC para output impedance is now an
       % RC series impedance, and coupled with the
       % coil series LC is an one zero, two pole RLC
       % circuit
       coilDriver.Zout = Rout + 1./(1i*w*Cout);
       
       coilDriver.dcZout = Rout; % [V/A]
       
       coilDriver.dcZc = Rc.cosem;
       
       coilDriver.Zc = Rc.cosem + 1i*w*Lc.cosem;
                                                   
       coilDriver.dcAdmittance = 1./(2*coilDriver.dcZout + coilDriver.dcZc); % [A_coil / V_out]    
       coilDriver.admittance = 1./(2*coilDriver.Zout + coilDriver.Zc);
       coilDriver.transConductance   = coilDriver.voltageGain .* ...     % [A_coil / V_in] 
                                         coilDriver.admittance; 
       
       %%%%%%%
       % TBD %
       %%%%%%%
       % Just putting in junk for now to make the code work:
       
       % Fills in the associated filter modules
       %Create the empty filter module strcutures
       filterModules = repmat(emptyFilter,1,10);
       
       %Fill in FM module 6
       zeros = [1000 1000]; % [Hz]
       poles = [0 0]; % [Hz]
       filterModules(6) = fillFilter(zeros,poles,['AntiAcq' isoStage]);
            
       %Fill in FM module 2 (plus one because of matlab indexing)
       zeros = [0 0];   % [Hz]
       poles = [1000 1000]; % [Hz]
       filterModules(2) = fillFilter(zeros,poles,['SimLP' isoStage]);

       %FM module 7
       zeros = [1000 1000];
       poles = [0 0];
       filterModules(7)  = fillFilter(zeros,poles,['AntiLP' isoStage]);

       % Sets total coil
       switch coilDriverState
           case coilDriver.modes(1)
               % ACQUIRE Mode (Relays NOT Energized, B I/O Switched to 0
               % (zero))
               modulesOn = [6];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(1)];
           case coilDriver.modes(2)
               % RUN Mode (Relays Energized, B I/O Switched to 1 (one))
               modulesOn = [6 7];
               coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn);              
               coilDriver.fmName = [isoStage ' ' coilDriver.modes(2)];
           case 'None'
               
           otherwise
               error(['Unrecognized coilDriverState "' coilDriverState '" for osemType "' osemType '"'])
       end
   coilDriver.Rc = Rc;
   coilDriver.Lc = Lc;
end

function satAmp = getsatamp(susType,isoStage)
global radFreq

    osemType = [susType isoStage];

    switch osemType
        % Satellite Amplifier, D0901284 (same for all suspensions, all types, A
        % or B OSEMs; thank god.)
        case {'QUADM0','QUADR0','QUADL1','QUADL2',...
              'BSFMM1','BSFMM2','HLTSM1','HLTSM2',...
              'HLTSM3','HSTSM1','HSTSM2','HSTSM3',...
              'TMTSM1','OMCSM1','HTTSM1','HAUXM1',...
              'HSDSM1','HPDSM1','HTDSM1','HDDSM1','HCDSM1',...
              'HXDSM1','OFISM1','OPOSM1',...
              'HSTSMODM1','HSTSMODM2','HSTSMODM3',...
              'HLTSMODM1','HLTSMODM2','HLTSMODM3',...
              'BSFMMODM1','BSFMMODM2','QUADMODL1',...
              'HTTSMODM1','HAUXMODM1',...
              'HSDSMODM1','HPDSMODM1','HTDSMODM1','HDDSMODM1',...
              'HXDSMODM1','OFISMODM1','OPOSMODM1',...
              'BBSSM1','BBSSM2','BBSSM3',...
              'HRTSM1',...
              'HATSM1','HATSM3',...
              'BHSSM1'}
            
            satAmp.dcGain = 240e3;             % [V/A] (Rx02 gain resistor * 2 for single-end to differential)
            
            satAmp.poles_Hz = 10;
            satAmp.zeros_Hz = 0.4;
            satAmp.poles_radps = -2*pi*satAmp.poles_Hz; % [Radians]
            satAmp.zeros_radps = -2*pi*satAmp.zeros_Hz;      % [Radians]
            satAmp.norm = prod(satAmp.poles_radps) / prod(satAmp.zeros_radps);
            
            
            satAmp.c = zpk(satAmp.zeros_radps,satAmp.poles_radps,satAmp.norm);
            satAmp.f = squeeze(freqresp(satAmp.c,radFreq));
            
        case {'QUADL3','BSFMM3','OMCSM2','TMTSM2','HDDSM2','HRTSM2','HATSM2','HCDSM2'}
            fprintf(['WARNING: susType ''%s'' and isoStage ''%s'' does not have OSEMs. Setting satAmp gain to 1.\n\n'],susType,isoStage)
            
            satAmp.dcGain = 1;
            
            satAmp.poles_Hz = [];
            satAmp.zeros_Hz = [];
            satAmp.poles_radps = -2*pi*satAmp.poles_Hz; % [Radians]
            satAmp.zeros_radps = -2*pi*satAmp.zeros_Hz;      % [Radians]
            satAmp.norm = prod(satAmp.poles_radps) / prod(satAmp.zeros_radps);
            
            
            satAmp.c = zpk(satAmp.zeros_radps,satAmp.poles_radps,satAmp.norm);
            satAmp.f = squeeze(freqresp(satAmp.c,radFreq));
        otherwise
            error(['susType ''%s'' and isoStage ''%s'' are not recognized to have the standard satamp.\n'],susType,isoStage)
            
    end

function out = getdacbitcount(susType,isoStage,in)
    osemType = [susType isoStage];
    switch osemType
        case {'QUADL3'}
            in.dac.nBits = 20;
        case {'HSDSM1','HPDSM1','HDDSM1','HTDSM1','HXDSM1','HTTSM1','OPOSM1','OFISM1',...
              'HSDSMODM1','HPDSMODM1','HDDSMODM1','HTDSMODM1','HXDSMODM1','HTTSMODM1','OPOSMODM1','OFISMODM1'}
            in.dac.nBits = 16;
        otherwise
            in.dac.nBits = 18;
    end
    
    in.dac.V_per_drivect = 20/2^(in.dac.nBits);
    
    out = in;

function filterDef = fillFilter(z,p,name)
    global emptyFilter;
    global digitizationType;
    global sampleTime;
    global radFreq;
    global plotFigs
    
    filterDef = emptyFilter;
    
    filterDef.name = name;
    filterDef.zeros_Hz = z; 
    filterDef.poles_Hz = p; 
    filterDef.zeros_radps = -2*pi*z; 
    filterDef.poles_radps = -2*pi*p; 
    filterDef.norm = prod(filterDef.poles_radps) / prod(filterDef.zeros_radps);
    filterDef.used = true;
    
    filterDef.c = zpk(filterDef.zeros_radps,...    %Convert to Radians from Hz
                        filterDef.poles_radps,...  %Convert to Radians from Hz
                            filterDef.norm);
                        
    filterDef.d = zpk(c2d(ss(filterDef.c),...
                        sampleTime,...
                        digitizationType));
    
    filterDef.f = squeeze(freqresp(filterDef.d,radFreq));
    
    if plotFigs
        figure
        bode(filterDef.c,filterDef.d)
    end
    
function coilDriver = calculateWithInverseModulesOn(coilDriver,filterModules,modulesOn)
    global radFreq
    global digitizationType
    global sampleTime
    
    coilDriver.zeros_Hz = [];
    coilDriver.poles_Hz = [];
    
    for ii = 1:10
        if (filterModules(ii).used == true) && ~isempty(intersect(ii,modulesOn))
            for jj = 1:length(filterModules(ii).zeros_Hz)
                coilDriver.zeros_Hz(length(coilDriver.zeros_Hz)+1) = filterModules(ii).poles_Hz(jj);
            end
            for jj = 1:length(filterModules(ii).poles_Hz)
                coilDriver.poles_Hz(length(coilDriver.poles_Hz)+1) = filterModules(ii).zeros_Hz(jj);
            end
        elseif isempty(modulesOn)
            coilDriver.zeros_Hz = [];
            coilDriver.poles_Hz = [];
        end
    end
    coilDriver.zeros_radps = -2*pi*coilDriver.zeros_Hz;
    coilDriver.poles_radps = -2*pi*coilDriver.poles_Hz;
    coilDriver.norm = (prod(coilDriver.poles_radps) / prod(coilDriver.zeros_radps));
    
    coilDriver.c = zpk(coilDriver.zeros_radps,coilDriver.poles_radps,coilDriver.norm);
    
    coilDriver.d = zpk(c2d(ss(coilDriver.c),...
                         sampleTime,...
                         digitizationType));
    
    coilDriver.f = squeeze(freqresp(coilDriver.d,radFreq));

    
    
    
