% histogram and percentile distribution of given channel
% forexample here: seism microseism

format long g
run /users/swinkels/deploy/MatlabVirgoTools/trunk/startup.m %%% this give access BAS's functions, including pole/zero TF fit
clear all
close all

cd '/users/fiori/home/AdV/Commissioning/clima/percentili_rms'

%chname='V1:ENV_CEB_SEIS_N_50Hz_rms_0.1_1Hz';
%chname='V1:ENV_CEB_SEIS_W_50Hz_rms_0.1_1Hz';
chname='V1:ENV_CEB_SEIS_V_50Hz_rms_0.1_1Hz';

fflT='/virgoData/ffl/trend.ffl';

gps=utc2gps(datenum(2015,1,25,0,0,0));

gpsend = utc2gps(datenum(2017,2,26,0,0,0));

len = gpsend - gps;
days = floor(len/86400);

if (chname(17) == 'V')
cleancut = 5E-8;  % V
else
cleancut = 8E-8;  % N,W
end

cebn=frgetvectN(fflT,chname,gps,len);
time = (1:1:length(cebn))'/86400/0.005;

figure(1)
plot(time,cebn); grid on; set(gca,'yscale','log')
cebnclean = cebn;
cebnclean(cebnclean<cleancut) = NaN;  % remove sensor-off data
hold on
plot(time,cebnclean,'r');

nbins = 100;
% X = linspace(min(cebn),max(cebn),nbins);
cebnlog = log10(cebnclean);

X = logspace(min(cebnlog),max(cebnlog),nbins);
N = histc(cebn,X);

ss = sum(N);
Nnorm = N./ss;

days = round(ss/86400/0.005)

figure(3)
subplot(1,3,1)
plot(X,Nnorm,'or-');
set(gca,'yscale','lin')
set(gca,'xscale','log')
ylabel('n.entries / Ntot')
xlabel('seism RMS 0.1-1Hz [m/s]')
grid on
tstring=strrep(chname,'_','\_')
title([tstring,' over ',num2str(days),' days'])
legend(['Ntot = ',num2str(ss)])

MED=mean(cebn(cebn>1E-7));
STD = std(cebn(cebn>1E-7));
ENTRIES = length(cebn(cebn>1E-7));

%annotation('textbox',...
%    [0.33 0.8 0.12 0.1],...
%    'String',{['ENTRIES:  ',num2str(ENTRIES)],['MEAN:  ',num2str(MED), ' m/s'],['STD:  ',num2str(STD), ' m/s']},...
%    'FontSize',12,...
%    'FontName','Arial',...
%    'LineStyle','-',...
%    'EdgeColor',[0 0 0],...
%    'LineWidth',1,...
%    'BackgroundColor',[1 1 1],...
%    'Color',[0 0 0]);

Nnorm = N./ss*100;
CUM = cumsum(Nnorm);

subplot(1,3,2)
plot(X,CUM,'b','LineWidth',2);
set(gca,'xscale','log')
set(gca,'yscale','log')
grid on
xlabel('seism RMS 0.1-1Hz [m/s]')
ylabel('Percentile Level [%]')
title('Cumulative distribution')
axis([1E-8 1E-5 1E-3 100])

subplot(1,3,3)
plot(X,CUM,'b','LineWidth',2);
set(gca,'xscale','log')
set(gca,'yscale','lin')
grid on
xlabel('seism RMS 0.1-1Hz [m/s]')
ylabel('Percentile Level [%]')
title('Cumulative distribution')
axis([5E-7 1E-5 90 100])



% write histograms on file
nomefile = [chname(4:end),'_histo.txt']
fid = fopen(nomefile,'w');

fprintf(fid,'%s\n',['% ',chname,])
fprintf(fid,'%s\n','% brms [m/s], entries [#]') 

for k=1:length(X)
fprintf(fid,'%6.2e  %d\n',X(k), N(k));
end
fclose(fid);


% write PL on file
nomefile = [chname(4:end),'_PL.txt']
fid = fopen(nomefile,'w');

fprintf(fid,'%s\n',['% ',chname,])
fprintf(fid,'%s\n','% brms [m/s], cumulative sum [%]') 

for k=1:length(X)
fprintf(fid,'%6.2e  %6.2e\n',X(k), CUM(k));
end
fclose(fid);



% extract and print some PL values
PL = [1,10,50,75,90,95,99];
pp = [];
for k=1:length(PL)
[tmp,index] = min(abs(CUM-PL(k)));
pptmp = X(index);
pp(k) = pptmp;
stringa = ['percentile ', num2str(PL(k)),': ',num2str(pp(k)),' m/s'];
fprintf(1,'%s\n',stringa)
end


