%cd '/Users/Jeff/Documents/DATA/KILAUEA2008/SEGY' % move to data directory
cd '/Users/andrea/Desktop/SEGY_tungu09'

%for day=0:1
    start_time = [169 13 52 12]; % rock fall
    %start_time = [172 8 28 31]; % LP
    %start_time = [170 18 12 57]; % lahaar
    %start_time = [168 8 49 5]; % explosion
    window = 3600*1/60; % 24 hours or 86400 s of data
    das = 'A860'; % name of seismic station to read in

    [seisV,hdr] = segycut(start_time,window,das,'1');
    ts = (1:length(seisV))/hdr.sps; % create time series vector
    
    % Convert to Volts to micro meter per sec
    seis_sens = 800; % sensitivity is 800 V/m/s
    atod = hdr.atod; % digitization V/count
    sens_gain = hdr.gain; % programmed gain setting
    seis_calib = atod/seis_sens/sens_gain*1e6; % total conversion of counts to mum/s
    %convert = seisV*seis_calib;
    
    %Detrend  removes the means
    
    detr_seis = detrend(seisV*seis_calib);
    
    % filtered seisV
   % sps = 100; % get sample rate from header
    npoles = 2; % specify order of filter
    low = 0.01; % corner frequency at 0.01 Hz
    high = 49;  % corner frequency at 25 Hz
    [B1,A1] = butter(npoles,[low high]/(hdr.sps/2));
    bandpassV = filter(B1,A1,detr_seis);
    
    time=linspace(0,3600*1/60,length(detr_seis));
    
    % Spectrogram
    nfft = 512;
    subplot(2,1,1);spectrogram(detr_seis,[] , [] ,nfft,hdr.sps,'yaxis');
    
    
    subplot(2,1,2),plot(time,detr_seis);
    xlim([0 7200])
    xlabel('time (secs)')
    ylabel('velocity (\mum/s)')