% LTPDA_MDC1_IFO_NOISE returns a model spectrum of the IFO sensing noise for MDC1.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % DESCRIPTION: LTPDA_MDC1_IFO_NOISE returns a model spectrum of the IFO
    %              sensing noise for MDC1.
    %
    % CALL:        b = ltpda_mdc1_ifo_noise(pl)
    %
    % PARAMETERS:
    % 
    %           'source'   - Choose source of noise [default: 'ifo']
    %                        'ifo' - IFO sensing noise [m^2/Hz]
    %                        'sc'  - Spacecraft force noise [m^4s^{-2}/Hz]
    %                        'tm'  - Test-mass force noise [m^4s^{-2}/Hz]
    %           '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_input_noises.m,v 1.2 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_ifo_noise(ao, 'Params')
    %
    % The following call returns a string that contains the routine CVS version:
    %
    % >> version = ltpda_mdc1_ifo_noise(ao,'Version')
    %
    % The following call returns a string that contains the routine category:
    %
    % >> category = ltpda_mdc1_ifo_noise(ao,'Category')
    %
    % HISTORY: 11-04-08 M Hewitson
    %             Creation
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function varargout = ltpda_mdc1_input_noises(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

    source = find(pl, 'source');

    %% Compute response for frequencies f

    switch source
      case 'ifo'
        a = getIfoNoise(f);
      case 'tm'
        a = getTMNoise(f);
      case 'sc'
        a = getSCNoise(f);
      otherwise
        error('### Unknown noise source requested.');
    end

    varargout{1} = a;
    end
    %--------------------------------------------------------------------------
    % Get IFO noise for each frequency
    function o = getIfoNoise(f)

      % Description of the noise
      fcn = '0.5*(5e-12)^2 * ( 1 + (1+1e8)./(1+f.^2/1e-12) .* (1 + (1+9e6)./(1+f.^2/1e-14) ))';
      % Pack in a plist
      pl  = plist('fsfcn', fcn, 'f', f);

      % Build AO from plist
      o = ao(pl);
      o.setName('o_nxx', 'internal');
      o.setXunits('Hz');
      o.setYunits('m^2/Hz');
    end
    %--------------------------------------------------------------------------
    % Get SC noise for each frequency
    function o = getSCNoise(f)

      % Zeros as constant factors
      z1 = 1 + (10e-3/1e-6)^2;
      z2 = 1 + (0.1e-3/0.1e-6)^2;
      % AOs for frequency-dependent parts
      p1 = ao(plist('fsfcn', '1 + (f./100).^2', 'f', f));
      p2 = ao(plist('fsfcn', '1 + (f/1e-6).^2', 'f', f));
      p3 = ao(plist('fsfcn', '1 + (f/0.1e-6).^2', 'f', f));

      % Combine the parts
      G = (0.1e-6/436)^2;
      o = G .* (1./p1 + (z1./p2).*(1 + z2./p3 ));

      % Set properties of final AO
      o.setName('A_nxx', 'internal');
      o.setXunits('Hz');
      o.setYunits('m^2s^{-4}/ Hz');

    end
    %--------------------------------------------------------------------------
    % Get TM noise for each frequency
    function o = getTMNoise(f)

      % Construct the various parts of the spectrum

      % zeros
      z1 = ao(plist('fsfcn', '(1 + (f./100).^2)', 'f', f));
      z2 = (1 + (3e-3/1e-6).^2);
      z3 = ao(plist('fsfcn', '(1 + (f./1e-6).^2)', 'f', f));
      z4 = (1 + (0.1e-3./0.1e-6).^2);
      z5 = ao(plist('fsfcn', '(1 + (f./0.1e-6).^2)', 'f', f));

      % poles
      p1 = (1 + (1e-3./100).^2);
      p2 = (1 + (3e-3./1e-6).^2);
      p3 = (1 + (1e-3./1e-6).^2);
      p4 = (1 + (0.1e-3./0.1e-6).^2);
      p5 = (1 + (1e-3./0.1e-6).^2);

      % Combine the parts
      o = (30e-15)^2 .* ( (1./z1 + z2./z3.*(1 + z4./z5))  ./ (1./p1 + p2./p3.*(1 + p4./p5)));

      % Set properties of final AO
      o.setName('A_1xx', 'internal');
      o.setXunits('Hz');
      o.setYunits('m^2s^{-4}/ Hz');

    end
    %% Default parameters
    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('source', 'ifo', 'f', [1]);
        case 'Range'
          plo = plist('source', 'ifo', ...
                      '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_input_noises.m,v 1.2 2008/08/08 13:35:23 anneke Exp $', sets, pl);
    end