Matlab

function [vol] = monte_carlo_torus(n)
x=4*rand(n,1);
y=4*rand(n,1);
z=4*rand(n,1);

 vol=0;
 for ii=1:n
     if (rand(2,1,1)==1)
         x(ii)=-x(ii);
     end
     if (rand(2,1,1)==1)
          y(ii)=-y(ii);
     end
     if (rand(2,1,1)==1)
           z(ii)=-z(ii);
     end
     cond1=z(ii)^2 + (sqrt(x(ii)^2+y(ii));
     cond2=x(ii)>=1;
     cond2=x(ii)>=-3;
     if (cond1)
         plot3(x(ii),y(ii),z(ii),':g.');
         hold on;
     elseif (cond2 && cond3)
         plot3(x(ii),y(ii),z(ii),':ro');
         xlabel('x axis');
         ylabel('y axis');
         zlabel('z axis');
     elseif (cond1 && cond2 && cond3)
         vol=vol + 1;
     else
           plot3(x(ii),y(ii),z(ii),':b');
     end

AngleDeg = 0:1:360;
AngleRad = AngleDeg *pi / 180;

%Aire totale délimitée
xmin2=-20;
xmax2=20;
ymin2=-20;
ymax2=20;
%xmin3=xmin2/2; %Cela permet au cercle de rester dans la surface
%xmax3=xmax2/2;
%ymin3=ymin2/2;
%ymax3=ymax2/2;
Aire_totale2=(ymax2-ymin2)*(xmax2-xmin2); %Aire totale agrandie pour contenir les deux cercles

%Création des cercles
Abscisse1=-5;
Ordonnee1=0;
Abscisse2=5;
Ordonnee2=0;
%Abscisse1=xmax3+(xmin3-xmax3)*rand;
%Ordonnee1=xmax3+(xmin3-xmax3)*rand;
%Abscisse2=xmax3+(xmin3-xmax3)*rand;
%Ordonnee2=xmax3+(xmin3-xmax3)*rand;
Rayon1=12;
Rayon2=12;
%Rayon1=xmax3*rand;
%Rayon2=xmax3*rand;

XCercle1 = Abscisse1 + Rayon1 * cos(AngleRad);
YCercle1 = Ordonnee1 + Rayon1 * sin(AngleRad);
XCercle2 = Abscisse2 + Rayon2 * cos(AngleRad);
YCercle2 = Ordonnee2 + Rayon2 * sin(AngleRad);
aire_cercle1= pi*Rayon1*Rayon1;
aire_cercle2= pi*Rayon2*Rayon2;

%Lecture du nombre de tirages à effectuer
nb_tirages=0;
while nb_tirages<=0
    nb_tirages=input('Combien de tirages souhaitez-vous effectuer ? \n');
    nb_tirages= abs(nb_tirages);
end

%CALCUL DE LA VALEUR DE L'AIRE A DÉTERMINER 

%Tirage aléatoire des points
PointX = xmax2+(xmin2-xmax2)*rand(1,nb_tirages);
PointY = ymax2+(ymin2-ymax2)*rand(1,nb_tirages);

points_intersection=0;
points_cercle1=0;
points_cercle2=0;
%aire_intersection_reel=0;
inter= 0; %Variable indiquant s'il y a une intersection entre les 2 cercles ou pas
kx=[];
ky=[];

%Détermination des points situés dans l'intersection
for N=1:1:nb_tirages
    if (PointX(N)-Abscisse1)^2+(PointY(N)-Ordonnee1)^2<=Rayon1^2
        points_cercle1=points_cercle1+1;
        kx(points_cercle1)=PointX(N);
        ky(points_cercle1)=PointY(N);
    end
    if (PointX(N)-Abscisse2)^2+(PointY(N)-Ordonnee2)^2<=Rayon2^2
        points_cercle2=points_cercle2+1;
        for m=1:1:(points_cercle1)
            if PointX(N)==kx(m) && PointY(N)==ky(m)
                inter= 1; 
                points_intersection= points_intersection+1;
            end 
        end
    end
end

aire_intersection_estime= (Aire_totale2*points_intersection)/nb_tirages;
%if inter==1
%    aire_intersection_reel= abs(aire_cercle1-aire_cercle2);
%end

%Affichage
figure (2)
plot(XCercle1,YCercle1,'r',XCercle2,YCercle2,'m',PointX,PointY,'b')
axis([xmin2 xmax2 ymin2 ymax2]);