clear all;
%close all;
clc;
 
load Final_TPicks.mat
 
decmin=X_final-floor(X_final);
second=decmin*60;
start_times2=[Start_time_final X_final'];
starttime=[start_times2(:,1:2) floor(X_final)' second'];
precursor=[0 0 0 -1];
window= 10;
das = 'A988'
%window data
for i=1:length(starttime')
    starttime_new(i,:)=starttime(i,:)+precursor;
    [infV(i,:),hdr]=segycut(starttime_new(i,:),window,das,'1');
end
 
 
 
%calibrate
conv_factor_1INCH_DM_V_V2PA= 1*248.843/(0.01*7/12);
infV_calib=1*conv_factor_1INCH_DM_V_V2PA*hdr.atod/32;
%infV_calib=hdr.atod;
infrasound_c=infV*infV_calib*2*pi*11.2*1000;
ts = (1:length(infV))/hdr.sps;

%flip data
infrasound=-infrasound_c;
 
% figure(1);
% plot(infrasound(3,:))

 ptt=size(infrasound);

 %Remove the DC Offset at the Beginning of the Trace.
 for kjh=1:length(X_final)
    DC_Off(kjh)=mean(infrasound(kjh,1:350));% 
    DC_Off2(kjh,:)=ones(length(infrasound),1)*DC_Off(kjh);
    infrasound2(kjh,:)=infrasound(kjh,:)-DC_Off2(kjh,:);
 end
 %Integrate to find (1)Mass Flux and (2)Volumetric Flux
for nn=1:ptt(1,1)
    signal1_q(nn,:)=cumsum(infrasound2(nn,500:3750),2);
    Mass_Flux_Tot(nn,1)=signal1_q(end)*(1/500);
    Mass_Flux(nn,:)=signal1_q(nn,:)*(1/500);
    Mass_Flux(nn,end)=0;
    I=min(find(Mass_Flux(nn,:)<0));
    Mass_Flux(nn,I:end)=0;
    signal2_q(nn,:)=cumsum(Mass_Flux(nn,:));
    Vol_Flux_Tot(nn,:)=signal2_q(end)*(1/500);
    Vol_Flux(nn,:)=signal2_q(nn,:)*(1/500)/(0.5);%divided by 0.5 to get cumulative volume 0.5 is density of air
end



%--------------------------------------------------------------------------
%Min and Max

%Infrasound Amp

for oo=1:nn
    Infra_min(oo)=min(infrasound2(oo,:));
    Infra_max(oo)=max(infrasound2(oo,:));
end
figure(1)
pk2pk=(abs(Infra_max-Infra_min));
hist(pk2pk,10)
xlabel 'Peak to Peak Amplitude (Pa)';
ylabel 'Number of Events';
title 'Chontilla Infrasound Array: Peak to Peak Amplitude Histogram';


% x = 0:.1:1; y = [x; exp(x)];
%         fid = fopen('exp.txt','wt');
%         fprintf(fid,'%6.2f  %12.8f\n',y);
%         fclose(fid);

%Mass Flux
figure(2)
plot(ts(1:length(Mass_Flux)),Mass_Flux)%Mass flux was collected sequentially, thus plotting it will show a sequential progression.
xlabel('Time (s)')
ylabel('Mass Flux (Pa*s)')
axis tight
title 'Chontilla Infrasound Array: Mass Flux vs. Time';

figure(3)
hist(Mass_Flux_Tot,10)
xlabel('Total Mass Flux (Pa*s)')
ylabel('Number of Events')
axis tight
title 'Chontilla Infrasound Array: Total Mass Flux Histogram';

%Volumetric Flux, Q

figure(4)
plot(ts(1:length(Vol_Flux)),Vol_Flux)%Like mass flux, vol flux shows a sequential progression.
xlabel('Time (s)');
ylabel('Cumulative Volume Flux (Pa*s*s)')
axis tight
title 'Chontilla Infrasound Array: Cumulative Volume Flux vs. Time';

figure(5)
hist(Vol_Flux_Tot)
xlabel('Total Cumulative Volume (Pa*s*s)')
ylabel('Number of Events')
axis tight
title 'Chontilla Infrasound Array: Total Cumulative Volume Flux Histogram';

% for kk=1:length(X_final)
%     plot(infrasound(kk,:))
%     title(['das: ' das ' start time: ' num2str(starttime_new(kk,:))])
%     pause
% end

figure(6)
hist(Start_time_final(:,1),14)
xlabel('2009 Julian Day')
ylabel('Number of Events')
axis tight
title 'Chontilla Infrasound Array: Number of Events per Day';