%% Preparation
run /users/swinkels/deploy/MatlabVirgoTools/trunk/startup.m

phase_ch = 'V1:SQZ_B1_7MHz_phi0';

ffl = '/virgoData/ffl/raw.ffl';

%% Parameters
% Previous attempts, used as tests
% t_start = utc2gps('2019-01-03 15:43:00');
% t_dur = 450;
% 
% t_start = utc2gps('2019-01-03 16:53:00');
% t_dur = 1500;
% 
% t_start = utc2gps('2019-01-03 20:08:00');
% t_dur = 1000;

% Modulation accidentally left on
ttext = 'Alignement modulation on';
t_base = [1230678907 1230678907+240]-3600; %interval to get baseline value before opening shutter
t_start = 1230679630-3600;
t_dur = 1230677150-t_start;

% % Modulation off
% ttext = 'Alignement modulation off';
% t_base = [utc2gps('2019-01-04 22:47:15') utc2gps('2019-01-04 22:48:30')]; %interval to get baseline value before opening shutter
% t_start = utc2gps('2019-01-04 22:49:05');
% t_dur = utc2gps('2019-01-04 23:10:30')-utc2gps('2019-01-04 22:49:05');
% %shutter closed at 23:10:40 and left quiet for at least 5 minutes

%TODO: currently can't change this because script relies on phase data
%vector being same length as BLRMS
BLRMSlen = 1;
BLRMSstep = 1;
BLRMSfl = 2600;
BLRMSfh = 2800;


%% Load data    
[phasedata, phasetime] = frgetvectN(ffl, phase_ch, t_base(1), t_start-t_base(1)+t_dur);
[BLRMStime, BLRMSdata] = fastBLRMS('V1:LSC_DARM', 10000, t_base(1), t_start-t_base(1)+t_dur, BLRMSfl, BLRMSfh, BLRMSlen, BLRMSstep, 0, 1);

[rangedata, rangetime] = frgetvectN(ffl, 'V1:Hrec_Range_BNS', t_base(1), t_start-t_base(1)+t_dur);
%% Analyze
timeidx = phasetime >= t_start-t_base(1) & BLRMStime <= t_start+t_dur-t_base(1);
phasevals = unique(phasedata(timeidx)); % Unique values of phase for which to computer averages
phasedata2 = ([0; phasedata(1:end-1)] + phasedata + [phasedata(2:end); 0])/3; %averaged phase data to remove transition points from computation
BLRMSmean = zeros(length(phasevals),1);
BLRMSstd = zeros(length(phasevals),1);
BLRMSstdm = zeros(length(phasevals),1);
for kk=1:length(phasevals)
    BLRMSvals = BLRMSdata(phasedata==phasevals(kk) & timeidx);
    BLRMSmean(kk) = sqrt(mean(BLRMSvals.^2)); % the mean is computed on the square of the BLRMS
    BLRMSstd(kk) = 1/(2*BLRMSmean(kk))*std(BLRMSvals.^2); %standard deviation of the square BLRMS, propagated through the sqrt
    BLRMSstdm(kk) = BLRMSstd(kk)/sqrt(length(BLRMSvals)); %standard deviation of the mean
end

baseBLRMSidx = BLRMStime <= diff(t_base);
baseBLRMSmean = sqrt(mean(BLRMSdata(baseBLRMSidx).^2));
baseBLRMSstd = 1/(2*baseBLRMSmean)*std(BLRMSdata(baseBLRMSidx).^2);
baseBLRMSstdm = baseBLRMSstd/sqrt(length(BLRMSdata(baseBLRMSidx)));

baserangeidx = rangetime <= diff(t_base) & rangedata~=0;
baserangemean = mean(rangedata(baserangeidx));
baserangestd = std(rangedata(baserangeidx));
baserangestdm = baserangestd/sqrt(length(rangedata(baserangeidx)));

%% Plots: phase and BLRM time series
figure
sb1 = subplot(3,1,1);
plot(phasetime,phasedata*180/pi,'.')
grid on
ylabel(['CC phase [deg]'])
title([ttext ' - GPS time: ' num2str(t_start)])
set(gca,'YTick',(0:45:225))

sb2 = subplot(3,1,2);
plot(BLRMStime, BLRMSdata)
hold all
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*baseBLRMSmean,'r-')
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*(baseBLRMSmean+baseBLRMSstd),'r--')
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*(baseBLRMSmean-baseBLRMSstd),'r--')
ylims=get(gca,'ylim');
plot([1 1],ylims,'r--','linewidth',3)
plot(diff(t_base)*[1 1],ylims,'r--','linewidth',3)
grid on
ylabel(['DARM BLRMS (' num2str(BLRMSfl) '-' num2str(BLRMSfh) ' Hz)'])
legend({'Measured data', 'Reference with shutter closed (\pm 1\sigma)'})

sb3 = subplot(3,1,3);
rangeidx = rangedata~=0;
plot(rangetime(rangeidx), rangedata(rangeidx))
hold all
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*baserangemean,'r-')
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*(baserangemean+baserangestd),'r--')
plot([BLRMStime(1), BLRMStime(end)],[1, 1]*(baserangemean-baserangestd),'r--')
ylims=get(gca,'ylim');
plot([1 1],ylims,'r--','linewidth',3)
plot(diff(t_base)*[1 1],ylims,'r--','linewidth',3)
grid on
ylabel('BNS range')
xlabel('Time [s]')


linkaxes([sb1, sb2, sb3],'x')

%% Plots: BLRMS vs phase
figure
errorbar(phasevals*180/pi, BLRMSmean, BLRMSstdm)
grid on
xlabel('CC phase [deg]')
set(gca,'XTick',(0:45:225))
ylabel(['DARM BLRMS (' num2str(BLRMSfl) '-' num2str(BLRMSfh) ' Hz)'])
title([ttext ' - GPS time: ' num2str(t_start)])

%% Plots: BLRMS vs phase in DB
figure
errorbar(phasevals*180/pi, 20*log10(BLRMSmean/baseBLRMSmean), 20/log(10)*(BLRMSstdm./BLRMSmean+baseBLRMSstdm./baseBLRMSmean))
grid on
xlabel('CC phase [deg]')
set(gca,'XTick',(0:45:225))
ylabel(['DARM BLRMS (' num2str(BLRMSfl) '-' num2str(BLRMSfh) ' Hz) ration wrt reference [dB]'])
title([ttext ' - GPS time: ' num2str(t_start)])
