% reachability_plot_Norcia_cluster:
%
% 1) computes the reachability plot of the Norcia cluster for Period I
%    and Period II
% 2) selects the events in the first two valleys (C5+C2)and in the first
%    three valleys (C5+C2+C1) of the reachability plot for Period I 
% 3) calls the function pca_valleys and checks if the clusters
%    can be associated to a planar surface;
% 4) saves in the matrix Ap the eigenvalues and the planar surface parameters 
%    (dip angle,strike angle,length,height).

close all


X1 = load('Norcia_I_period_scaling_updated2.txt');
X2 = load('Norcia_II_period_scaling_updated2.txt');


X1c=[X1(:,6) X1(:,7) X1(:,8)];

X2c=[X2(:,6) X2(:,7) X2(:,8)];


maxEpsilon = 1.6;
minNumPoints = 100;

[order1, reachdist1]=clusterDBSCAN.discoverClusters(X1c,maxEpsilon,minNumPoints);
[order2, reachdist2]=clusterDBSCAN.discoverClusters(X2c,maxEpsilon,minNumPoints);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(1,2,1)
bar(reachdist1(order1))
yline(0.85)
hold on 
%yline(0.4)
axis([0 length(X1c(:,1)) 0.4 1.6])
%title('cluster 1 ')
subplot(1,2,2)
bar(reachdist2(order2))
hold on 
%yline(0.4)
axis([0 length(X2c(:,1)) 0.4 1.6])
%title('cluster 2 ')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


X1cm =[X1c X1(:,10)];
red =[X1cm(order1,:) order1' reachdist1(order1)' (1:length(order1))'];
[pks,locs, w, p] = findpeaks(red(:,6));
figure
findpeaks(red(:,6),red(:,7),'MinPeakProminence',0.2,'Annotate','extents')

redvalleys=zeros(length(X1cm(:)),5);
red_valley1 = red(1:locs(3),1:4);       % C5 and C2
red_valley2 = red(locs(1):locs(4),1:4);   % C5, C2 and C1


[dip1,strike1,L1,H1,p1,xp1,yp1,zp1,e11,e12,e13] = pca_valleys(red_valley1);
[dip2,strike2,L2,H2,p2,xp2,yp2,zp2,e21,e22,e23] = pca_valleys(red_valley2);


%% structural parameters
Ap = [e11,e12,e13,dip1,strike1,L1,H1;...
    e21,e22,e23,dip2,strike2,L2,H2];