% Noise budget  PDHGR test script
%
% $Id$

%% Add paths
clear all
close all
addpath('../');
findNbSVNroot;
addpath(genpath([NbSVNroot 'Common/Utils']));
addpath(genpath([NbSVNroot 'Dev/Utils/']));

%% Define some filters and plot expected OLTF

%%%%%%%parameter%%%%%%%%%
c = 299792458;
L = 3000;
Finesse = 10;
lambda = 1.064e-6;
m2Hz = 2/lambda * c/2/L;
%%%%%%%%%%%%%%%%%%%%%%%%%

SHG_1 = 2;
Cavity_Pole = zpk([],[-2*pi*c/4/L/Finesse],2*pi*c/4/L/Finesse);
Sen = Finesse/(c/2/L);
Filt_GR = zpk([-2*pi*2.5e3],[0],0.1);

% RF
Filt_GR1 = zpk([0],[-10*2*pi],1);
Act_AOM = -1e6;
AOM = 1;

% Low frequency
Filt_mass = zpk([-10*2*pi,-10*2*pi,-10*2*pi],[0,0,-100*2*pi,-500*2*pi,-500*2*pi],-200*500*2*pi*500*2*pi);
num = [(-5.0e-6)*2*m2Hz*2*pi*2*pi];
den = [1 2*pi/10 4*pi*pi];
Act_mass = tf(num,den);


%Cavity
Changing_unit = c/L/532e-9;
num1 = [0.01*pi*pi];
den1 = [1 0.1*pi 0.01*pi*pi];

num2 = [4*pi*pi];
den2 = [1 2*pi/10 4*pi*pi];

inverted_pend = tf(num1,den1);
Pendulum =  tf(num2,den2);

Gol = Cavity_Pole*Sen*Filt_GR*(Filt_GR1*Act_AOM*AOM-Filt_mass*Act_mass);

%%%% to create a bode plot with the freqeuncy axis in Hz
set(cstprefs.tbxprefs,'FrequencyUnits','Hz')

figure(100)
bode(Gol)
grid on

%% Define some parameters and get live parts parameters
freq = logspace(-12,4,1000);
liveModel = 'TEST_PDH';
dof = 'TST';    % name of DOF to plot NB
 startTime = 1078250000;   % start GPS time
 durationTime = 512;
 IFO = 'H1';
 site = 'LHO';

% Try setting different NDS server if you couldn't get data
% setenv('LIGONDSIP','h1nds1:8088');
% mdv_config;

% load cached outputs
loadFunctionCache()

% get live parts parameters
%liveParts(liveModel, startTime, durationTime, freq)

%% Compute noises and save cache
% Compute noises
[noises, sys] = nbFromSimulink(liveModel, freq, 'dof', dof);

% save cached outputs
saveFunctionCache();

%% Make a quick NB plot
disp('Plotting noises')
nb = nbGroupNoises(liveModel, noises, sys);

% Get noise data from DAQ. Put NdNoiseSource block with DAQ channel
% specified. Put something (e.g. 1) in ASD parameter of that block.
%nb = nbAcquireData(liveModel, sys, nb, startTime, durationTime);

nb.sortModel();
matlabNoisePlot(nb);
figure(1)
%ylim([1e-10,1])

%% plot expected curve from calculation

%[mag,ph]=bode((Gol)/(1+Gol),2*pi*freq);
%loglog(freq,squeeze(mag),'b.')
%hold on

%[mag,ph]=bode(1/Sen,2*pi*freq);
%loglog(freq,freq./freq/Sen,'r.')