run /users/swinkels/deploy/MatlabVirgoTools/trunk/startup.m

% g = 1147800377; % friday excited mirror
g = 1148055660; % monday

[dc, ac, vel, t] = getChannels( ...
    {'LSC_B7_DC', 'LSC_NArm_6MHz_ROT_CH', 'LSC_GL_vel_est'}, g, 200);

%% find peaks
peak = dc > 0.0012;
[istart, istop] = block_istart_istop(dc > 0.0014);
ipeak = nan(size(istart));

for i=1:length(istart)
    d = dc(istart(i):istop(i));
    [~, imx] = max(d);
    ipeak(i) = istart(i) + imx - 1;
end

%% hack: kick points that are too close

nkick = nnz(diff(ipeak) < 300);

for i=1:nkick
    [~, imn] = min(diff(ipeak)); % find closest spaced peaks
    if dc(ipeak(imn)) < dc(ipeak(imn+1))  % remove lowest of two peaks
        ipeak(imn) = [];
    else
        ipeak(imn + 1) = [];
    end
end
        
    
%%

% hack, beginning of ringing effect in ac signal at moment of peak is clear
% indication of sign of velocity
v_sign = sign(ac(ipeak));

% count fringes up and down
nfringe = cumsum([0; 0.5 * (v_sign(1:end-1) + v_sign(2:end))]);

% cubic spline interpolation 
fit = spline(t(ipeak), nfringe);
nfringe_int = ppval(fit, t);

% calculate cavity length and velocity
z_int = 1.064e-6 / 2 * nfringe_int;
v_int = [0; (z_int(3:end) - z_int(1:end-2)) / (t(3) - t(1)); 0];


%%
multiPlot({t, dc, t(ipeak), dc(ipeak), '.r'}, ...
    {t, 1e6*z_int, t(ipeak), 1e6*z_int(ipeak), '.r'}, ...
    {t, 1e6 * v_int}, ...
    'xlabel', sprintf('Time (s) after %s', gps2str(g)), ...
    'title', {'LSC\_B7\_DC', 'interpolated cavity length', 'estimated velocity'}, ...
    'ylabel', {'W', 'um', 'um/s'})

