
close all
clear all
%% get data
addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/
addpath(genpath('/ligo/apps/linux-x86_64/mDV/extra')) %path for tconvert

%clear all
quads={'ITMX'};
quad_names=['ITMX'];
%H1.channels={};
for ii=1:length(quads)
    chan_names.M0(ii)=strcat('H1:SUS-',quads(ii),'_M0_DAMP_P_INMON.mean,m-trend');
    chan_names.R0(ii)=strcat('H1:SUS-',quads(ii),'_R0_DAMP_P_INMON.mean,m-trend');
    chan_names.RF1(ii)=strcat('H1:SUS-',quads(ii),'_R0_COILOUTF_F1_OUT16.mean,m-trend');
    chan_names.L1(ii)=strcat('H1:SUS-',quads(ii),'_L1_WIT_PMON.mean,m-trend');
    chan_names.L2(ii)=strcat('H1:SUS-',quads(ii),'_L2_WIT_PMON.mean,m-trend');
    chan_names.L3(ii)=strcat('H1:SUS-', quads(ii),'_L3_OPLEV_PIT_OUT16.mean,m-trend');
    chan_names.F1(ii)=strcat('H1:SUS-',quads(ii),'_M0_COILOUTF_F1_OUT16.mean,m-trend');
    chan_names.ISI_state(ii)=strcat('H1:ISI-',quads(ii),'_ST2_ISO_X_GAIN.mean,m-trend');
    chan_names.BIAS(ii)=strcat('H1:SUS-',quads(ii),'_M0_OPTICALIGN_P_OFFSET.mean,m-trend');
    %chan_names.OL_sum(ii)=strcat('H1:SUS-', quads(ii),'_L3_OPLEV_SUM_OUT16.mean,m-trend');
end

start_time=tconvert('29 July 2017 00:00:00');
stop_time=tconvert('2 August 2017 10:00:00');
dt=stop_time-start_time;

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

chans=[];
stage_names=fieldnames(chan_names);
for ii=1:length(stage_names)
    chans=[chans, chan_names.(stage_names{ii})];
end
%add a few channels that aren't sus
chans=[chans,'H1:ISI-GND_STS_ETMX_Z_BLRMS_30M_100M.mean,m-trend'];

 tic
 [slow_data, t, info] = gwd.fetch(start_time, dt, chans);
 toc
 
 %get the data into a structure that is convient, there must be a better
 %way
data.M0=slow_data(:,1);
data.R0=slow_data(:,2);
data.RF1=slow_data(:,3);
data.L1=slow_data(:,4);
data.L2=slow_data(:,5);
data.L3=slow_data(:,6);
data.F1=slow_data(:,7);
data.ISI=slow_data(:,8);
data.BIAS=slow_data(:,9);

sei_data=slow_data(:,10);

%calibrate F1 counts into urad using alignment slider and Eul2Osem matrix
F1_cnts2urad=[-1/(12.82*15.589), -1/(12.82*19.92), -1/(12.82*15.344), -1/(12.82*19.341)];
RF1_cnts2urad=[-1/(12.82*23.219), -1/(12.82*23.219), -1/(12.82*23.219), -1/(12.82*23.2191)];

%make it easier to plot all data on same (similar) scales, , calibrate forces into urad
for ii=1:1
    %calbrate applied forces into urad using bias sliders
    data.F1(:,ii)=data.F1(:,ii).*F1_cnts2urad(ii);
    data.RF1(:,ii)=data.RF1(:,ii).*RF1_cnts2urad(ii);
    %detrend some channels (not oplevs since there is time when there is no light on diode
%and numbers are huge)
    data.M0(:,ii)=data.M0(:,ii)-mean(data.M0(:,ii));
    data.R0(:,ii)=data.R0(:,ii)-mean(data.R0(:,ii));
    data.L1(:,ii)=data.L1(:,ii)-mean(data.L1(:,ii));
    data.L2(:,ii)=data.L2(:,ii)-mean(data.L2(:,ii));
    data.F1(:,ii)=data.F1(:,ii)-mean(data.F1(:,ii));
    data.RF1(:,ii)=data.RF1(:,ii)-mean(data.RF1(:,ii));
    %data.L3(:,ii)=data.L3(:,ii)-mean(data.L3(:,ii));

end

%top mass diff is M0-R0
top_mass_diff=data.M0-data.R0; %how well does this predict L1?
%put time into minutes (only do this once each time you download data)
t=t./(60);

%% make plots
EQ_thresh=3e2; %level of ground motion that we use to recognize EQ time
EQ_indx=5000;
post_EQ=(t>t(EQ_indx)); %vector which is 0 before EQ, 1 after



for quad=1:1
    ISI_iso(:,quad)=(data.ISI(:,quad)==1);
    iso_post(:,quad)=post_EQ.*ISI_iso(:,quad);
    reiso(quad)=find(iso_post,1);%index of first time ISI was isolated after the EQ
end
    %realign=find(t>1.11e4, 1);
    reiso(quad)=5000;
    realign=5500;

figure(1)
title('ITMX Status')
subplot(5,1,1)
    semilogy(t, sei_data)
    grid on
    legend('ITMX Z BLRMS 30M 100M')
subplot(5,1,2)
    plot(t, iso_post(:,1))
    legend('ISI Isolation Status')
    grid on
subplot(5,1,3)
    plot(t,data.L3(:,1),'k')
    legend( 'Optical Lever')%'all','before EQ', 'during', 'after',
    grid on
subplot(5,1,4)
    plot(t,data.M0(:,1))
    legend('M0 osem readbacks')
    grid on
subplot(5,1,5)
    plot(t,data.F1(:,1))
    legend('F1 Force')
    grid on
    xlabel('Time (Minutes) from 2017-07-29 00:00:00')

oplev_upper_lim=[20,20,20,20];
oplev_lower_lim=[-20,-20,-20,-20];

stage_str={'M0','L1','L2','F1'};
for quad=1:1
    fig=figure(10+quad);
    for ii=1:4
        subplot(2,2,ii)
            plot(data.L3(1:EQ_indx,quad),data.(stage_str{ii})(1:EQ_indx,quad), '.r',...
                data.L3(reiso(quad):realign,quad),data.(stage_str{ii})(reiso(quad):realign,quad), '.g',...
                data.L3(realign:end,quad),data.(stage_str{ii})(realign:end,quad), '.b')
            xlabel('Op Lev Pit (urad)')
            xlim([oplev_lower_lim(ii), oplev_upper_lim(ii)])
            ylim([-20,20])
            ylabel(strcat(stage_str{ii},' Osems (urad)'))
            grid on
    end
    subplot(2,2,3)
    legend('Before Test', 'Approx During Test', 'After Test')
    subplot(2,2,4)
    ylabel('pit force (M0 F1 in test mass urad)')
    ylim([-400,400])
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,strcat('\bf ',quads(quad)),'HorizontalAlignment','center','VerticalAlignment', 'top')
end


for quad=1:1
    figure(40+quad)
    title(quads(quad))
        subplot(4,1,1)
            plot(data.M0(1:EQ_indx,quad),data.F1(1:EQ_indx,quad), '.r',...
                data.M0(reiso(quad):realign,quad),data.F1(reiso(quad):realign,quad), '.g',...
                data.M0(realign:end,quad),data.F1(realign:end,quad), '.b')
            xlabel('M0 Osems readback (urad)')
            %xlim([-500,500])
            ylabel('M0 Osems force (urad)')
            %ylim([-600,600])
            xlim([-80,80])
            ylim([-25,25])
            grid on
            legend('Before Test', 'Approx During Test', 'After Test')
        subplot(4,1,2)
            plot(data.L1(1:EQ_indx,quad),data.M0(1:EQ_indx,quad), '.r',...
                data.L1(reiso(quad):realign,quad),data.M0(reiso(quad):realign,quad), '.g',...
                data.L1(realign+1:end,quad),data.M0(realign+1:end,quad), '.b')
            xlabel('L1 Osems (urad)')
            %xlim([-500,500])
            ylabel('M0 Osems (urad)')
            ylim([-1000,1000])
            grid on
        subplot(4,1,3)
            plot(data.L2(1:EQ_indx,quad),data.L1(1:EQ_indx,quad), '.r',...
                data.L2(reiso(quad):realign,quad),data.L1(reiso(quad):realign,quad), '.g',...
                data.L2(realign+1:end,quad),data.L1(realign+1:end,quad), '.b')
            xlabel('L2 Osems (urad)')
            %xlim([-250, 250])
            ylabel('L1 Osems (urad)')
            %ylim([-500,500])
            grid on
        subplot(4,1,4)
            plot(data.L3(1:EQ_indx,quad),data.L2(1:EQ_indx,quad), '.r',...
                data.L3(reiso(quad):realign,quad),data.L2(reiso(quad):realign,quad), '.g',...
                data.L3(realign+1:end,quad),data.L2(realign+1:end,quad), '.b')
            xlabel('Op Lev Pit (urad)')
            xlim([-2,5])
            ylabel('L2 Osems (urad)')
            ylim([-3,3])
            grid on
            
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,strcat('\bf ',quads(quad)),'HorizontalAlignment','center','VerticalAlignment', 'top')
    
end

figure(20)
for quad=1:1
    subplot(1,1,quad)
            plot(data.M0(1:EQ_indx,quad),data.F1(1:EQ_indx,quad), '.r',...
                data.M0(reiso(quad):realign,quad),data.F1(reiso(quad):realign,quad), '.g',...
                data.M0(realign:end,quad),data.F1(realign:end,quad), '.b')
            xlabel('M0 Osems readback (urad)')
            xlim([-60,60])
            ylabel('M0 Osems force (urad)')
            ylim([-40,40])
            grid on
            title(quads(quad))
end
    legend('Before Test', 'Approx During Test', 'After Test')
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,'\bf Linear relationship between top mass force and top mass position offset by EQ','HorizontalAlignment','center','VerticalAlignment', 'top')

    
print('-f20','M0_Force_vs_position','-dpdf')

figure(21)
for quad=1:1
    subplot(1,1,quad)
            plot(data.R0(1:EQ_indx,quad),data.RF1(1:EQ_indx,quad), '.r',...
                data.R0(reiso(quad):realign,quad),data.RF1(reiso(quad):realign,quad), '.g',...
                data.R0(realign:end,quad),data.RF1(realign:end,quad), '.b')
            xlabel('R0 Osems readback (urad)')
            xlim([-1,1])
            ylabel('R0 Osems force (urad)')
            ylim([-.1,.1])
            grid on
            title(quads(quad))
end

    legend('Before Test', 'Approx During Test', 'After Test')
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,'\bf R0 force vs position','HorizontalAlignment','center','VerticalAlignment', 'top')
    
print('-f21','R0_Force_vs_position','-dpdf')

figure(22)
for quad=1:1
    subplot(1,1,quad)
            plot(top_mass_diff(1:EQ_indx,quad),data.L1(1:EQ_indx,quad), '.r',...
                top_mass_diff(reiso(quad):realign,quad),data.L1(reiso(quad):realign,quad), '.g',...
                top_mass_diff(realign:end,quad),data.L1(realign:end,quad), '.b')
            xlabel('R0-M0 Osems readback (urad)')
            %xlabel('M0 force (urad)')
            xlim([-300,300])
            ylabel('L1 Osems (urad)')
            ylim([-300,300])
            grid on
            title(quads(quad))
            legend('Before Test', 'Approx During Test', 'After Test')
            
end
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    %text(0.5, 1,'\bf Linear relationship between top mass force and top mass position mostly unchanged by EQ','HorizontalAlignment','center','VerticalAlignment', 'top')
    
print('-f22','M0R0_vs_L1','-dpdf')

figure(23)
for quad=1:1
    subplot(1,1,quad)
            plot(data.M0(1:EQ_indx,quad),data.L1(1:EQ_indx,quad), '.r',...
                data.M0(reiso(quad):realign,quad),data.L1(reiso(quad):realign,quad), '.g',...
                data.M0(realign:end,quad),data.L1(realign:end,quad), '.b')
            xlabel('M0 Osems readback (urad)')
            xlim([-300,300])
            ylabel('L1 Osems (urad)')
            ylim([-300,300])
            grid on
            legend('Before Test', 'Approx During Test', 'After Test')
            title(quads(quad))
end
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    %text(0.5, 1,'\bf Linear relationship between top mass force and top mass position mostly unchanged by EQ','HorizontalAlignment','center','VerticalAlignment', 'top')
    
    
figure(24)
for quad=1:1
    subplot(1,1,quad)
            plot(data.L1(1:EQ_indx,quad),data.L2(1:EQ_indx,quad), '.r',...
                data.L1(reiso(quad):realign,quad),data.L2(reiso(quad):realign,quad), '.g',...
                data.L1(realign:end,quad),data.L2(realign:end,quad), '.b')
            xlabel('L1 Osems (urad)')
            xlim([-300,300])
            ylabel('L2 Osems (urad)')
            ylim([-300,300])
            grid on
            title(quads(quad))
end
    legend('Before Test', 'Approx During Test', 'After Test')
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    %text(0.5, 1,'\bf Linear relationship between top mass force and top mass position mostly unchanged by EQ','HorizontalAlignment','center','VerticalAlignment', 'top')

    
 figure(25)
for quad=1:1
    subplot(1,1,quad)
            plot(data.L2(1:EQ_indx,quad),data.L3(1:EQ_indx,quad), '.r',...
                data.L2(reiso(quad):realign,quad),data.L3(reiso(quad):realign,quad), '.g',...
                data.L2(realign:end,quad),data.L3(realign:end,quad), '.b')
            xlabel('L2 osems (urad)')
            xlim([-30,30])
            ylabel('oplevs (urad)')
            ylim([-30,30])
            grid on
            title(quads(quad))
end
    legend('Before Test', 'Approx During Test', 'After Test')
    ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    %text(0.5, 1,'\bf Linear relationship between top mass force and top mass position mostly unchanged by EQ','HorizontalAlignment','center','VerticalAlignment', 'top')

%% find shifts
before.M0F=[-545, 186,-64.45,237.2];
after.M0F=[-426,33,-7.426, 282];
before.M0R=[-79,280,-92,-21];
after.M0R=[88,-324,134,44];
before.L1=[-26,236,-113,38];
after.L1=[39,-262,134,-24];
before.L2=[-2,21,17,66];
after.L2=[9.3,-31,-13,-52];

shift.M0F=after.M0F-before.M0F;
shift.M0R=after.M0R-before.M0R;
shift.L2=after.L2-before.L2;
shift.L1=after.L1-before.L1;

%%
