% LTPDA_MDC1_IFO2ACC_FS calculates the external acceleration in the frequency-domain.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LTPDA_MDC1_IFO2ACC_FS calculates the external acceleration in 
%              the frequency-domain for MDC1.
%
% It computes: a = [DS^-1 + C]o
% 
% CALL:        b = ltpda_mdc1_ifo2acc_fd(pl)
%
% PARAMETERS:
% 
%           'Omega1'   - the square of the stiffness term for the dynamical response
%                        of test-mass 1 coupling to SC [default: 1.3e-6]
%           'Omega3'   - the square of the stiffness term for the dynamical response
%                        of test-mass 2 coupling to SC [default: 2e-6]
%           'delta'    - the cross-coupling factor of o1 into o12 [default: -1e-4]
%           'o1xx'     - spectral estimate of the IFO output o1  [default: empty ao]
%           'o12xx'    - spectral estimate of the IFO output o12 [default: empty ao]
%           'o112xx'   - cross-spectral estimate of the IFO output o1 and o12 [default: empty ao]
% 
% VERSION:     $Id: ltpda_mdc1_ifo2acc_fd.m,v 1.4 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_ifo2acc_fd(ao, 'Params')
%
% The following call returns a string that contains the routine CVS version:
%
% >> version = ltpda_mdc1_ifo2acc_fd(ao,'Version')
%
% The following call returns a string that contains the routine category:
%
% >> category = ltpda_mdc1_ifo2acc_fd(ao,'Category')
%
% HISTORY: 11-04-08 M Hewitson
%             Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = ltpda_mdc1_ifo2acc_fd(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);

%% Extract parameters from plist

% Get parameters out
w1     = find(pl, 'Omega1');
w3     = find(pl, 'Omega3');
delta  = find(pl, 'delta');
o1xx   = find(pl, 'o1xx');
o12xx  = find(pl, 'o12xx');
o112xx = find(pl, 'o112xx');

if isempty(o1xx.data) || isempty(o12xx.data) || isempty(o112xx.data)
  error('### Please provide three input spectra.');
end
if ~isequal(o1xx.data.x, o12xx.data.x, o112xx.data.x)
  error('### The two input spectra should be computed at the same frequencies.');
end

% Get frequency vector from one of the input spectra
f = o1xx.data.x;

%% Compute response for frequencies f

[a11xx, a22xx] = computeAcc(f, w1, w3, delta, o1xx, o12xx, o112xx);

varargout{1} = a11xx;
varargout{2} = a22xx;
end
%--------------------------------------------------------------------------
% Get DF controller for each frequency
function [a11xx, a22xx] = computeAcc(f, w1, w3, delta, o1xx, o12xx, o112xx)

% A model of the downsampling from 100 to 10Hz
pl = plist('gain',1, 'poles', [pz(10), pz(10), pz(10)], 'zeros', []);
pzm = pzmodel(pl);
ds  = resp(pzm, plist('f', f));

% Drag-free
Cdf  = ltpda_mdc1_C(plist('f', f, 'Controller', 'df'));
Adf  = ltpda_mdc1_actuator(plist('f', f, 'Actuator', 'df'));
Cdf  = ds.*Cdf.*Adf;
Sw1  = ltpda_mdc1_dynamics(plist('f', f, 'Omega2', w1.data.y.^2));
Sw1.setName('Sw1', 'internal');

% Suspension
Csus = ltpda_mdc1_C(plist('f', f, 'Controller', 'sus'));
Asus = ltpda_mdc1_actuator(plist('f', f, 'Actuator', 'sus'));
Csus = ds.*Csus.*Asus;
Sw3  = ltpda_mdc1_dynamics(plist('f', f, 'Omega2', w3.data.y.^2));
Sw3.setName('Sw3', 'internal');

% Square terms for dealing with PSDs
Csusxx = abs(Csus).^2;
Cdfxx  = abs(Cdf).^2;
Sw1xx  = abs(Sw1).^2;
Sw3xx  = abs(Sw3).^2;

wd   = w3.^2-w1.^2;
beta = wd - delta.*Sw3;
bxx  = abs(beta).^2;

% Calibrate back to a1

a11xx = o1xx .* (Sw1xx + Cdfxx - Cdf.*conj(Sw1) - conj(Cdf).*Sw1);
a11xx.setName('PSD(a11)', 'internal');

% Calibrate back to a2

e1 = o1xx.*bxx;
e2 = sqrt(o12xx).*Sw3;
e3  = sqrt(o12xx).*Csus;
a22t1 = e1 + abs(e2 - e3).^2;

a22t2 = beta.*(conj(Sw3)+conj(Csus)).*o112xx;
a22t3 = conj(beta).*(Sw3+Csus).*o112xx;
a22xx =  abs(a22t1 ...
       + a22t2  ...
       + a22t3);
a22xx.setName('PSD(a22)', 'internal');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Local Functions                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getInfo
%
% DESCRIPTION: Get Info Object
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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_ifo2acc_fd.m,v 1.4 2008/08/08 13:35:23 anneke Exp $', sets, pl);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getDefaultPlist
%
% DESCRIPTION: Get Default Plist
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plo = getDefaultPlist()
  plo = plist('Omega1', 1.3e-6, ...
            'Omega3', 2e-6, ...
            'delta', -1e-4, ...
            'o1xx', ao, ...
            'o12xx', ao, ...
            'o112xx', ao);
end