% LTPDA_MDC1_C returns a frequency-domain model of the controllers for MDC1.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_MDC1_C returns a frequency-domain model of the controllers 
%              for MDC1.
%
% CALL:        b = ltpda_mdc1_C(pl)
%
% PARAMETERS:
% 
%           'Controller' - Choose controller [default: 'DF']
%                        'DF'   - Drag-free controller (including delay)
%                        'SUS'  - Suspension controller (including delay)
%           'f'        - a vector of frequencies [default: 1]
%    or
%           'f1'       - start frequency [default: 1e-6]
%           'f2'       - stop frequency  [default: 5]
%           'nf'       - number of frequency points [default: 1000]
%           'scale'    - frequency spacing, 'lin' or 'log' [default: 'log']
% 
% VERSION:     $Id: ltpda_mdc1_C.m,v 1.3 2008/08/08 13:35:23 anneke Exp $
%
% The following call returns a parameter list object that contains the
% default parameter values:
%
% >> pl = ltpda_mdc1_C(ao, 'Params')
%
% The following call returns a string that contains the routine CVS version:
%
% >> version = ltpda_mdc1_C(ao,'Version')
%
% The following call returns a string that contains the routine category:
%
% >> category = ltpda_mdc1_C(ao,'Category')
%
% HISTORY: 11-04-08 M Hewitson
%             Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = ltpda_mdc1_C(varargin)
    %%% Check if this is a call for parameters
        if utils.helper.isinfocall(varargin{:})
          varargout{1} = getInfo(varargin{3});
          return
        end

        %%% Collect input variable names
        in_names = cell(size(varargin));
        for ii = 1:nargin,in_names{ii} = inputname(ii);end

        pli             = utils.helper.collect_objects(varargin(:), 'plist', in_names);

        %%% Decide on a deep copy or a modify
        %%% REMARK: If you create a new AO (call the constructor) then
        %%%         it is not necessay to copy the input-AOs !!!!!!!!!!!!!!!!!!!!!!!!!


        %%% Combine plists
        pl = combine(pli, getDefaultPlist('Range'));


    %% Extract parameters from plist

    f = find(pl, 'f');
    if isempty(f)
    f1    = find(pl, 'f1');
    f2    = find(pl, 'f2');
    nf    = find(pl, 'nf');
    scale = find(pl, 'scale');

    switch scale
      case 'lin'
        f   = linspace(f1, f2, nf);
      case 'log'
        f = logspace(log10(f1), log10(f2), nf);
    end
    end

    cont = find(pl, 'Controller');

    %% Compute response for frequencies f

    switch upper(cont)
    case 'DF'
      a = getDFcontroller(f);
    case 'SUS'
      a = getSUScontroller(f);
    otherwise
      error('### Unknown controller requested.');
    end

    varargout{1} = a;
end
%--------------------------------------------------------------------------
% Get DF controller for each frequency
function o = getSUScontroller(f)

  % Compute desired response
  T = 0.1;
  z = exp(1i*2*pi*f*T);
  zi = 1./z;

  f1 = -2.4365e-4./(1-exp(-1.0772e-2).*zi);
  f2 = complex(1.2115e-4, 3.0511e-5)./(1-exp(complex(-5.5589e-3,-6.4448e-3)).*zi);
  f3 = complex(1.2115e-4, -3.0511e-5)./(1-exp(complex(-5.5589e-3, 6.4448e-3)).*zi);
  f4 = 1.0741e-6./(1 - 1.*zi);

  % Controller delay
  s = 1i*2*pi*f;
  dt = 305/1000;
  D = (exp(-s.*(T+dt)) .* (-1+exp(s.*T))  )./(s.*T);

  H = 0.1.*( f1 + f2 + f3 + f4).*D;

  % Create fsdata object
  fsd = fsdata(f, H);

  % Create an AO
  o = ao(fsd);
  o.setXunits('Hz');
  o.setYunits('N/V');
  o.setName('C_{sus}');
end
%--------------------------------------------------------------------------
% Get DF controller for each frequency
function o = getDFcontroller(f)

  % Transform from f to z
  T = 0.1;
  z = exp(1i*2*pi*f*T);
  zi = 1./z;

  % Compute individual terms in the frequency series
  f1 = -1.61./(1-exp(-0.30943).*zi);
  f2 = complex(0.86807, 0.12606)./(1-exp(complex(-9.753e-2,-1.6175e-1)).*zi);
  f3 = complex(0.86807, -0.12606)./(1-exp(complex(-9.753e-2,1.6175e-1)).*zi);
  f4 = -complex(6.4699,70.435)./(1-exp(complex(-5.5289e-5,-4.8946e-6)).*zi);
  f5 = -complex(6.4699,-70.435)./(1-exp(complex(-5.5289e-5,4.8946e-6)).*zi);
  f6 = 12.946./(1 - 1.*zi);

  % Controller delay
  s = 1i*2*pi*f;
  dt = 315/1000;
  D = (exp(-s.*(T+dt)) .* (-1+exp(s.*T))  )./(s.*T);

  % Combine all terms into a single response
  H = 0.1.*( f1 + f2 + f3 + f4 + f5 + f6 ) .* D;

  % Build a frequency-series data object
  fsd = fsdata(f, H);
  
  % build the AO and set its name
  o = ao(fsd);
  % Set X-units, Y-units and name to something meaningful
  o.setXunits('Hz');
  o.setYunits('N/V');
  o.setName('C_{df}');
end

%--------------------------------------------------------------------------
function plo = getDefaultPlist(varargin)
  % List of available parameter sets
  sets = {'List', 'Range'};

  if nargin == 0
    plo = sets;
    return
  end

  set = varargin{1};

  switch set
    case 'List'
      plo = plist('Controller', 'DF', 'f', [1]);
    case 'Range'
      plo = plist('Controller', 'DF', ...
                  'f1', 1e-6,...
                  'f2', 5,...
                  'nf', 1000,...                
                  'scale', 'log');
    otherwise
    plo = plist();
  end
end
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'CLASS', '', 'CATEGORY', '$Id: ltpda_mdc1_C.m,v 1.3 2008/08/08 13:35:23 anneke Exp $', sets, pl);
end