% Load the lithological column
!awk '{print $1, $2, $4}' ../Data_clean/U1309D_lithunit.dat > lithomat.dat
load lithomat.dat
l=lithomat;clear lithomat;

figure(fignum);
clf

% load the cataclastic data
a=dci;

display('start arrays')
load U1309D_arrays1cm
L=Hlith;A=Hcata;
clear U1309D_arrays1cm Hvein Hopen Hdril Hserp Hcata Hmagm Hplas Hlith

% Define intervals of lithologies and of values for averages
lg=[1  14;    19 25;    25  31;    31  45;    45  65];
in=[0 .5 1 1.5 2 3 5];

N=zeros(length(in)-1,length(lg(:,1)));

display ('loop over lithologies')
for j=1:length(lg(:,1));
    for i=1:length(in)-1;
        ig=find(L(:,2)>lg(j,1) & L(:,2)<lg(j,2) & A(:,2)>=in(i) & A(:,2)<in(i+1));
        N(i,j)=length(ig);
    end
end

dz=.01;
% diagram with actual cored distances
N=N*dz;nN=N;
clf
subplot(2,1,1)
bar(N','stacked'), colormap(cool)
for i=1:length(N(1,:));
    nN(:,i)=100*nN(:,i)./sum(nN(:,i));
axis([0 6 0 620])
ylabel('m')
title('Cataclasis vs lithology, col [0 .5 1 1.5 2 3 5]')
text(3,500,'1 - Diabase')
text(3,465,'2 - Gabbro')
text(3,430,'3 - Oxide gabbro')
text(3,395,'4 - Olivine & troctolitic gabbro')
text(3,360,'5 - Troctolite & dunitic troctolite')
end

% diagram with % for each lithology group
subplot(2,1,2)
bar(nN','stacked'), colormap(hot)
axis([0 6 00 100])
ylabel('%')

% 1 basalt & diabase
% 2 gabbro
% 3 Ox gabbro
%4 Ol & troct gabbro
% 5 Troctolite
% 6 wehr - dun - haz
