% addpath /ligo/svncommon/NbSVN/aligonoisebudget/trunk/Externals/SimulinkNb/Utils/

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


%% Fetch the data in chunks, concatenate it for use.
channels = {'H1:OMC-PI_DCPD_64KHZ_AHF_DQ'}; 


time_excite         = tconvert('3 July 2019 16:38:00'); % Turn damping gain on at 16:38:00 UTC, off at 16:41:15
time_ringdown_start = tconvert('3 July 2019 16:41:15');
time_ringdown_end   = tconvert('3 July 2019 17:12:33'); % started power up at 17:12:33

pretime = 400; % start data taking a little beforehand, so can window / ignore it.
times_fetch = (time_excite-pretime):500:time_ringdown_end;

duration = 500; % 500 sec okay, 1000 sec not okay at one time.  Can concat them after fetching.

srate    = 65536;
DCPDdata = zeros(1,srate*duration*length(times_fetch));
times    = zeros(1,srate*duration*length(times_fetch));


for ii = 1:length(times_fetch)
    
    [data, t, ~] = gwd.fetch(times_fetch(ii), duration, channels);

    
    DCPDdata(1+(ii-1)*duration*srate:(ii)*duration*srate) = data;
    times(1+(ii-1)*duration*srate:(ii)*duration*srate)    = t+times_fetch(ii)-time_ringdown_start;
    
end



%% Use GWFilt to bandpass the data

BPdata = GWFilt.bandpass(DCPDdata, 10429.5, 10431.5, srate, 8);

clear data t DCPDdata % Get back some memory

% figure(1)
% plot(times, BPdata)


%% Try a Hilbert Transform

hilb = imag(hilbert(BPdata));

envelope = sqrt(BPdata.^2 + hilb.^2);

clear hilb % get back some memory

% Plot to confirm that the envelope matches the ringdown shape
% figure(3) 
% plot(times, BPdata,...
%     times, envelope)

%% Decimate the data

envelope_downsampled = decimate(envelope, 64);
times_downsampled = downsample(times, 64);

keepIDX = find(times_downsampled > 20 & times_downsampled < max(times)-100);


%%

fun = @(x,timepoints)x(1)*exp(x(2)*timepoints); % define my function

x0 = [4.8,-0.0039]; % initial guesses of parameters (from a fit of the log of the data)
x = lsqcurvefit(fun,x0,times_downsampled(keepIDX),envelope_downsampled(keepIDX));

%%
plot(times_downsampled,envelope_downsampled,...
    times_downsampled(keepIDX),fun(x,times_downsampled(keepIDX)))
ylim([0,150])
xlim([-100,max(times)])

xlabel('Time since excitation stopped [sec]')
ylabel('Magnitude of bandpassed data [arb]')
title('OMC DCPD 64kHz channel, bandpassed, downsampled')

legend('Envelope magnitude', ['Fit: ' num2str(x(1),3) '*exp(' num2str(x(2),3) '*t)'])

%%



plot(times, BPdata,...
    times_downsampled(keepIDX), fun(x,times_downsampled(keepIDX)),...
    times_downsampled(keepIDX), -fun(x,times_downsampled(keepIDX)))

ylim([-150,150])
xlim([-100,max(times)])

xlabel('Time since excitation stopped [sec]')
ylabel('Bandpassed OMC DCPD [arb]')
title('OMC DCPD 64kHz channel, bandpassed 2 Hz around 10430.5 Hz')

legend('Bandpassed DCPD', ['Fit: ' num2str(x(1),3) '*exp(' num2str(x(2),3) '*t)'], ['-1*Fit: -' num2str(x(1),3) '*exp(' num2str(x(2),3) '*t)'])

SaveFigureFunction('RingdownFit_10430Hz')
