% RSAM Script
% Station TLAV channel 1
% Written by Brando and Carolyn
% June 22, 2009 NM Tech Field Course
%%
%  RSAM Notes - Real-time Seismic Amplitude
% 	synthesizing amount of energy released over time
% 	good diagnostic of variation in seismic energy over time
% 	qualitative tool	
%%
% The Goal of this script is to:
% 	Step 1 : Remove the mean of trace or apply high pass filter to zero center trace
% 	Step 2 : Compute abs value of waveform
% 	Step 3 : Compute average amplitude of waveform for overlapping 10 min intervals ("running average") (Why 10 min intervals?  historical)
% 	Step 4 : Sum individual time windows

%% To improve this script:

% use a function other than reshape to obtain overlapping 10 min intervals
% ("running average")

%expand to all days instead of just 7
%further expand to all DAS units on all days

%%
clear all
k = 1; % Applies to creation of RSAM matrix (more infor? See second to last cell)

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


for i=166:172 % This is a for loop for Julian days 166:172

start_time = [i 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

ts = (1:length(wiggleV))/hdr.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

%% RESHAPE THE MATRIX, REMOVE THE MEAN FROM THE SIGNAL, TAKE ABS VALUE
% "RESHAPE" = divides data into 10 minute chunks. Why 10 min? Standard
% methodology

% Is "Detrend" command better here?

wiggleV2 = wiggleV(1:round(length(wiggleV)/60000)*60000,1); % This was done b/c on Julian day 172 there were only 92 ten minute intervals
wiggle_10 = reshape(wiggleV2,60000,length(wiggleV2)/60000); % 60000 rows of samples (60sec*10min*100samples per second), 144 columns (24hrs*6 10 min intervals per hour)..  i.e. 1 10 minute interval = 1 column
wiggle_mean = mean(wiggle_10,1); % takes the mean of each column, returning 1 row, all columns as means

wiggle_removed = zeros(size(wiggle_10)); % this creates a matrix to be filled in the next for loop

for j = 1:length(wiggle_mean)
wiggle_removed(:,j) = wiggle_10(:,j) - wiggle_mean (1, j); % a : means all rows or all columns pending which part of the matrix you call
% The mean of the wiggle is taken to center the wiggles at zero (so that
% baseline does not have a positive amplitude)
% This can also be accomplished by taking a high pass filter
end

wiggle_abs = abs(wiggle_removed); % This takes the abs value of the wiggle so all amplitudes are positive

%% Average amplitude of waveform for overlapping 10 min intervals

if length(wiggleV2)/60000 == 144  % NOTE: 24 hrs = length of 144
    RSAM_all(k,:) = sum(wiggle_abs,1)/1.0e+07; 
else
    RSAM_all(k,:) = NaN; % NaNs dont plot, just fill
    % This creates a matrix filled with NaNs, many or all of which are
    % replaced in the next line:
    RSAM_all(k,1:length(wiggleV2)/60000) = sum(wiggle_abs,1)/1.0e+07;

end

k = k+1;
end
 
%% FIGURES!

time24 = 0:1/6:23.84; % Done for x-axis
figure(1); clf;

subplot(4,2,1) % subplot(#rows, # columns, position)
plot(time24, RSAM_all(1,:))
xlabel('time (JDAY 166)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM'); %TLAV = name of station associated with DAS 9026
axis([0 24 0 10])

subplot(4,2,2)
plot(time24, RSAM_all(2,:))
xlabel('time (JDAY 167)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])

subplot(4,2,3)
plot(time24, RSAM_all(3,:))
xlabel('time (JDAY 168)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])

subplot(4,2,4)
plot(time24, RSAM_all(4,:))
xlabel('time (JDAY 169)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])

subplot(4,2,5)
plot(time24, RSAM_all(5,:))
xlabel('time (JDAY 170)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])

subplot(4,2,6)
plot(time24, RSAM_all(6,:))
xlabel(' time (JDAY 171)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])

subplot(4,2,7)
plot(time24, RSAM_all(7,:))
xlabel('time (JDAY 172)');
ylabel('RSAM x1.0E07');
title('TLAV Vertical RSAM');
axis([0 24 0 10])