% Helicorder script for TLAV Vert Comp
% Modified by Branden Christensen from
% June 22, 2009 for NM Tech Field Course

%% Notes
 % For infrasound high pass filter data (above 5s or 0.2 Hz)
%% ISSUES

% Possible to change color of every other line?

%% READ IN DATA
cd '/Users/geop572/Desktop/SEGY/' % move to data directory

start_time = [170 0 0 0]; % [jday hour min sec]
window = 3600*24; % seconds of data
das = '9026'; % name of seismic station datalogger to read in. das = datalogger = 9026
Station_Name = 'TLAV';

[wiggleV,hdr] = segycut(start_time,window,das,'1'); % start_time,window,das,channel

sps = hdr.sps
ts = (1:length(wiggleV))/sps; % create time series vector
seis_sens = 800; % sensitivity of Guralp T40 is 800 V/m/s
atod = hdr.atod; % Anolog to digital. converts counts to volts. 
sens_gain = hdr.gain; % programmed gain setting
seis_calib = atod/seis_sens/sens_gain*1e6; % total conversion of counts to um/s

%% Apply Filter 0.01-10 Hz

npoles = 2
low = 5; % corner frequency at 5 Hz
high = 0.01; % corner frequency at 1 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);

%% now plot the calibrated and filtered data trace for 24 hours
figure(1); clf
decimation = 10; %resamples data as not to include all data points. 
plot(ts(1:decimation:end)/60/60,bandpass_wiggle(1:decimation:end)*seis_calib, 'k') % plot decimated data trace for 24 hours (k = blck)
ylabel('velocity (\mum/s)')
xlabel('time (hours)')
title(['TLAV Vertical Component' '  | start time: ' num2str(start_time)])
XLim([0 24])


%% now plot a traditional helicorder record of the same data
figure(2); clf
decimation = 2;
ts_jday = hdr.day + ts/60/60/24; % time axis for the ogram script must be decimal julian day
scaling = 8; % scaling is set to the window, larger numbers amplify the traces by more
ogram_revised(ts_jday(1:decimation:end),bandpass_wiggle(1:decimation:end),scaling) % this is J. Johnson's script
ylabel('time(hours)')
xlabel('minutes')
title(['TLAV Vertical Component' '  | start time: ' num2str(start_time)])