s%% Script written by Omar and Hunter
% Revised by Brando and Ezer
% June 24 2009 NM Tech Field Course

% This code will load an hour of data from station X, allow the user to
% pick the onset of an event and write that pick to a ascii text file. 

% Threshold 10 Pa

%% READ IN DATA
clear all 
close all

cd '/Users/geop572/Desktop/SEGY/'
j = input('What Julian day?');
h = 1;
l = 1;

for i = 1:24
% NOTE: data for 161 starts at 21:00
% NOTE: for 164 16th hour does not exist
% NOTE: resume picking on 168 at 19th hr
start_time(l,:) = [j i-1 0 0]; % [jday hour min sec] of start
window = 3600*1; % 1 hour of data
das = 'A988'; % name of seismic station to read in
Station_Name = 'TINF';
% A988 = TINF (Infrasound only array)

[wiggleV,hdr] = segycut(start_time(l,:),window,das,'6');

% Using channel 6 of TINF b/c best signal at infrasound only array

sps = hdr.sps;
ts = (1:length(wiggleV))/sps; % create time series vector
infra_sens = 1*248.843/(0.01*7/12); 
% Sensitivity of Jeff's sensor = 23 microVolts/ Pa | 7 Voltage control
% oscillator
atodtoPa = wiggleV*infra_sens*hdr.atod/32; %atodtoPa analog_to_dig_to_Pascals
infra_gain = hdr.gain; % programmed gain setting
wiggleV_detrend=detrend(atodtoPa);

%% Butter filter data from 0.2 to 10 Hz

npoles = 2;
low = 10; % corner frequency at 10 Hz
high = 0.2; % corner frequency at 0.2 Hz
% signal processing tool box with butter is necessary for following operations
% first do bandpass

[B1,A1] = butter(npoles,[high low]/(sps/2)); % The order of high and low MATTERS. lowest number must come first.
% now do highpass

bandpass_wiggle = filter(B1,A1,wiggleV_detrend);


%% Plot the calibrated data trace for 1 Hour
figure(1); clf
subplot(2,1,1);
plot(ts,bandpass_wiggle,'b') % plot decimated data trace for 1 hour
ylabel('Pa')
xlabel('time (s)')
title([Station_Name 'Ch.6' ' start time: ' num2str(start_time(l,:))])

% Spectrogram for the calibrated data trace for 1 Hour
subplot(2,1,2);
nfft = 512;
spectrogram(bandpass_wiggle,10000 , [] ,nfft,sps,'yaxis');
xlabel('time(s)')
ylabel('Hz')
YLim([0 50])

%% Pick the onset of the event(s)
num_events = input('How many events are evident in this Hour Trace?');
% Pick the event using ginput. Don't forget to use the zoom.
    for k = 1:num_events
        [x(h),y(h)] = ginput(1)
    if num_events == 0
        x(h) = NaN;
    end
    pause
    h = h+1;
    l = l+1;
    if num_events > 1
        start_time(l,:) = start_time(l-1,:);
    end
    end
    
end 
start_timef=start_time(1:length(x),:);

% Save your picks in a file called Tpicks.mat
save Tpicks.mat start_timef x

% Method:
% 1. Look at number of events in window 
% 2. Zoom into event
% 3. Enter number of events from 1. in Command Window
% 4. Hit space bar to select next event.
% Once all events have been entered hit Enter