% LTPDA_MDC1_DYNAMICS returns a frequency-domain model of the dynamics for MDC1.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_MDC1_DYNAMICS returns a frequency-domain model of the dynamics
%              for MDC1.
%
% CALL:        b = ltpda_mdc1_dynamics(pl)
%
% PARAMETERS:
% 
%           'Omega2'   - the square of the stiffness term for the dynamical response
%                        [default: 1.3e-6]
%           '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_dynamics.m,v 1.4 2008/08/12 12:33:25 anneke Exp $
%
% The following call returns a parameter list object that contains the
% default parameter values:
%
% >> pl = ltpda_mdc1_dynamics(ao, 'Params')
%
% The following call returns a string that contains the routine CVS version:
%
% >> version = ltpda_mdc1_dynamics(ao,'Version')
%
% The following call returns a string that contains the routine category:
%
% >> category = ltpda_mdc1_dynamics(ao,'Category')
%
% HISTORY: 11-04-08 M Hewitson
%             Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = ltpda_mdc1_dynamics(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

omega2 = find(pl, 'Omega2');

%% Compute response for frequencies f

varargout{1} = getDynamics(f, omega2);
end
%--------------------------------------------------------------------------
% Get DF controller for each frequency
function o = getDynamics(f, w2)

  % We take abs() here and use +- omega later
  w1 = sqrt(abs(w2));

  % DF dynamics
  pl  = plist('gain',(w1*w1), 'poles', [], 'zeros', [pz(w1/2/pi) pz(-w1/2/pi)]);
  dfm  = pzmodel(pl);

  % Make an AO from the pzmodel
  o    = resp(dfm, plist('f', f))  ;
  o.setName('s^2+\omega^2');
  o.addHistory(getInfo,plist('Omega2', w2),[],o.hist);
end

function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pls   = [];
  elseif nargin == 1&& ~isempty(varargin{1}) && ischar(varargin{1})
    sets{1} = varargin{1};
    pls = getDefaultPlist(sets{1});
  else
    sets = {'List', 'Range'};
    pls = [];
    for kk=1:numel(sets)
      pls = [pls getDefaultPlist(sets{kk})];
    end
  end
  % Build info object
  ii = minfo(mfilename, 'MDC1', '', 'Signal Processing', '$Id: ltpda_mdc1_dynamics.m,v 1.4 2008/08/12 12:33:25 anneke Exp $', sets, pls);
end

function plo = getDefaultPlist(set)
  switch set
    case 'List'
      plo = plist('Omega2', 1.3e-6, 'f', [1]);
    case 'Range'
      plo = plist('Omega2', 1.3e-6, ...
        'f1', 1e-6,...
        'f2', 5,...
        'nf', 1000,...
        'scale', 'log');
    otherwise
      plo = plist();
  end
end