%%
clear all;
close all;
%%
H1.channels = {'H1:SUS-ITMX_L3_OPLEV_PIT_OUTPUT', 'H1:SUS-ITMX_L3_OPLEV_YAW_OUTPUT', ...
    'H1:SUS-ITMX_L2_WIT_PMON', 'H1:SUS-ITMX_L2_WIT_YMON', ... 
    'H1:SUS-ITMX_M0_TEST_P_OUT16',  'H1:SUS-ITMX_M0_TEST_Y_OUT16', ...
    'H1:SUS-ITMX_M0_TEST_P_GAIN', 'H1:SUS-ITMX_M0_TEST_Y_GAIN', ...
    'H1:SUS-ITMX_R0_DAMP_V_INMON', ...
    'H1:SUS-ITMY_L3_OPLEV_PIT_OUTPUT', 'H1:SUS-ITMY_L3_OPLEV_YAW_OUTPUT', ...
    'H1:SUS-ITMY_L2_WIT_PMON', 'H1:SUS-ITMY_L2_WIT_YMON', ... 
    'H1:SUS-ITMY_M0_TEST_P_OUT16',  'H1:SUS-ITMY_M0_TEST_Y_OUT16', ...
    'H1:SUS-ITMY_M0_TEST_P_GAIN', 'H1:SUS-ITMY_M0_TEST_Y_GAIN', ...
    'H1:SUS-ITMY_M0_DAMP_V_INMON', ...
    'H1:CDS-SENSMON_CAL_SNSW_EFFECTIVE_RANGE_MPC'};
for ii=1:length(H1.channels)
    H1.channels{ii}=[H1.channels{ii}, '.mean,m-trend'];
end

RANGEIDX=19;
OPTICSIDX=[1, 10];
VOFF=8;
WITOFF=2;
TESTOFSTOFF=4;
TESTGAINOFF=6;
%%
if ~exist('mintrend.mat', 'file')
    gwd = GWData();
    gwd.make_kerberos_ready;
    gwd.site_info(2).port = 31200;
    gwd.site_info(2).server = 'nds.ligo-wa.caltech.edu';
    ts = tconvert('9 Jan 2017 00:00:00');
    te = tconvert('23 Jan 2017 00:00:00');
    [data, t, ~] = gwd.fetch(ts, te-ts, H1.channels);
    save('mintrend.mat', 'data', 't', 'H1');
else
    load('mintrend.mat');
end

%%
MPC=data(:, RANGEIDX);
OPTICSNAME={'ITMX', 'ITMY'};
PYNAME={'PIT', 'YAW'};
OLsensitivities=[];
WTsensitivities=[];
for OPTICS=1:2
    for PY=1:2
        OPLEVIDX=OPTICSIDX(OPTICS)+PY-1;
        OPLEV=data(:, OPLEVIDX);
        WIT=data(:, OPLEVIDX+WITOFF);
        TESTOFST=data(:, OPLEVIDX+TESTOFSTOFF)./data(:, OPLEVIDX+TESTGAINOFF);
        delta=TESTOFST(2:end)-TESTOFST(1:end-1); % change in the minute average trend in microrad
        isMoving=find(abs(delta)>1e-3); % changes larger than 1nrad.
        ii=1;
        movingSeg=[];
        while ii<length(isMoving)
            sIdx=isMoving(ii);
            eIdx=sIdx;
            % find continuously "moving" segments
            for jj=1:(length(isMoving)-ii)
                if sIdx+jj==isMoving(ii+jj)
                    eIdx=sIdx+jj;
                else
                    movingSeg=[movingSeg; sIdx, eIdx];
                    ii=ii+jj;
                    break;
                end
            end
            if ii+jj==length(isMoving)
                break
            end
        end
        % first row of movingSeg = before any change happened.
        % second row = last minute of consecutive motion.
        % Now look at oplev response to the TEST offset step.
        
        OLresp=[];
        WTresp=[];
        respT=[];
        for ii=1:length(movingSeg)
            OFSTSTEP=TESTOFST(movingSeg(ii, 2)+1)-TESTOFST(movingSeg(ii, 1));
            if OFSTSTEP==0
                continue;
            else
                respT=[respT; (t(movingSeg(ii, 2)+1)+t(movingSeg(ii, 1)))/2];
                OLSTEP=OPLEV(movingSeg(ii, 2)+1)-OPLEV(movingSeg(ii, 1));
                OLresp=[OLresp; OLSTEP/OFSTSTEP];
                WTSTEP=WIT(movingSeg(ii, 2)+1)-WIT(movingSeg(ii, 1));
                WTresp=[WTresp; WTSTEP/OFSTSTEP];
            end
        end
        OLsensitivities{(OPTICS-1)*2+PY}=[respT, OLresp];
        WTsensitivities{(OPTICS-1)*2+PY}=[respT, WTresp];

        figure
        plot(t/3600/24, TESTOFST, '.', t(movingSeg(:, 2)+1)/3600/24, TESTOFST(movingSeg(:, 2)+1), 'ro', t(movingSeg(:, 1))/3600/24, TESTOFST(movingSeg(:, 1)), 'go');
        grid on
        xlabel('days since Jan 09 2017 00:00 UTC')
        ylabel('urad')
        title([OPTICSNAME{OPTICS}, ' ', PYNAME{PY}, ' M0 TEST offset']);
        legend('mean minute trend', 'end of constant offset', 'start of constant offset')

        figure
        subplot 311
        plot(t/3600/24, TESTOFST)
        grid on
        legend('M0 TEST OFFSET')
        ylabel('urad')
        title([OPTICSNAME{OPTICS}, ' ', PYNAME{PY}, ' raw data']);
        subplot 312
        plot(t/3600/24, OPLEV)
        grid on
        legend('OPLEV')
        ylabel('urad')
        subplot 313
        plot(t/3600/24, WIT)
        grid on
        legend('L2 WIT')
        ylabel('urad')
        xlabel('days since Jan 09 2017 00:00 UTC')
        
        figure

        subplot 211
        plot(t/3600/24, MPC)
        grid on
        legend('BNS range')
        ylabel('Mpc')
        xlim([min(t), max(t)]/3600/24);
        subplot 212
        plot(respT/3600/24, OLresp, 'o');
        plot(respT/3600/24, WTresp, 'x')
        grid on
        xlim([min(t), max(t)]/3600/24);
        ylim([0.7, 1.2])
        title([OPTICSNAME{OPTICS}, ' sensitivity against M0 TEST ', PYNAME{PY}, ' offset']);
        legend('OPLEV');
        legend('L2 WIT')
        ylabel('urad/urad')
        xlabel('days since Jan 09 2017 00:00 UTC')
    end
end
%%
figure
subplot 411
plot(t/24/3600, MPC)
grid on
ylabel('Mpc')
subplot 412
plot(WTsensitivities{1}(:,1)/3600/24, WTsensitivities{1}(:,2), 'ro', ...
   WTsensitivities{2}(:,1)/3600/24,  WTsensitivities{2}(:,2), 'bo', ...
   WTsensitivities{3}(:,1)/3600/24,  WTsensitivities{3}(:,2), 'gx', ...
   WTsensitivities{4}(:,1)/3600/24,  WTsensitivities{4}(:,2), 'kx');
ylim([0.5, 1.2]);
grid on
legend('IX PIT', 'IX YAW', 'IY PIT', 'IY YAW');
ylabel('L2 WIT/M0 TEST, urad/urad')
subplot 413
plot(t/24/3600, data(:, OPTICSIDX(2)+VOFF), [t(1), t(end)]/24/3600, [-65.73, -65.73])
legend('IY M0 V position', 'threshold')
grid on
ylabel('um')
subplot 414
plot(t/24/3600, data(:, OPTICSIDX(1)+VOFF));
grid on
legend('IX R0 V position')
ylabel('um')
xlabel('days since Jan 09 2017 00:00 UTC')
      

figure
plot(WTsensitivities{1}(:,1)/3600/24, WTsensitivities{1}(:,2), 'ro', ...
   WTsensitivities{2}(:,1)/3600/24,  WTsensitivities{2}(:,2), 'bo', ...
   WTsensitivities{3}(:,1)/3600/24,  WTsensitivities{3}(:,2), 'gx', ...
   WTsensitivities{4}(:,1)/3600/24,  WTsensitivities{4}(:,2), 'kx');
ylim([0.8, 1.1]);
grid on
legend('IX PIT', 'IX YAW', 'IY PIT', 'IY YAW');
ylabel('L2 WIT/M0 TEST, urad/urad')
xlabel('days since Jan 09 2017 00:00 UTC')
