printFigs = false;
f1=transpose([0.001 : 0.0001 : 1]);
fetchdata = false;
sampf = 512; % sampling frequency in Hz

% Define and fetch Channels
if fetchdata
    Channels = {'H1:NGN-CS_CBRS_RY_PZT1_MASTER_OUT_DQ',...
                'H1:NGN-CS_CBRS_RY_PZT1_CTRL_IN1_DQ',...
                'H1:NGN-CS_CBRS_RY_PZT2_MASTER_OUT_DQ',...
                'H1:NGN-CS_CBRS_RY_PZT2_CTRL_IN1_DQ',...
                'H1:ISI-GND_STS_ITMY_Z_DQ',...
                'H1:NGN-CS_L4C_Z_5_OUT_DQ',...
                'H1:NGN-CS_L4C_Z_7_OUT_DQ',...
            };

    Starttime =  1154969000; % TF measurement
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    duration = 1600;      
%Fetch data
    buffers = conn.fetch(Starttime, Starttime+duration, Channels);  

    rawPZT1 = decimate(double(buffers(1).getData()), buffers(1).getChannel.getSampleRate/sampf);
    rawPDIN1 = decimate(double(buffers(2).getData()), buffers(2).getChannel.getSampleRate/sampf);
    rawPZT2 = decimate(double(buffers(3).getData()), buffers(3).getChannel.getSampleRate/sampf);
    rawPDIN2 = decimate(double(buffers(4).getData()), buffers(4).getChannel.getSampleRate/sampf);
    rawSTSZ = decimate(double(buffers(5).getData()), buffers(5).getChannel.getSampleRate/sampf);
    rawL4C5 = decimate(double(buffers(6).getData()), buffers(6).getChannel.getSampleRate/sampf);
    rawL4C7 = decimate(double(buffers(7).getData()), buffers(7).getChannel.getSampleRate/sampf);
    
    mydata = [rawPZT1 rawPDIN1 rawPZT2 rawPDIN2 rawSTSZ rawL4C5 rawL4C7];
    save('GPS1154969000_cBRS_Quiet.mat','mydata')   % Saving data for future reference
end    

%% File Saving & Retrieving  (Comment out if not needed)

if (not(exist('mydata','var')))
    if exist('GPS1154969000_cBRS_Quiet.mat','file')
        myfile = load('GPS1154969000_cBRS_Quiet.mat');
        mydata = myfile.mydata;
        rawPZT1 = mydata(:,1);
        rawPDIN1 = mydata(:,2);
        rawPZT2 = mydata(:,3);
        rawPDIN2 = mydata(:,4);
        rawSTSZ = mydata(:,5);
        rawL4C5 = mydata(:,6);
        rawL4C7 = mydata(:,7);
    end
end

%% Time selection

Sttime = 01 * sampf;
Totaltime = 1590; %duraiton in sec
Endtime = Sttime + Totaltime * sampf;

%% Calibrations

PZT1 = 1.1e-7*5/0.15*10/(131072) * rawPZT1(Sttime:Endtime);  % volts to radians
PDIN1 = 1.6e-10 * rawPDIN1(Sttime:Endtime);  % volts to radians
PZT2 = 1.1e-7*5/0.15*10/(131072) * rawPZT2(Sttime:Endtime);  % volts to radians
PDIN2 = 1.23e-10 * rawPDIN2(Sttime:Endtime);  % volts to radians
STSZ = 1e-9 * rawSTSZ(Sttime:Endtime);
L4C5 = 1.7e-7 * rawL4C5(Sttime:Endtime);
L4C7 = 1.7e-7 * rawL4C7(Sttime:Endtime);

%% Calculate ASDs
Rot_time = transpose(1/sampf * (0:1:length(PZT1)-1));
Navg = 33;

PDSUM = (PDIN1+PDIN2)/2;
PDDIFF = (PDIN1-PDIN2)/2;

[APZT1, ~] = asd2(PZT1,1/sampf, Navg, 1, @hann);
[APZT2, ~] = asd2(PZT2,1/sampf, Navg, 1, @hann);
[APD1, ~] = asd2(PDIN1,1/sampf, Navg, 1, @hann);
[APD2, freq] = asd2(PDIN2,1/sampf, Navg, 1, @hann);
[ASUM, ~] = asd2(PDSUM,1/sampf, Navg, 1, @hann);
[ADIFF, ~] = asd2(PDDIFF,1/sampf, Navg, 1, @hann);
[ASTSZ, ~] = asd2(STSZ,1/sampf, Navg, 1, @hann);
[AL4C5, ~] = asd2(L4C5,1/sampf, Navg, 1, @hann);
[AL4C7, ~] = asd2(L4C7,1/sampf, Navg, 1, @hann);

%% Calculate Coherences
Navg2 = 29;

[CohPZTs, ~] = coh2(PZT1,PZT2,1/sampf, Navg2, 1, @hann);
[CohPDs, freq2] = coh2(PDIN1,PDIN2,1/sampf, Navg2, 1, @hann);
[CohPDSUM, ~] = coh2(PDIN1,PDSUM,1/sampf, Navg2, 1, @hann);
[CohPDDIFF, ~] = coh2(PDIN1,PDDIFF,1/sampf, Navg2, 1, @hann);
[CohL4C5DIFF, ~] = coh2(PDDIFF,L4C5,1/sampf, Navg2, 1, @hann);
[CohSTSZDIFF, ~] = coh2(PDDIFF,STSZ,1/sampf, Navg2, 1, @hann);

%% Plots

figure(1)
l=plot(Rot_time, detrend(-PZT1,'constant'), Rot_time, detrend(PDIN1,'constant'), Rot_time, detrend(-PZT2,'constant'), Rot_time, detrend(PDIN2,'constant'));
set(l,'LineWidth',2);
set(gca,'FontSize',18)
xlabel('time (s)');
ylabel('PZT signal (rad)');
legend('PZT 1','PD 1','PZT 2','PD 2','Location','northeast');
grid on

figure(2)
ll=loglog(freq, APZT1,freq, APD1,freq, APZT2,freq, APD2, freq, ASTSZ);
set(ll,'LineWidth',2);
set(gca,'FontSize',18)
xlim([1e-3 200])
ylim([1e-12 1e-5])
%set(gca,'YTick',10.^(-10:-5))
xlabel('frequency (Hz)');
ylabel('ASD (rad/rt(Hz) or m/s/rt(Hz))');
legend('PZT 1','PD 1','PZT 2','PD 2','ITMY STS Z','Location','northeast');
grid on

figure(3)
    subplot(311)
    lll=semilogx(freq2, CohPDs,freq2, CohPZTs,'r');
    set(lll,'LineWidth',2);
    set(gca,'FontSize',18)
    xlim([1e-2 200])
    %ylim([1e-10 1e-5])
    %set(gca,'YTick',10.^(-10:-5))
    xlabel('frequency (Hz)');
    ylabel('Coherence');
    legend('PZT 1 and PZT 2','PD 1 and PD 2','Location','northeast');
    grid on

    subplot(312)
    lll=semilogx(freq2, CohPDSUM, freq2, CohPDDIFF,'r');
    set(lll,'LineWidth',2);
    set(gca,'FontSize',18)
    xlim([1e-2 200])
    %ylim([1e-10 1e-5])
    %set(gca,'YTick',10.^(-10:-5))
    xlabel('frequency (Hz)');
    ylabel('Coherence');
    legend('PD1 and PDSUM','PD1 and PDDIFF','Location','northeast');
    grid on

    subplot(313)
    lll=semilogx(freq2, CohL4C5DIFF, freq2, CohSTSZDIFF,'r');
    set(lll,'LineWidth',2);
    set(gca,'FontSize',18)
    xlim([1e-2 200])
    %ylim([1e-10 1e-5])
    %set(gca,'YTick',10.^(-10:-5))
    xlabel('frequency (Hz)');
    ylabel('Coherence');
    legend('L4C and PDDIFF','STSZ and PDDIFF','Location','northeast');
    grid on

figure(4)
%llll=loglog(freq, ASUM,freq, ADIFF,freq, APD1,freq, APD2);
llll=loglog(freq, ASUM,freq, ADIFF,'r', freq, ASTSZ, freq, AL4C5);
set(llll,'LineWidth',2);
set(gca,'FontSize',18)
xlim([1e-1 250])
ylim([1e-12 1e-6])
%set(gca,'YTick',10.^(-10:-5))
xlabel('frequency (Hz)');
ylabel('ASD (rad/rt(Hz) or m/s/rt(Hz))');
legend('PDSUM','PDDIFF (Actual floor tilt) ','ITMY STS Z','L4C 5','Location','northeast');
grid on

%% printing

if printFigs
%     figure(1)
%     FillPage('w')
%     saveas(gcf,'cBRS2_raw_angle.pdf')
        
%     figure(2)
%     FillPage('w')
%     saveas(gcf,'cBRS2_raw_ASD_Night.pdf')
%     
%     figure(3)
%     FillPage('w')
%     saveas(gcf,'cBRS2_Coherences.pdf')
%         
%     figure(4)
%     FillPage('w')
%     saveas(gcf,'cBRS_Sensitivity2.pdf')
end