% load list of times
glitches = load('glitch_times.txt');
% select only class 1 and 2
glitches = glitches(glitches(:,2) == 1 | glitches(:,2) == 2, :);

%% get one minute of data around each drop (for some reason I can only access the first 97 glitches at Caltech, and the others at Hanford...)
conn = nds2.connection('nds.ligo.caltech.edu', 31200);    
for i=1:97
    fprintf(2, '%d\n', i);
    % get DARM data
    data = conn.fetch(glitches(i,1)-1, glitches(i,1)+1, {'H1:LSC-DARM_IN1_DQ'});
    e = data(1).getData();
    fs = length(e)/2;
    tt = (0:length(e)-1)/fs;
    
    % compute BLRMS
    [B,A] = butter(4, [30, 300]/(fs/2), 'bandpass');
    ee = filtfilt(B, A, double(e));
    ee = ee.^2;
    % find the right time
    [size_blrms(i), j] = max(ee(fs/8:end-fs/8));
    j = j + fs/8;
    
    % high pass
    [B,A] = butter(4, 50/(fs/2), 'high');
    eh = filtfilt(B, A, double(e));
    
    % notch some lines
    fnotch = [60, 502, 992, 1000, 1462, 2450];
    Qz = 1000;
    Qp = 10;
    for k=1:numel(fnotch)
        num = [1, 2*pi*fnotch(k)/Qz, (2*pi*fnotch(k))^2];
        den = [1, 2*pi*fnotch(k)/Qp, (2*pi*fnotch(k))^2];
        [B,A] = bilinear(num, den, fs, fnotch(k));
        eh = filtfilt(B,A,eh);
    end
    
    shape = @(p, t) (t>p(1) & t<p(1)+p(2)) .* (t-p(1))/p(2)*p(3) + (t>p(1)+p(2) & t<p(1)+p(2)+p(4)) .* (p(3) - (t-p(1)-p(2))/p(4)*p(3)) + ...
         - (t>p(1)+p(2)+p(4) & t<p(1)+p(2)+p(4)+p(5)) .* (t-p(1)-p(2)-p(4))/p(5)*p(6) + ...
         (t>p(1)+p(2)+p(4)+p(5) & t<p(1)+p(2)+p(4)+p(5)+p(7)) .* (-p(6) + (t-p(1)-p(2)-p(4)-p(5))/p(7)*p(6));

    % p(1) = start time
    % p(2) = rising time
    % p(3) = amplitude of first peak x 1e9
    % p(4) = decay time
    % p(5) = second rising time
    % p(6) = (negative) amplitude of second peak x 1e9
    % p(7) = second decay time
     
    % fit with shape (restrict to 0.1 seconds around glitch time
    idx = j-round(fs*0.05):j+round(fs*0.05);
    % initial guesses
    [~,jj] = max(abs(eh(idx)));
    peak1 = eh(idx(jj));
    if peak1 > 0
        [peak2,kk] = min(eh(idx));
    else
        [peak2,kk] = max(eh(idx));
    end
    p0 = [tt(j)-0.001, 0.001, 1e9*peak1, 0.001, 0.004, -1e9*peak2, 0.008];
    p(i,:) = lsqcurvefit(shape, p0, tt(idx), 1e9*eh(idx)');
    
    % plot
    scrsz = get(0,'ScreenSize');
    f2 = figure('Position',[1 scrsz(4)/2 scrsz(3)/3 scrsz(4)/1.2]);
    plot(tt(idx), 1e9*eh(idx), 'b', tt(idx), shape(p(i,:), tt(idx)), 'r', 'LineWidth', 2)
    xlabel('Time [s]')
    ylabel('DARM\_IN1 [cts]')
    title('High pass 100 Hz - notched')
    
    img = getframe(f2);
    imwrite(img.cdata, sprintf('glitch_plots/fit_%d_darm.png', i))
    pause
    close(f2)
    
    save amplitudes.mat size_blrms p 
end

%% get one minute of data around each drop (for some reason I can only access the first 97 glitches at Caltech, and the others at Hanford...)
conn = nds2.connection('nds.ligo-wa.caltech.edu', 31200);    
for i=98:size(glitches,1)
    fprintf(2, '%d\n', i);
    % get DARM data
    data = conn.fetch(glitches(i,1)-1, glitches(i,1)+1, {'H1:LSC-DARM_IN1_DQ'});
    e = data(1).getData();
    fs = length(e)/2;
    tt = (0:length(e)-1)/fs;
    
    % compute BLRMS
    [B,A] = butter(4, [30, 300]/(fs/2), 'bandpass');
    ee = filtfilt(B, A, double(e));
    ee = ee.^2;
    % find the right time
    [size_blrms(i), j] = max(ee(fs/8:end-fs/8));
    j = j + fs/8;
    
    % high pass
    [B,A] = butter(4, 50/(fs/2), 'high');
    eh = filtfilt(B, A, double(e));
    
    % notch some lines
    fnotch = [60, 502, 992, 1000, 1462, 2450];
    Qz = 1000;
    Qp = 10;
    for k=1:numel(fnotch)https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=
        num = [1, 2*pi*fnotch(k)/Qz, (2*pi*fnotch(k))^2];
        den = [1, 2*pi*fnotch(k)/Qp, (2*pi*fnotch(k))^2];
        [B,A] = bilinear(num, den, fs, fnotch(k));
        eh = filtfilt(B,A,eh);
    end
    
    shape = @(p, t) (t>p(1) & t<=p(1)+p(2)) .* (t-p(1))/p(2)*p(3) + (t>p(1)+p(2) & t<=p(1)+p(2)+p(4)) .* (p(3) - (t-p(1)-p(2))/p(4)*p(3)) + ...
         - (t>p(1)+p(2)+p(4) & t<=p(1)+p(2)+p(4)+p(5)) .* (t-p(1)-p(2)-p(4))/p(5)*p(6) + ...
         (t>p(1)+p(2)+p(4)+p(5) & t<=p(1)+p(2)+p(4)+p(5)+p(7)) .* (-p(6) + (t-p(1)-p(2)-p(4)-p(5))/p(7)*p(6));

    % p(1) = start time
    % p(2) = rising time
    % p(3) = amplitude of first peak x 1e9
    % p(4) = decay time
    % p(5) = second rising time
    % p(6) = (negative) amplitude of second peak x 1e9
    % p(7) = second decay time
     
    % fit with shape (restrict to 0.1 seconds around glitch time
    idx = j-round(fs*0.05):j+round(fs*0.05);
    % initial guesseshttps://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=
    [~,jj] = max(abs(eh(idx)));
    peak1 = eh(idx(jj));
    if peak1 > 0
        [peak2,kk] = min(eh(idx));
    else
        [peak2,kk] = max(eh(idx));
    end
    p0 = [tt(j)-0.001, 0.001, 1e9*peak1, 0.001, 0.004, -1e9*peak2, 0.008];
    p(i,:) = lsqcurvefit(shape, p0, tt(idx), 1e9*eh(idx)');
    
    % plot
    scrsz = get(0,'ScreenSize');
    f2 = figure('Position',[1 scrsz(4)/2 scrsz(3)/3 scrsz(4)/1.2]);
    plot(tt(idx), 1e9*eh(idx), 'b', tt(idx), shape(p(i,:), tt(idx)), 'r', 'LineWidth', 2)
    xlabel('Time [s]')
    ylabel('DARM\_IN1 [cts]')
    title('High pass 100 Hz - notched')
    
    img = getframe(f2);
    imwrite(img.cdata, sprintf('glitch_plots/fit_%d_darm.png', i))
    pause
    close(f2)
    
    save amplitudes.mat size_blrms p 
end


% glitches that are of dubious origin
db = [16, 44, 60, 85, 87, 95, 104, 117, 121, 123, 126, 130, 140];

% select only the good ones
idx = setdiff(1:size(glitches,1), db);

glitches = glitches(idx,:);
p = p(idx,:);

X(:,1) = glitches(:,1);   % gps time
X(:,2) = p(:,3);          % amplitude of the first peak
X(:,3) = abs(p(:,2));     % first peak rising time
X(:,4) = abs(p(:,4));     % first peak decay time
X(:,5) = p(:,6);          % amplitude of second peak
X(:,6) = abs(p(:,5));     % rising time of the second peak
X(:,7) = abs(p(:,7));     % decay time of the second peak

f = fopen('glitch_shapes.txt', 'w');
fprintf(f, '# GPS, FirstPeakAmplitude, FirstPeakRisingTime_ms, FirstPeakDecayTime_ms, SecondPeakAmplitude, SecondPeakRisingTime_ms, SecondPeakDecayTime_ms\n');
for i=1:size(X,1)
    fprintf(f, '%d %d %f %f %d %f %f\n', floor(X(i,1)), 1e-9*X(i,2), 1e3*X(i,3), 1e3*X(i,4), 1e-9*X(i,5), 1e3*X(i,6), 1e3*X(i,7));
end
fclose(f);


% histograms
x = -25:0.15:5;
[n,x] = hist(X(:,2), x);
plot(x,n, 'o-', 'LineWidth', 2)
grid
xlabel('Amplitude of the first peak in DARM\_IN1 (x1e-9)')
ylabel('Counts')


[n1,x1] = hist(X(:,3), 100);
[n2,x2] = hist(X(:,4), x1);
subplot(211)
plot(x1,n1, 'o-', 'LineWidth', 2)
grid
ylabel('Counts')
title('First peak rising time')
subplot(212)
plot(x2,n2, 'o-', 'LineWidth', 2)
grid
xlabel('Time [ms]')
ylabel('Counts')
title('First peak decay time')


[n,x] = hist(X(:,5), 100);
plot(x,n, 'o-', 'LineWidth', 2)
grid
xlabel('Amplitude of the second peak in DARM\_IN1 (x1e-9)')
ylabel('Counts')


[n1,x1] = hist(X(:,6), 100);
[n2,x2] = hist(X(:,7), x1);
subplot(211)
plot(x1,n1, 'o-', 'LineWidth', 2)
grid
ylabel('Counts')
title('Second peak rising time')
subplot(212)
plot(x2,n2, 'o-', 'LineWidth', 2)
grid
xlabel('Time [ms]')
ylabel('Counts')
title('Second peak decay time')

