clear

%%%%%%%%%%%%% paramètres physiques %%%%%%%%%%%%%%%%%%%
St=5.1e14;
S1=5.76e14;
Sm=5.75e14; 
S2=5.78e14;     % surface de la sphère définie par l'orbite du satellite en m2
A=3.16e7;       % un an exprimé en secondes  
T=5.77;	       % durée d'exposition du satellite	en années.
teta=S2*A/1e9;  %permet de convertir pour Brownlee un flux en g/cm2 en un flux sur toute la terre, exprimé en milliers de tonnes 
rho=2.5;        % densité des particules en g/cm3

ST=17;          %Surface de collection de Taylor en m2;
AT=430;         %Durée de la collection de Taylor en années.
tetat=St*A/1e9; %facteur de conversion pour Taylor.


trou=4;         %Correction apportée a posteriori à la durée de vie des trous à cryoconite.
tetar=St*A/1e9*trou; %facteur de conversion pour Robin.



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

%Brownlee.

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

npasb=1000;

xb1=-8.22; % x1 et x2 sont les bornes entre lesquelles varie le log de la masse du micrométéoride
xb2=-4.12;
dxb=(xb2-xb1)/(npasb-1);
%x est le vecteur log de la masse du micrométéoride (en g)

xb=xb1:dxb:xb2;
mb=10.^(xb);

%npas est le nombre de points contenus dans le vecteur x.

npasb=1+(xb2-xb1)/dxb;

%y est le logarithme du diamètre de la particule (en cm)
%yy est le logarithme du diamètre de la particule (en µm)
%yyy est le diamètre de la particule (en µm)

y=xb/3-0.04;
yy=y+4;
yyy=10.^yy;


%z est le logarithme du diamètre du cratère (en cm).
 %zz est le logarithme du diamètre du cratère (en µm).

z=1.056*y+0.69;
zz=4+z;

pb(6)=8.8749;
pb(5)=-25.439;
pb(4)=28.054;
pb(3)=-15.536;
pb(2)=4.1732;
pb(1)=-0.44175;

alpha=polyval(pb,zz);
%alpha=ab+bb*zz+cb*zz.^2+db*zz.^3+eb*zz.^4+fb*zz.^5;

for j=1:npasb
   ncb(j)=(10.^alpha(j))/T/A*10000; %ncb flux cumulé Brownlee en nombre de météorides par m2 et par seconde.
end



%FLUX DIFFERENTIEL DE Brownlee


nb(1,npasb)=0;


for j=1:npasb-1   
   nb(j)=ncb(j)-ncb(j+1);
   phib(j)=teta*mb(j+1)*nb(j)/dxb;
   
end

phib(1,npasb)=phib(1,npasb-1);

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

%TAYLOR.

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





u=['load c:\matthieu\matlab\datataylor.dat;';] % données Taylor.
eval(u);

xt=datataylor(:,2); %les données sont rangées à l'envers.
yt=datataylor(:,1);  %Le flux cumulatif de Taylor est donné en nombre de sphérules.
yt=yt-log10(ST*AT*A); %Il faut le transformer en nombre de sphérules par m2 et par seconde.
%yt0=yt0-10.51+0.48;

%flux différentiel  à partir du polynome.

dxt=dxb;

xt1=-6.79;
xt2=-3.35;
xtt=xt1:dxt:xt2;
mt=10.^xtt;

npast=1+round((xt2-xt1)/dxt);
pt=polyfit(xt,yt,5);
lnct=polyval(pt,xtt);

for j=1:npast
   nct(j)=(10.^lnct(j));
   end

for j=1:npast-1   
   nt(j)=nct(j)-nct(j+1);
   phit(j)=tetat*mt(j+1)*nt(j)/dxt;
end

phit(1,npast)=phit(1,npast-1);



%flux différentiel brut.

mt0=10.^xt;
npast0=size(xt);

for j=1:npast0
   nct0(j)=(10.^yt(j));
end

nt0(1,npast)=0;

delta=4;

for j=1:(npast0-1)/delta
    dxt0(j)=-xt(1+(j-1)*delta)+xt(j*delta+1);
    xtt0(j)=(xt(1+(j-1)*delta)+xt(j*delta+1))/2;
    mmt0(j)=mt0(1+j*delta-(delta)/2);
    nt0(j)=nct0(1+(j-1)*delta)-nct0(j*delta+1);
    phit0(j)=tetat*mmt0(j)*nt0(j)/dxt0(j);
end




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

%Courbe de robin%

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


%Les données recopiées dans la thèse de Robin sont des particules par m2 et par an.

u=['load c:\matthieu\matlab\datarobin.dat;';] % données Robin pour toutes les micrométéorites.
eval(u);
xr=datarobin(:,1);
yr=datarobin(:,2);
yr=yr-log10(A); %Il faut transformer en nombre de particules par m2 et par seconde.
%yr0=yr0-6.76+0.02-0.145;


u=['load c:\matthieu\matlab\datarobspS.dat;';] % données sphérules S Robin.
eval(u);
xrS=datarobspS(:,1);
yrS=datarobspS(:,2);
yrS=yrS-log10(A); %Il faut transformer en nombre de particules par m2 et par seconde.

%yrS=yrS-6.45;

u=['load c:\matthieu\matlab\datarobspV.dat;';] % données sphérules V Robin.
eval(u);
xrV=datarobspV(:,1);
yrV=datarobspV(:,2);
yrV=yrV-log10(A); %Il faut transformer en nombre de particules par m2 et par seconde.

%yrV=yrV-5.60;

yrA=log10((10.^yrS+10.^yrV)); %somme des 2 types de sphérules.

%yrA=yrA-0.35;

%xrS et xrA et xr sont presque identiques : la différence est due à l'erreur du scanner.

%FLUX DIFFERENTIEL ROBIN%


% Fit polynomial du flux cumulatif de micrométéorites de Robin.

dxr=dxb;

xr1=-7;
xr2=-3.4;

xrr=xr1:dxr:xr2;

pr=polyfit(xr,yrA,5);
lncr=polyval(pr,xrr);


mr=10.^xrr;
npasr=1+round((xr2-xr1)/dxr);

for j=1:npasr
   ncr(j)=(10.^lncr(j));
end

nr(1,npasr)=0;


for j=1:npasr-1   
   nr(j)=ncr(j)-ncr(j+1);
   phir(j)=tetar*mr(j+1)*nr(j)/dxr;
end

phir(1,npasr)=phir(1,npasr-1);

%Flux différentiel brut.

mr0=10.^xr;
npasr0=size(xr);

for j=1:npasr0
   ncr0(j)=(10.^yrA(j));
end

nr0(1,npasr0)=0;
delta=2;


 for j=1:(npasr0-2)/delta
     dxr0(j)=-xr(1+(j-1)*delta)+xr(j*delta+1);
     xrr0(j)=(xr(1+(j-1)*delta)+xr(j*delta+1))/2;
     mmr0(j)=mr0(1+j*delta-(delta)/2);
     nr0(j)=ncr0(1+(j-1)*delta)-ncr0(j*delta+1);
     phir0(j)=tetar*mmr0(j)*nr0(j)/dxr0(j);
end


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

%Antarctique 1998.

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

%flux cumulatif des sphérules 1998

u=['load c:\matthieu\matlab\spcom.dat;';] %mes données de la collecte 1998.
eval(u);

dp98=spcom; %liste des diamètres des sphérules en µm  de la collecte 1998.

ddp98=dp98*1e-4; %liste des diamètres des sphérules en cm de la collecte 1998.

mi=3.14/6*rho.*ddp98.*ddp98.*ddp98; %liste de la masse des sphérules en grammes .
                                    %mi masse intermédiaire.

ms98=sort(mi); %liste ordonnée de la masse des sphérules en grammes 

[npas980,ni]=size(ms98); %np98 nombre de sphérules ni n'a aucun sens.

ms98=ms98';

x98=(log10(ms98));


nc980=npas980:-1:1; %ncum98 nombre cumulatif de sphérules.
nc980=2.32e-8*nc980; %comme nc980 est déterminé à une constante près...


%Fit polynomial avec un potentiel de Woods et Saxon.

dx98=dxb;
x981=-7.42;
x982=-4.87;
xx98=x981:dx98:x982;

%Potentiel de Woods et Saxon.

f1=1.23e-6;
a1=0.15;
x1=-6.1;
xw1=-7.5:0.01:-4.5;
ff=-f1./(1+exp((-xx98+x1)./a1))+1.25e-6;

npas98=1+round((x982-x981)/dx98);

m98=10.^xx98;

for j=1:npas98-1
   n98(j)=ff(j)-ff(j+1);
   phi98(j)=teta*m98(j)*n98(j)/dx98; %On multiplie par teta mais ce n'a pas de sens puisqu'on ne connait pas le temps d'exposition.
end
phi98(1,npas98)=phi98(1,npas98-1);

%Différenciation à partir des données brutes.

n980(1,npas980)=0;

delta=4;

for j=1:(npas980-2)/delta
   dx980(j)=-x98(1+(j-1)*delta)+x98(j*delta+1);
   j
    xx980(j)=(x98(1+(j-1)*delta)+x98(j*delta+1))/2;
    mm980(j)=ms98(1+j*delta-(delta)/2);
    n980(j)=nc980(1+(j-1)*delta)-nc980(j*delta+1);
    phi980(j)=teta*mm980(j)*n980(j)/dx980(j);
end




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

%points de l'Antarctique 1994 routine%%%

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

xa=[-7.69 -6.78 -5.88 -4.07];
ya=[-6.44 -6.50 -6.68 -9];
ya=ya+0.385;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FLUX DIFFERENTIEL%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit polynomial du flux cumulatif de l'Antarctique 1994.


ma=10.^xa;
pa=polyfit(xa,ya,3);


xa1=-7.69;
xa2=-4.07;

dxa=dxb;

xaa=xa1:dxa:xa2;
mba=10.^xaa;

npasa=1+round((xa2-xa1)/dxa);

lnca=polyval(pa,xaa);

for j=1:npasa
   nca(j)=(10.^lnca(j));
end

na(1,npasa)=0;




for j=1:npasa-1   
   na(j)=nca(j)-nca(j+1);
   phia(j)=teta*mba(j+1)*na(j)/dxa;
end

phia(1,npasa)=phia(1,npasa-1);







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



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

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

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

%FLUX INTEGRAL
%Renormalisation de Antarctique 1998.
%[k0,N98]=size(phiw98)
%int98(1,1)=phiw98(1,1)*dxw98;

for j=1:N98-1

%int98(1,j+1)=int98(1,j)+phiw98(1,j+1)*dxw98;
end

N=int98(1,N98)
Phi98=phiw98;


%Calcul du flux intégral de Brownlee.

intb(1,1)=phib(1,1)*dxb;
 
for j=1:npasb-1
   intb(1,j+1)=intb(1,j)+phib(1,j+1)*dxb;
end
Phib=phib/intb(1,npasb)*N;
%B=intb(1,npasb)
%Renormalisation de Taylor.

intt(1,1)=phit(1,1)*dxt;

for j=1:npast-1
   intt(1,j+1)=intt(1,j)+phit(1,j+1)*dxt;
end

Phit=phit/intt(1,npast)*N;


%Renormalisation de Robin.

intr(1,1)=phir(1,1)*dxr;

for j=1:npasr-1
   intr(1,j+1)=intr(1,j)+phir(1,j+1)*dxr;
end

Phir=phir/intr(1,npasr)*N;


%Renormalisation de Antarctique 1994.

inta(1,1)=phia(1,1)*dxa;

for j=1:npasa-1
   inta(1,j+1)=inta(1,j)+phia(1,j+1)*dxa;
end

Phia=phia/inta(1,npasa)*N;

%renormalisation avec Brownlee comme enveloppe.

PhiB=Phib;
PhiA=Phia/13.6*9.7;
PhiT=Phit/13.25*9.5;
PhiR=Phir/7*9.9;
PhI98=(Phi98)/33.4*5.7;

%calcul du flux intégral relatif au flux intégral total :rel
%calcul du flux renormalisé à celui de Brownlee.

for j=1:npasb
   relb(1,j)=intb(1,j)/intb(1,npasb);
   brob(1,j)=40*phib(1,j)/intb(1,npasb);
 end
 
 
 for j=1:npasr
   relr(1,j)=intr(1,j)/intr(1,npasr);
    end

for j=1:npast
   relt(1,j)=intt(1,j)/intt(1,npast);
end

for j=1:npasa
   rela(1,j)=inta(1,j)/inta(1,npasa);
    end
    
    for j=1:N98
       rel98(1,j)=int98(1,j)/int98(1,N98);
    end

