clear

%%%%%%%%%%%%% paramètres physiques %%%%%%%%%%%%%%%%%%%
S1=5.76e14;
Sm=5.75e14; 
S2=5.78e14; % surface de la distance où se trouve le satellite en m2
A=3.16e7;     % un an exprimé en secondes  
T=5.77;	     % durée d'exposition du satellite	
rho=2.5;      % densité des particules en g/cm3

%da=input('paramètre')
% dx est le pas de variation du log de la masse

dx=0.005;

% x1 et x2 sont les bornes entre lesquelles varie le log de la masse du micrométéoride.

x1=-9;
x2=-4.10;

%x est le vecteur log de la masse du micrométéoride (en g)

x=x1:dx:x2;
mb=10.^x;

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

npas=1+(x2-x1)/dx;

%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=(x-log10(pi/6*rho))/3;
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;


%Paramètres du fit polynomial effectué par Brownlee et al.

a=8.8749;
b=-25.439;
c=28.054;
d=-15.536;
e=4.1732;
f=-0.44175;

%Le fit polynomial lui-même et sa dérivée.

alpha=a+b*zz+c*zz.^2+d*zz.^3+e*zz.^4+f*zz.^5;

dalpha=b+2*c*zz+3*d*zz.^2+4*e*zz.^3+5*f*zz.^4;

%flux cumulatif de météorides /m2/s

for j=1:npas
   cnm(j)=(10.^alpha(j))/T/A*10000;
end



%calcul pas à pas du flux différentiel

np(1,npas)=0;

teta=S2*A/1e9; %résultat en milliers de tonnes sur toute la terre

phip(1,npas)=teta*mb(1,npas)*np(1,npas);


dxb=0.5;
npasb=1+(x2-x1)/dxb;
for j=1:npas-1   
   np(j)=cnm(j)-cnm(j+1);
   phip(j)=teta*mb(j+1)*np(j)/dx;
end


%On a ppelle phi le flux différentiel.
eta=S2/T*10000; %eta facteur géométrique

phi0=eta*1.056/3*log(10);

phi=-(phi0.*10.^(alpha+x-9)).*dalpha; %expresion en milliers de tonnes par an pour toute la terre.

%Calcul du flux intégral.

int(1,1)=phi(1,1)*dx;
   
for j=1:npas-1
   int(1,j+1)=int(1,j)+phi(1,j+1)*dx;
end

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

for j=1:npas
   rel(1,j)=int(1,j)/int(1,npas);
   bro(1,j)=40*phi(1,j)/int(1,npas);
 end


x3=[75 125 175 225 275 350 400];
B=[-6.74 -6.70 -6.70 -6.79 -6.96 -7.25 -7.60];
l3=log10(x3);

plot(x,phi)
int(1,npas)