%% time period
gps0 = 1122454860 + 1320;
gps1 = gps0 + 503*60;

%% get range
conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
data = conn.fetch(gps0, gps1, {'H1:DMT-SNSW_EFFECTIVE_RANGE_MPC.mean,m-trend'});
range = data(1).getData();
fs = data(1).getLength() / (gps1-gps0);
t = gps0 + (0:length(range)-1)/fs;

%% find glitches
% detrend
np = 60;
wind = hanning(np);
wind = wind / sum(wind);
ranget = conv([ones(np,1) * range(1) ; range; ones(np,1) * range(end)], wind, 'same');
ranget = ranget(np+1:end-np);

% select largest drops
idx = find(range>0 & (range-ranget)<-6);

f1 = figure();
plot(t - gps0, range, t(idx) - gps0, range(idx), 'or')
xlabel(sprintf('Time starting from GPS %d [s]', gps0)gps0 = 1122454860 + 1320;
gps1 = gps0 + 503*60;
)
ylabel('Range [Mpc]')
legend('Range', 'Glitches')gps0 = 1122454860 + 1320;
gps1 = gps0 + 503*60;

grid
img = getframe(f1);
imwrite(img.cdata, 'range_drops.png')

fprintf(2, 'Found %d glitches\n', length(idx));
pause

%% get three minutes of data around each drop
g = GWData;
for i=1:length(idx)
    % get DARM data
    [e, tt] = g.fetch(t(idx(i))- 60, t(idx(i)) + 120, {'H1:LSC-DARM_IN1_DQ'});
    fs = length(e)/180;

    % compute spectra 
    np = 10 * fs;
    [sp1, fr] = pwelch(e(1:fs*60), hanning(np), np/2, np, fs);
    [sp2, fr] = pwelch(e(fs*60+1:fs*120), hanning(np), np/2, np, fs);
    [sp3, fr] = pwelch(e(fs*120+1:fs*180), hanning(np), np/2, np, fs);
    
    % plot them
    scrsz = get(0,'ScreenSize');
    f2 = figure('Position',[1 scrsz(4)/1.8 scrsz(3)/2 scrsz(4)/2.4]);
    subplot(211)
    loglog(fr, sqrt(sp1), fr, sqrt(sp2), fr, sqrt(sp3))
    title(sprintf('Glitch #%d GPS %d', i, t(idx(i))))
    legend('60 s before', '60s during', '60s after')
    grid
    xlabel('Frequency [Hz]')
    ylabel('DARM\_IN1 [cts]')
    xlim([10 2000])
    ylim([1e-14, 1e-11])
    
    % find the exact time using maximum of BLRMS
    [B,A] = butter(4, [30, 300]/(fs/2), 'bandpass');
    ee = filtfilt(B, A, double(e));
    ee = ee.^2;
    [~, j] = max(ee(fs:end-fs));
    j = j + fs;
    
    % GPS time of the glitch
    tglitch = t(idx(i)) - 60 + tt(j);
    
    % plot time series (+- deltaT seconds)
    deltaT = 2;
    subplot(212)
    if j-fs*deltaT > 0 
        plot(tt(j-fs*deltaT:j+fs*deltaT), e(j-fs*deltaT:j+fs*deltaT), tt(j), e(j), 'or')
    else
        plot(tt(1:fs*2*deltaT), e(1:fs*2*deltaT), tt(j), e(j), 'or')
    end
    title(sprintf('Glitch #%d GPS %d', i, t(idx(i))))
    xlabel('Time [s]')
    ylabel('DARM\_IN1 [cts]')
    grid
    img = getframe(f2);
    imwrite(img.cdata, sprintf('glitch_%d_darm.png', i))

    % and check saturations (ETMX)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-ETMX*MASTER*DQ');
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);
    
    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    tshift = round(tglitch) - t(idx(i)) - deltaT + 60;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(5,5,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_etmx.png', i))

    % and check saturations (ETMY)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-ETMY*MASTER*DQ');
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);
    
    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(5,5,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_etmy.png', i))

    
    % and check saturations (ITMX)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-ITMX*MASTER*DQ');
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);
    
    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(5,5,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_itmx.png', i))    
    
    
    % and check saturations (ITMY)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-ITMY*MASTER*DQ');
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);
    
    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(5,5,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_itmy.png', i))      
    
    % and check saturations (PRM)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-PRM*MASTER*DQ');     
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);

    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(4,6,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_prm.png', i))

    % and check saturations (PR2)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-PR2*MASTER*DQ');     
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);

    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(4,6,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_pr2.png', i))    
    
    % and check saturations (SRM)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-SRM*MASTER*DQ');     
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);

    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(4,4,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_srm.png', i))
    
    % and check saturations (SR2)
    conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);
    ch = conn.findChannels('H1:SUS-SR2*MASTER*DQ');     
    channels = arrayfun(@(x) x.getName(), ch, 'UniformOutput', false);

    [data,ta] = g.fetch(round(tglitch)-deltaT, round(tglitch)+deltaT, channels);
    fsa = size(data,1)/10;
    
    scrsz = get(0,'ScreenSize');
    f3 = figure('Position',[1 0 scrsz(3)/2 scrsz(4)/2.35]);
    for k=1:size(data,2)
        subplot(4,6,k)
        
        ii = find(ta > tglitch - round(tglitch) + deltaT, 1);
        plot(ta + tshift, data(:,k), ta(ii) + tshift, data(ii,k), 'ro')
        title(strrep(channels{k}, '_', '\_'))
    end

    img = getframe(f3);
    imwrite(img.cdata, sprintf('glitch_%d_sr2.png', i))       
    
    close all
end
