% test csd for ao/noisegen1D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% DESCRIPTION: Output a model for mdc1 cross spectral density
% 
% CALL:
%       [PSD,a1,a2,a3,a4,a5,av,am,plstd] = get_test_obj_ao_noisegen1D()
% 
%
% OUTPUTS:
%              CSD     -  a 2x2 matrix of analysis objects containing CSD
%                       models
%              TF      -  a 2x2 matrix of analysis objects containing TF
%                       models
%              a#      -  are white noise aos
%              av      -  a vector of wn aos
%              am      -  a matrix of wn aos
%              plst    -  a standard plist to call noisegen2D
% 
%
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [PSD,a1,a2,a3,a4,a5,av,am,plstd,plstd2] = get_test_obj_ao_noisegen1D()
  
  nsecs = 100;
  fs = 10; % Hz
  f = logspace(-6,log10(fs/2),100).'; % vector of frequencies Hz

  % Models response calculation
  tf11 = calcTF11(f,fs);
  tf12 = calcTF12(f,fs);

  % PSD calculation
  PSD = tf11.*conj(tf11)+tf12.*conj(tf12);
  PSD.setName('PSD');
  
  
  % Building test objects
  a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
  a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs));
  a3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'V'));
  a4 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', ''));
  a5 = ao(plist('tsfcn', 'randn(size(t)).''', 'fs', fs, 'nsecs', nsecs, 'yunits', 'fm s^-2'));
  av = [a1 a2 a3];
  am = [a1 a2; a3 a4];
  
  a1.setName;
  a2.setName;
  a3.setName;
  a4.setName;
  a5.setName;
  
  % building standard plist
  plstd = plist(...
      'model', PSD, ...
      'MaxIter', 60, ...
      'PoleType', 2, ...
      'MinOrder', 15, ...
      'MaxOrder', 45, ...
      'Weights', 3, ...
      'Plot', false,...
      'Disp', false,...
      'MSEVARTOL', 0.1,...
      'FITTOL', 0.01);
    
    plstd2 = plist(...
      'fs', fs, ...
      'Iunits', '', ...
      'Ounits', 'm', ...
      'MaxIter', 60, ...
      'PoleType', 2, ...
      'MinOrder', 15, ...
      'MaxOrder', 45, ...
      'Weights', 3, ...
      'Plot', false,...
      'Disp', false,...
      'MSEVARTOL', 0.1,...
      'FITTOL', 0.01);

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOCAL FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf11 = calcTF11(f,fs)
  
  
  
  % Model Stefano TF11 coefficients
  dRes11 = [2.44554138162509e-011 - 1.79482547894083e-011i;
  2.44554138162509e-011 + 1.79482547894083e-011i;
  2.66402334803101e-009 +  1.1025122049153e-009i;
  2.66402334803101e-009 -  1.1025122049153e-009i;
  -7.3560293387644e-009;
  -1.82811618589835e-009 - 1.21803627800855e-009i;
  -1.82811618589835e-009 + 1.21803627800855e-009i;
  1.16258677367555e-009;
  1.65216557639319e-016;                         
  -1.78092396888606e-016;
  -2.80420398962379e-017;
  9.21305973049041e-013 - 8.24686706827269e-014i;
  9.21305973049041e-013 + 8.24686706827269e-014i;
  5.10730060739905e-010 - 3.76571756625722e-011i;
  5.10730060739905e-010 + 3.76571756625722e-011i;
  3.45893698149735e-009;
  3.98139182134446e-014 - 8.25503935419059e-014i;
  3.98139182134446e-014 + 8.25503935419059e-014i;
  -1.40595719147164e-011];

  dPoles11 = [0.843464045655194 - 0.0959986292915475i;
  0.843464045655194 + 0.0959986292915475i;
  0.953187595424927 - 0.0190043625473383i;
  0.953187595424927 + 0.0190043625473383i;
  0.967176277937188;
  0.995012027005247 - 0.00268322602801729i;
  0.995012027005247 + 0.00268322602801729i;
  0.996564761885673;
  0.999999366165445;
  0.999981722418555;
  0.999921882627659;
  0.999624431675213 - 0.000813407848742761i;
  0.999624431675213 + 0.000813407848742761i;
  0.997312006278751 - 0.00265611346834941i;
  0.997312006278751 + 0.00265611346834941i;
  0.990516544257531;
  0.477796923118318 - 0.311064085401834i;
  0.477796923118318 + 0.311064085401834i;
  0];

  dDTerms11 = 0;
  
  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes11;
  pfparams.pol = dPoles11;
  pfparams.dterm = dDTerms11;
  pfr = utils.math.pfresp(pfparams);
  mtf11 = pfr.resp;
  
  % building AOs
  tf11 = ao(plist('xvals', f, 'yvals', mtf11, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf11.setName('TF11');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf12 = calcTF12(f,fs)
  
  % Model Stefano TF12
  dRes12 = [1.44258422208796e-017 + 7.07359428613009e-019i;
  1.44258422208796e-017 - 7.07359428613009e-019i;
  -3.4918408053655e-021 - 1.05662874569329e-021i;
  -3.4918408053655e-021 + 1.05662874569329e-021i;
  -7.61773292876976e-021;
  4.84357724603939e-020 + 2.38824204294595e-019i;
  4.84357724603939e-020 - 2.38824204294595e-019i;
  -4.07088520945753e-020 - 2.31474543846105e-019i;
  -4.07088520945753e-020 + 2.31474543846105e-019i;
  8.73316588658882e-023;
  -5.21840635377469e-020;
  1.8461911504859e-023;                         
  5.20105247464461e-020;
  -4.68960092394415e-022;
  -1.44261407664171e-017 +  6.8922564526833e-019i;
  -1.44261407664171e-017 -  6.8922564526833e-019i;
  3.13688133935426e-022];

  dPoles12 = [0.477546340377332 -     0.310830571032376i;
            0.477546340377332 +     0.310830571032376i;
             0.99790715414307 -    0.0028490561287024i;
             0.99790715414307 +    0.0028490561287024i;
            0.998014205354671                         ;
            0.999585354543332 -  0.000780408757425194i;
            0.999585354543332 +  0.000780408757425194i;
             0.99966003029931 -  0.000830944038363768i;
             0.99966003029931 +  0.000830944038363768i;
            0.999962770401331                         ;
            0.999981881865521                         ;
            0.999999365763457                         ;
            0.999981706320212                         ;
             0.99992421574188                         ;
            0.477898460791003 +     0.311001926610074i;
            0.477898460791003 -     0.311001926610074i;
                            0];
  dDTerms12 = 0;
  
  % response calculation
  pfparams.type = 'disc';
  pfparams.freq = f;
  pfparams.fs = fs;
  pfparams.res = dRes12;
  pfparams.pol = dPoles12;
  pfparams.dterm = dDTerms12;
  pfr = utils.math.pfresp(pfparams);
  mtf12 = pfr.resp;
  
  % building AOs
  tf12 = ao(plist('xvals', f, 'yvals', mtf12, 'fs', fs, 'type', 'fsdata'));
  % name for this object
  tf12.setName('TF12');

end