--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_dispersion.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,94 @@
+% 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