% Noise budget script for FSS
%2017/06/10 by M.Nakano
%
clear all
close all
%%  add path
svnpath = '/Users/Masayuki/Documents/01.Study/05.svn/';
working_dir = [svnpath 'noisebudget/FSS/bKAGRA_phase1/'];
cd(working_dir)
datapath = [working_dir 'data'];
functionpath = [working_dir 'functions'];
ttfsspath = [svnpath 'ioo/TTFSS'];
addpath(datapath);
addpath(functionpath);
addpath(ttfsspath);
make_TFfile;
load('TFandparameters.mat');
model_name = 'FSS_NB';
%% import noises and transfer functions
import_data

%% compensation of the error signal for the SR560 used in measurement
SR560gain = 100;% the gain of SR560 used in measurement
err = err_compensation(err_no_compensation,SR560gain,ff_tf_OUT1C,abs_tf_OUT1C,ff_tf_OUT2C,abs_tf_OUT2C);

%% compensation of the olg for the gain between out1 and out2
abs_olg = olg_compensation(abs_olg_without_OUT12OUT2_10dB,ff_tf_OUT1C,abs_tf_OUT1C,ff_tf_OUT2C,abs_tf_OUT2C);

% mnbode(ff_olg,abs_olg,ang_olg)
%% cavity pole estimation
RefCav_L = 10e-2;
RefCav_cavitypole = estimate_cavitypole( ff_cavpole, abs_cavpole, ang_cavpole, false);
RefCav_finesse = 3e8/(4*RefCav_L*RefCav_cavitypole);
RefCav_r = (sqrt(pi^2+4*RefCav_finesse^2)-pi)/2/RefCav_finesse;
RefCav_T = 0.000093189;%from JGW-T1503493
RefCav_L = 1-RefCav_r^2-RefCav_T;

%% NPRO PZT actuator efficiency
NPRO_PZT = 51.75e6/35.0;%klog2853

%% optical gain estimation  @ error signal : 86mVpp
exc_peak = 10^(-52.073/20);
err_peak = 10^(-49.748/20);
RefCav = make_tfRefCav(exc_peak,err_peak,NPRO_PZT,RefCav_cavitypole);

%% EOM efficiency estimation from the measurement on 2017/06/06
EOM_eff = estimate_EOMeff(ff_cross,abs_cross,RefCav,NPRO_PZT,ttfsspath);

%% model servo
ss_ttfss = TTFSS_ss('/Users/Masayuki/Dropbox/share ioo/document/circuit/TTFSS',...
    'CG',334,'FG',270,'RL','Local','RT','Off','AO','Off',...
    'Exc','Off','Servo','Inv','FastExc','Off','FastSign','-');
ss_ttfss(2) = ss_ttfss(2)*5;%SR560

%% noise budget config
ff_NB = ff_err;
% noise spectrum
n_TTFSSSLOW = zeros(size(ff_NB));
n_TTFSSFAST = 0*make_ttfss_fast_noise(ff_NB,ff_s_noise,s_noise,ttfsspath,false);
n_TTFSSEOM = zeros(size(ff_NB));
n_RefCav = zeros(size(ff_NB));
n_RFPD = zeros(size(ff_NB));
n_NPRO = 10000./ff_NB;

% actuator efficiency
TF_NPRO_temp = -3e9;
TF_NPRO_PZT = NPRO_PZT;
TF_BEOM = EOM_eff*tf(1,[1/(100e6*2*pi) 1]);

% optical gain
TF_RefCav = RefCav*2.3; 
% 測定されたoptical gainだとこの日のopen loopと合わない
% 調査が必要


FSSSW = 1;
%% Open loop gain plot
load_system(model_name)
close all
[~,modelss_olg] = olgFromSimulink(model_name,ff_olg,'FSSOLGall');
[model_olg,~]=mnbode_list(modelss_olg,ff_olg,'c');
dt = 4e-7;
model_abs_olg = abs(model_olg);
model_ang_olg = angle(model_olg.*exp(-1i*dt*ff_olg*2*pi))*180/pi;
mnbode(ff_olg,abs_olg,ang_olg,ff_olg,model_abs_olg,model_ang_olg,'Title','Open loop gain','Legend',{'measured','modeled'})

%%
% Compute noises and save cache
% Compute nois  es
FSSerr = err;
dof = 'FSSERRcal';    % name of DOF to plot NB
[noiseserrFSS, sys] = nbFromSimulink(model_name, ff_NB, 'dof', dof);

% save cached outputs
saveFunctionCache();

%
% Make a sum noise
disp('Plotting noises')
nberrFSS = nbGroupNoises(model_name, noisesFSS, sys);%%
% Compute noises and save cache
% Compute nois  es
FSSfb = fb;
dof = 'FSSFBcal';    % name of DOF to plot NB
[noisesfbFSS, sys] = nbFromSimulink(model_name, ff_NB, 'dof', dof);

% save cached outputs
saveFunctionCache();

%
% Make a sum noise
disp('Plotting noises')
nbfbFSS = nbGroupNoises(model_name, noisesFSS, sys);

%%
close all
loglog(ff_NB,10000./ff_NB,noisesfbFSS{1}.f,noisesfbFSS{1}.asd,noiseserrFSS{1}.f,noiseserrFSS{1}.asd)
grid on
xlabel('frequency[Hz]');
ylabel('free run noise[Hz/rtHz]')
title('Free run noise')
legend('typical NPRO noise','estimated by fb signal','estimated by err signal')

%% inloop noise