% plot_data.m
%# Copyright (C) 2012  Alan Weinstein
%#
%# This program is free software; you can redistribute it and/or modify it
%# under the terms of the GNU General Public License as published by the
%# Free Software Foundation; either version 2 of the License, or (at your
%# option) any later version.
%#
%# This program is distributed in the hope that it will be useful, but
%# WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
%# Public License for more details.
%#
%# A copy of the GNU General Public License may be found at
%# http://www.gnu.org/copyleft/gpl.html
%# or write to the Free Software Foundation, Inc.,
%# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

% This is example MATLAB code to read in detector strain data,
% and plot it in various ways.

%clear all
%close all

% loop over detectors:
IFOs =  {'H2' 'L1'};
for I = 1:2,

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % First, the 256s of data around GRB051103
   % sampling frequency in Hz and duration in seconds
   fs = 4096;
   Duration = 256;
   GPSstart = 815045078;

   % open the file and get the data (strain time series):
   IFO = char(IFOs(I));
   fn = [IFO,'-','STRAIN_',num2str(fs),'Hz-',num2str(GPSstart),...
       '-',num2str(Duration),'.hdf5'];
   dataset = [IFO,':LSC-STRAIN'];
   data = hdf5read(fn,dataset);
   data = double(data);

   % print out some basic data checks
   display(['Contents of ',fn]) 
   display(['IFO, Npoints = ',IFO,', ',num2str(length(data))])
   display(['Duration = ',num2str(length(data)/(fs*Duration)), ...
       ' * ',num2str(fs),' * ',num2str(Duration)])
   display(['min,max,mean,std = ',num2str(min(data)),', ', ...
      num2str(max(data)),', ',num2str(mean(data)),', ',num2str(std(data))])

   %% histogram the raw data:
   [hn,hbins] = hist(data, 100);
   % make a corresponding Gaussian:
   N = max(hn);
   mu = mean(data);
   sigma = std(data);
   y = N*exp(-((hbins-mu)/sigma).^2/2.0);
   % plot it:
   figure
   bar(hbins,hn)
   hold on
   plot(hbins,y,'r')
   xlabel('data value (strain)')
   ylabel('N')
   titl =  ['N = ',num2str(length(data)),', mean = ',num2str(mu), ...
       ', std = ',num2str(sigma)];
   title(titl)
   print('-dpng',[IFO,'_hhist.png'])

   % plot the first bunch of data points 
   Npoints = 10000;
   datn = data(1:Npoints);
   % make a time vector
   timn = (0:(Npoints-1))/fs;
   figure
   plot(timn,datn,'b')
   xlabel('time (s)')
   ylabel([IFO,'strain data value (strain)'])
   titl =  ['N = ',num2str(length(data)),', mean = ',num2str(mu), ...
       ', std = ',num2str(sigma)];
   title(titl)
   print('-dpng',[IFO,'_hoft.png'])

   % make the asd
   nfft = fs;
   [psdm,freqm] = pwelch(data,hann(nfft),nfft/2,nfft,fs);
   asdm = sqrt(psdm);

   % plot the asd
   figure
   loglog(freqm, asdm,'b')
   xlim([20, 4096]) 
   ylim([1.e-23, 1.e-18]) 
   grid on 
   xlabel('f (Hz)')
   ylabel('ASD')
   title([IFO,' noise ASD around GRB051103'])
   print('-dpng',[IFO,'_asd.png'])

   % make and plot the spectrogram
   figure
   nfft = 1024;
   [y,f,t,p] = spectrogram(data,hann(nfft),nfft/2,nfft,fs);
   surf(t,f,10*log10(abs(p)),'EdgeColor','none')
   axis xy; axis tight; colormap(jet); view(0,90);
   xlabel('Time since 815045078 (s)') 
   ylabel('frequency (Hz)') 
   title([IFO,' spectrogram around GRB051103'])
   %print('-dpng',[IFO,'_spectrogram.png'])
   

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Next, the 2190s of data around GRB051103
   % sampling frequency in Hz and duration in seconds
   fs = 16384;
   Duration = 2190;
   GPSstart = 815043278;

   % open the file and get the data (strain time series):
   IFO = char(IFOs(I));
   fn = [IFO,'-','STRAIN_',num2str(fs),'Hz-',num2str(GPSstart),...
       '-',num2str(Duration),'.hdf5'];
   dataset = [IFO,':LSC-STRAIN'];
   data = hdf5read(fn,dataset);

   % print out some basic data checks
   display(['Contents of ',fn]) 
   display(['IFO, Npoints = ',IFO,', ',num2str(length(data))])
   display(['Duration = ',num2str(length(data)/(fs*Duration)), ...
       ' * ',num2str(fs),' * ',num2str(Duration)])
   display(['min,max,mean,std = ',num2str(min(data)),', ', ...
      num2str(max(data)),', ',num2str(mean(data)),', ',num2str(std(data))])

   
   %% histogram the raw data:
   [hn,hbins] = hist(data, 100);
   % make a corresponding Gaussian:
   N = max(hn);
   mu = mean(data);
   sigma = std(data);
   y = N*exp(-((hbins-mu)/sigma).^2/2.0);
   % plot it:
   figure
   bar(hbins,hn)
   hold on
   plot(hbins,y,'r')
   xlabel('data value (strain)')
   ylabel('N')
   titl =  ['N = ',num2str(length(data)),', mean = ',num2str(mu), ...
       ', std = ',num2str(sigma)];
   title(titl)
   print('-dpng',[IFO,'_hhistf.png'])

   % make the asd
   nfft = fs;
   [psdf,freqf] = pwelch(data,hann(nfft),nfft/2,nfft,fs);
   asdf = sqrt(psdf);

   % now read in the spectrum from the distributed file:
   fn = [IFO,'-','SPEC-',num2str(GPSstart),'-',num2str(Duration),'.hdf5'];
   dataset = [IFO,':SPECTRUM'];
   data = hdf5read(fn,dataset);
   asdd = data(2,:);
   freqd = data(1,:);

   % plot the asd, this time with both the short and long data stretches
   figure
   loglog(freqm, asdm,'b',freqf, asdf,'r',freqd, asdd,'g--')
   xlim([20, 4096]) 
   ylim([1.e-23, 1.e-18]) 
   grid on 
   xlabel('f (Hz)')
   ylabel('ASD')
   legend('256 s @ 4096 Hz','2190 s @ 16384 Hz','SPEC file')
   title([IFO,' noise ASD around GRB051103'])
   print('-dpng',[IFO,'_asd.png'])

end
