clear

q=input('Avez vous indiqué au programme le numéro d''ordre des éléments ?','s');

nO=1;
%nso=8;
nmg=2;
nal=3;
nsi=4;
nph=5;
%le soufre s'appelait vif-argent
nva=6;
%npo=1;
nca=7;
nti=8;
ncr=9;
nmn=10;
nfe=11;
nni=12;

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

% Chargement des données

u=['load c:\matthieu\matlab\d89g03.dat']
eval(u);

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

% Si il s'agit de microsonde, il faut mutiplier les données par 100 pour
% les lire en %

wo=d89g03;
% Parfois les donnees ne sont pas rangées sous le bon format
%wo=wo';

[np,ne]=size(wo)

% np : nombre de points d'analyse
% ne : nombre d'éléments chimiques retenus

% A partir de cette ligne, les donnees doivent etre presentees convenablement.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Définition des masses moléculaires.

%mmol(1,npo)=94.20;	%K2O
mmol(1,nca)=56.08;	%CAO
mmol(1,nti)=79.88;	%TIO2
mmol(1,ncr)=151.99;	%CR2O3
mmol(1,nmn)=70.94;	%MNO
mmol(1,nfe)=71.85;	%FEO
mmol(1,nni)=74.69;	%NIO
%mmol(1,nso)=61.98;	%NAO
mmol(1,nsi)=60.08;	%SIO2
mmol(1,nal)=101.96;	%AL2O3
mmol(1,nmg)=40.30;	%MGO
mmol(1,nva)=32.06;	%S
mmol(1,nph)=141.94;	%P2O5

mmol(1,nO)=1 %O??

mmol=mmol';


for i=1:np
for j=1:ne

rm(i,j)=wo(i,j)./mmol(j);

end
end

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

% Pour calculer le nombre d'O il faut connaitre le nombre d'atomes d'O pour 
% chaque molécule.

%mine(1,npo)=1;
mine(1,nca)=1;
mine(1,nti)=2;
mine(1,ncr)=3;
mine(1,nmn)=1;
mine(1,nfe)=1;
mine(1,nni)=1;
%mine(1,nso)=1;
mine(1,nsi)=2;
mine(1,nal)=3;
mine(1,nmg)=1;
mine(1,nva)=0;
mine(1,nph)=5;

mine=mine';

for i=1:np
for j=1:ne

rmO(i,j)=rm(i,j).*mine(j);

end
end

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

% Calcul du total de l'analyse et du nombre total d'O déterminé par 
% stoichiométrie.

somE=sum(wo');

somO=sum(rmO');

%size(somO)


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

% Calcul de la stoichiométrie.

for i=1:np
for j=1:ne
st(i,j)=rm(i,j)*24/somO(i);
end
end

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

% somcat est la somme de la stoichiométrie de tous les cations.
% rcasi est le rapport de la somme de tous les cations sur le silicium.

somcat=sum(st');

for i=1:np
rcasi(i)=(somcat(i)-st(i,nsi))/st(i,nsi);
end

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

% On fait un tri des points pour déterminer qui est olivine et qui est 
% pyroxène.
% A la fin de cette procédure, la matrice ol contient 1 quand le point 
% analysé est une olivine. Idem pour px. Pour vérifier que ce n'est pas 
% stupide on peut faire ol+px (on ne peut être à la fois olivine et pyroxène)

ol=ones(1,np);
px=ones(1,np);

for i=1:np



if somE(i)<97.5
px(i)=px(i)*0;
ol(i)=ol(i)*0;
end

if somE(i)>102.5
px(i)=px(i)*0;
ol(i)=ol(i)*0;
end

if somcat(i)>16.2
px(i)=px(i)*0;
end

if somcat(i)<15.8
px(i)=px(i)*0;
end

if somcat(i)>18.2
ol(i)=ol(i)*0;
end

if somcat(i)<17.8
ol(i)=ol(i)*0;
end

if rcasi(i)<0.9
px(i)=px(i)*0;
end

if rcasi(i)>1.1
px(i)=px(i)*0;
end

if rcasi(i)<1.8
ol(i)=ol(i)*0;
end

if rcasi(i)>2.2
ol(i)=ol(i)*0;
end

end

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

% On redéfinit les matrices de données en se débarassant des points qui ne
% sont ni pyroxènes ni olivines.

ko=1;
kp=1;

for i=1:np

if ol(i)==1
wool(ko,:)=wo(i,:);
rmol(ko,:)=rm(i,:);
fa(ko)=100*rmol(ko,nfe)/(rmol(ko,nfe)+rmol(ko,nmg));
ko=ko+1;
end

if px(i)==1
wopx(kp,:)=wo(i,:);
rmpx(kp,:)=rm(i,:);
fs(kp)=100*rmpx(kp,nfe)/(rmpx(kp,nfe)+rmpx(kp,nmg)+rmpx(kp,nca));
en(kp)=100*rmpx(kp,nmg)/(rmpx(kp,nfe)+rmpx(kp,nmg)+rmpx(kp,nca));
kp=kp+1;
end

end

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

% Exemple de tri : on fabrique des matrices qui ne contiennent que des grains
% de taille 25-50 et que des grains de taille 50-100
% taillepx et taille ol permettent de trier 1 code les petits grains
%					    0 code les gros grains
%					    -1 code un grain aberrant
		

%taillepx=[-1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1]

%tailleol=[1 1 1 0 0 0 0 0 0 0 1]


%npx=sum(px)
%nol=sum(ol)

% npx donne le nombre total de pyroxènes
% nol donne le nombre total d'olivines

%clear ko
%clear kp

%ko=1;
%kp=1;

%for i=1:nol

%wool1(ko,:)=wool(i,:);
%rmol1(ko,:)=rmol(i,:);
%fa1(ko)=100*rmol(ko,nfe)/(rmol(ko,nfe)+rmol(ko,nmg));
%ko=ko+1;
%end

%if px(i)==1
%wopx(kp,:)=wtox(i,:);
%rmpx(kp,:)=relmol(i,:);
%fs(kp)=100*rmpx(kp,nfe)/(rmpx(kp,nfe)+rmpx(kp,nmg)+rmpx(kp,nca));
%kp=kp+1;
%end







