clear

// **************
// *            *
// * Constantes *
// *            *
// **************

norme2=1/%pi; // integrale de exp(-r^2) (fonction de Gauss non normee) sur le plan
R=3; // rayon 3*sigma (sigma=1)
pas=0.01; // pas d'integration


// *************
// *           *
// * Fonctions *
// *           *
// *************


// Fonction de Gauss non normee

function [y]=gauss(x)
   y=exp(-x^2);
endfunction


// Fonction cartesienne du cercle
// de centre O et de rayon r

function x=cercle(y,r)
   if abs(y)>abs(r) then
      x=0;
   else
      x=sqrt(r^2-y^2);
   end
endfunction


// Integrale d'une fonction radiale sur une surface
// sur une bande entre (x1;y-pas/2) et (x2;y+pas/2)

function [somme]=integrale_rad2_ligne(f,x1,x2,y1)

   somme=integrate('f(sqrt(x^2+y1^2))','x',x1,x2);
   somme=somme*pas;

endfunction


// Pour une case du tableau :
// Integrale de la gaussienne normee
// sur un disque de rayon r0 et de centre (x0,0)

function [proba]=case(x0,r0)
   proba=0;
// sur un demi disque, sans compter deux fois la bande centrale
   for y=0+pas/2:pas:r0
     x1=cercle(y,r0);
     proba=proba+integrale_rad2_ligne(gauss,-x1-x0,x1-x0,y);
   end
   proba = round(2*norme2*proba*100)
endfunction


// ***********************
// *                     *
// * Programme principal *
// *                     *
// ***********************


// Variables ajustables

M=zeros(11,10); // tableau initial rempli de 0

// Calcul

for i=1:11
   for j=1:10
      // ligne i : l'ecart vaut 0,1x(i-1)xR (donc de 0 a R)
      // colonne j : le rayon du disque vaut (1,1-0,1j)xR (donc de R a 0,1xR)
      M(i,j)=case((i-1)*0.1*R,R*(1.1-0.1*j));
   end
end

M