addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/

gwd = GWData();
gwd.make_kerberos_ready;
gwd.site_info(2).port = 31200;
gwd.site_info(2).server = 'nds.ligo-wa.caltech.edu';


%%


time_fetch_SRCLinj = tconvert('17 Sept 2016 22:24:30');
duration = 60 * 4;

channels = {'H1:LSC-SRCL_OUT_DQ',...
            'H1:LSC-DARM_IN1_DQ'};

[data_SRCLinj, t_SRCLinj, ~] = gwd.fetch(time_fetch_SRCLinj, duration, channels);
srate_SRCLinj = 1/(t_SRCLinj(2) - t_SRCLinj(1));




%%
time_fetch_FFinj = tconvert('17 Sept 2016 22:13:20');
duration = 60 * 4;

% SRCLFF is only in Commissioning Frames, so must be accessed from local nds
channels = {'H1:LSC-SRCLFF1_OUT_DQ',...
            'H1:LSC-DARM_IN1_DQ'};

[data_FFinj, t_FFinj, ~] = gwd.fetch(time_fetch_FFinj, duration, channels);
srate_FFinj = 1/(t_FFinj(2) - t_FFinj(1));

%%

newSrate = min(srate_SRCLinj, srate_FFinj);

srcl_SRCLinj = detrend(decimate(data_SRCLinj(:,1),srate_SRCLinj/newSrate));
darm_SRCLinj = detrend(decimate(data_SRCLinj(:,2),srate_SRCLinj/newSrate));

srcl_FFinj = detrend(decimate(data_FFinj(:,1),srate_FFinj/newSrate));
darm_FFinj = detrend(decimate(data_FFinj(:,2),srate_FFinj/newSrate));


save('SRCL_FF_injection_data', 'srcl_SRCLinj', 'darm_SRCLinj', 'srcl_FFinj', 'darm_FFinj', 'newSrate' ,'time_fetch_SRCLinj', 'time_fetch_FFinj')

nfft = newSrate * 2;

%%


[srcspec, freq] = pwelch(srcl_SRCLinj, hanning(nfft), nfft/2, nfft,newSrate);
srcspec = sqrt(srcspec);
[darmspec, freq] = pwelch(darm_SRCLinj, hanning(nfft), nfft/2, nfft,newSrate);
darmspec = sqrt(darmspec);

[srcffspec, freq] = pwelch(srcl_FFinj, hanning(nfft), nfft/2, nfft,newSrate);
srcffspec = sqrt(srcffspec);
[darmffspec, freq] = pwelch(darm_FFinj, hanning(nfft), nfft/2, nfft,newSrate);
darmffspec = sqrt(darmffspec);

% figure(10)
% subplot(211)
% loglog(freq,darmspec./srcspec)
% xlim([10,200])
% subplot(212)
% loglog(freq,darmffspec./srcffspec)
% xlim([10,200])
% 
% figure(20)
% loglog(freq, (darmspec./srcspec)./(darmffspec./srcffspec))
% xlim([10,200])

%%

% "actuator" is darm_in1/srclff_out
% coupling is darm_in1/srcl_out

% ff filter should be coupling/actuator



%% 

[coupTF, freq] = tfestimate(srcl_SRCLinj, darm_SRCLinj, hanning(nfft), nfft/2, nfft,newSrate);
[ffTF, freq]   = tfestimate(srcl_FFinj,   darm_FFinj,   hanning(nfft), nfft/2, nfft,newSrate); 

ffFilt = coupTF ./ ffTF;

% figure(11)
% subplot(211)
% loglog(freq, abs(coupTF))
% xlim([10,200])
% subplot(212)
% loglog(freq, abs(ffTF))
% xlim([10,200])


figure(21)
subplot(211)
loglog(freq, abs(ffFilt))
xlim([10,200])

subplot(212)
semilogx(freq, 180/pi * angle(ffFilt))
xlim([10,200])


save('SRCL_FF_Filter', 'freq', 'ffFilt')









