// Initialisation

clear;

// constantes

n = 1e6; // nombre d'événements

Resp = 269; // espérance de la résistance
Rect = 17; // écart type ""
Sesp = 214; // espérance de la contrainte
Sect = 14 // écart type ""

sigma = 172:320; // ensemble de définition
taille = size(sigma, "c");
sigma3d = 172:5:320;
taille3d = size(sigma3d, "c");
sigmaniv = 172:2:320;
tailleniv = size(sigmaniv, "c");

// Fonctions

densiteR = pdfnormal(sigma, Resp, Rect);
densiteS = pdfnormal(sigma, Sesp, Sect);

densiteR3d = pdfnormal(sigma3d, Resp, Rect);
densiteS3d = pdfnormal(sigma3d, Sesp, Sect);

densiteRniv = pdfnormal(sigmaniv, Resp, Rect);
densiteSniv = pdfnormal(sigmaniv, Sesp, Sect);

densiteegalite = densiteR.*densiteS;
for i=1:taille3d
    densite3d(i, :) = densiteR3d(i)*densiteS3d;
    // x = R, y = S
end

for i=1:tailleniv
    densiteniv(i, :) = densiteRniv(i)*densiteSniv;
end

scf(0);
clf;

subplot(2,1,1)
plot(sigma, 100*densiteS, "r");
plot(sigma, 100*densiteR, "b");

xtitle("Méthode contrainte-résistance",...
    "Contrainte (MPa)", "Densité de probabilité f (%/MPa)");
legend("Contrainte S", "Résistance R");

subplot(2,1,2)
plot(sigma, 100*densiteegalite, "k");

xtitle("Probabilité d''égalité",...
    "Contrainte (MPa)", "Densité de probabilité f (%/MPa)");

scf(1);
clf;
plot3d1(sigma3d, sigma3d, 1e5*densite3d)
h3d = gce();
h3d.thickness = 0;
cmap = 1-0.7*graycolormap(100);
xset("colormap", cmap);

param3d(172*ones(sigma), sigma, 5e3*densiteS);

param3d(sigma, 172*ones(sigma), 5e3*densiteR)

param3d(sigma, sigma, 1e5*densiteegalite, 135)
a = gca();
a.tight_limits="on";

xtitle("Probabilité de rupture", "Résistance R (MPa)", "Contrainte S (MPa)",...
    "Densité de probabilité (u.a.)")

R = 269 + 17*rand(n, 1, "normal"); // résistance
S = 214 + 14*rand(n, 1, "normal"); // contrainte
Mbool = S>R; // comparaison des valeurs (matrice de booléens)
moyenne = sum(Mbool)/n; // les "vrai" comptent pour 1 et les "faux" pour 0
disp(string(100*moyenne)+" %") // affichage du résultat

scf(2)
clf;
grayplot(sigmaniv, sigmaniv, 1e5*densiteniv)
xpoly([172, 320], [172, 320]);
xset("colormap", cmap);
xtitle("Densité de probabilité",...
    "Résistance R (MPa)", "Contrainte S (MPa)")