%% Load data
close all;
Vpp=load('/Users/geop572/Desktop/GEOP572_MATLAB/data.txt');

VppInSeEv=Vpp(:,1);
VppInA=Vpp(:,3) * (248.843/(0.01*7/12));
VppSeA=Vpp(:,6) ./800 * 100000;
VppInD=Vpp(:,2);
VppSeD=Vpp(:,5);
VppInF=Vpp(:,4);
VppSeF=Vpp(:,7);

Jday=[164 165 166 167 168 170];

%% Peak to Peak ratio
figure(1); clf;
scatter(VppSeA(1:11,:),VppInA(1:11,:),'r*');
hold on
scatter(VppSeA(12:26,:),VppInA(12:26,:),'b*');
scatter(VppSeA(27:30,:),VppInA(27:30,:),'c*');
scatter(VppSeA(31:38,:),VppInA(31:38,:),'m*');
scatter(VppSeA(39:47,:),VppInA(39:47,:),'k*');
scatter(VppSeA(48:51,:),VppInA(48:51,:),'g*');
title('Peak-to-Peak Amplitude Ratios');
xlabel('Seismic (\mum/s)');
ylabel('Infrasound (Pa)');
legend('164','165','166','167','168','170');

%% Statistical histograms

%Seismic mean and standard deviation calculations
M=mean(VppSeA);
SeMean=strcat('MEAN = ',num2str(M));


StD=std(VppSeA);
SeStanDev=strcat('Stan Dev = ',num2str(StD));


%Infrasound mean and standard deviation calculations
M=mean(VppInA);
InMean=strcat('MEAN = ',num2str(M));


StD=std(VppInA);
InStanDev=strcat('Stan Dev = ',num2str(StD));



figure(2); clf;

subplot(2,1,1)
hist(VppSeA,10)
h = findobj(gca,'Type','patch');
set(h,'FaceColor','c','EdgeColor','k')
hold on
title('Seismic Peak-to-Peak Velocity Distribution');
xlabel('Velocity (\mum/s)');
ylabel('# of Events');
annotation('textbox',[0.7 0.8 .1 .1], 'string',SeMean, 'LineStyle', 'none')
annotation('textbox',[0.7 0.75 .1 .1],'string',SeStanDev,  'LineStyle', 'none')

subplot(2,1,2)
hist(VppInA)
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','k')
title('Infrasonic Peak-to-Peak Pressure Distribution');
xlabel('Pressure (Pa)');
ylabel('# of Events');
annotation('textbox',[0.7 0.3 .1 .1], 'string',InMean, 'LineStyle', 'none')
annotation('textbox',[0.7 0.25 .1 .1],'string',InStanDev,  'LineStyle', 'none')
% patch or fill

%% frequency vs duration
figure(3);clf;
subplot(2,1,1)
plot(VppInSeEv,VppInF,'*')
title('Infrasonic Frequency Variation by Day');
xlabel('Julian Day');
ylabel('Frequency (Hz)');


subplot(2,1,2)
plot(VppInSeEv,VppSeF,'*')
title('Seismic Frequency Variation by Day');
xlabel('Julian Day');
ylabel('Frequency (Hz)');


%% Events and Cummulative Amplitude per day
%events
counts=[11 14 4 8 9 4];

%cummulative amplitude
CummulIn=[sum(Vpp(1:11,3)), sum(Vpp(12:26,3)), sum(Vpp(27:31,3)), sum(Vpp(32:38,3)), sum(Vpp(39:47,3)), sum(Vpp(48:51,3))];
CummulSe=[sum(Vpp(1:11,6)), sum(Vpp(12:26,6)), sum(Vpp(27:31,6)), sum(Vpp(32:38,6)), sum(Vpp(39:47,6)), sum(Vpp(48:51,6))];


figure(4);
subplot(2,2,1:2)
bar(Jday,counts, 'group')
h = findobj(gca,'Type','patch');
set(h,'FaceColor','k','EdgeColor','k')

title('# Events per Day');
xlabel('Julian Day');
ylabel('# of Events');
XLim([163 171])


subplot(2,2,3)
bar(Jday,CummulIn,'group')
title 'Group'
XLim([163 171])
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','k')
title('Infrasound Cummulative Amplitude per Day');
xlabel('Julian Day');
ylabel('Total Amplitude (volts)');
XLim([163 171])

subplot(2,2,4)
bar(Jday,CummulSe, 'group')
h = findobj(gca,'Type','patch');
set(h,'FaceColor','c','EdgeColor','k')

title('Seismic Cummulative Amplitude per Day');
xlabel('Julian Day');
ylabel('Total Amplitude (volts)');
XLim([163 171])

%%
% % % Seismic Traces
% % high Seismic/ low infrasound
% % addpath '/Users/geop572/Desktop/TDATA2/SEGY/' % add to data directory
% % 
% % 
% % start_time = [165 14 29 0]; % [jday hour min sec]
% % :)
% % window = 90; % seconds of data
% % das = 'A871'; % name of seismic station datalogger to read in. das = datalogger = 9026
% % 
% % [wiggleV1,hdr1] = segycut(start_time,window,das,'1'); % start_time,window,das,channel
% % [wiggleV4,hdr4] = segycut(start_time,window,das,'4'); % start_time,window,das,channel
% % 
% % ts1 = (1:length(wiggleV1))/hdr.sps; % create time series vector
% % seis_sens = 800; % sensitivity is 800 V/m/s
% % atod1 = hdr.atod; % Anolog to digital. converts counts to volts. 
% % sens_gain1 = hdr.gain; % programmed gain setting
% % seis_calib1 = atod/seis_sens/sens_gain*1e6; % total conversion of counts to um/s
% % 
% % ts4 = (1:length(wiggleV))/hdr.sps; % create time series vector
% % infra_sens = 800; % sensitivity is 800 V/m/s
% % atod4 = hdr.atod; % Anolog to digital. converts counts to volts. 
% % sens_gain4 = hdr.gain; % programmed gain setting
% % infra_calib = atod/infra_sens/sens_gain4*1e6; % total conversion of counts to um/s
% % 
% % figure(5); clf;
% % plot(wiggleV1)