% A test script for ao/dispersion
% 
% M Nofrarias 07-05-09
% 
% $Id: test_ao_dispersion.m,v 1.2 2009/05/08 13:15:36 miquel Exp $
% 

function test_ao_dispersion()

fs =10;
f = logspace(-5,log10(fs/2),100)';

%% nominal plist
pl = plist( 'f', f,...
'fs', fs,...
'wp1',0.0010488*1i,...  % sqrt(-11e-7)
'wp2',0.0014832*1i,...  % sqrt(-22e-7)
'Cdf',0.8,...  
'Clfs',1.15,...  
'dS',1.35e-4,...
'p_ifo1',[3.6e-12 10e-3 4.2 1.8e-3 8],...
'p_tm1',[7e-15 5e-3 3 4e-4 9],...
'p_ifo2', [3.6e-12 10e-3 4.2 1.8e-3 8],...
'p_tm2', [7e-15 5e-3 3 4e-4 9],...
'p_sc', [2.5e-10 12e-3 3.8 1e-3 8],...
'nsecs', 1e4,...
'fsgen', 50);

%% Compute transfer function

tf1 = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',1),pl);
% tf2 = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',3),pl);
% tf3 = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',4),pl);

% matched stiffness case
% 
% plx = pset(pl,param('wp1',0.00154919*i)) % sqrt(-24e-7)
% plx = pset(plx,param('wp2',0.00154919*i)) % sqrt(-24e-7)
% 
% tf1x = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',1),plx);
% tf2x = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',3),plx);
% tf3x = ao(plist('built-in','mdc2r2_fd_ifo2ifo','index',4),plx);


a = ao(fsdata(f,(f'.^2).*ones(1,length(f))));
a2 = ao(fsdata(f,(1./f'.^2).*ones(1,length(f))));

%% dispersion function

disp = dispersion(a,plist('model','mdc2r2','channel','i1o1'),pl);
disp2 = dispersion(a2,plist('model','mdc2r2','channel','i1o1'),pl);


%% plot

iplot(abs(tf1),disp,disp2)

%% compute transfer function error

noise_in = length(f)*fs/2*noise_param(f,7e-15,5e-3,3,4e-4,9)
noise_out = length(f)*fs/2*noise_param(f,3.6e-12,10e-3,4.2,1.8e-3,8)


varTF = disp.*(noise_in.*abs(tf1) + noise_out)
varTF.setName
varTF2 = disp2.*(noise_in.*abs(tf1) + noise_out)
varTF2.setName

% iplot(abs(tf1),plist('YErrU',sqrt(varTF)))
iplot(abs(tf1),sqrt(varTF),sqrt(varTF2))

end


%--------------------------------------------------------------------------
% Parametric function for noise contributions
%--------------------------------------------------------------------------

function o = noise_param(varargin)
f = varargin{1};
p1 = varargin{2};
p2 = varargin{3};
p3 = varargin{4};
p4 = varargin{5};
p5 = varargin{6};

% Description of the noise
fcn = p1^2*(1 + 1./(f./p2).^p3 + 1./(f./p4).^p5).^(1/2) ;

% Build AO from plist
o = ao(fsdata(f,fcn));


end