clc;
clear all; 
close all;
p1=20; 
p2=15; 
r1=30;
r2=50;
l=100;
E=390000;
mi=0.24;
m=7.19;
sigma0=473.8;
lambda=(mi*E)/((1+mi)*(1-2*mi));
G=E/(2*(1+mi));


%riešenie sústavy lin. rovníc
Ax=[1,-1/(r1^2),lambda;1,-1/(r2^2),lambda;2*mi,0,(2*G+lambda)];
b=[0;-p1;-p2];
x=inv(Ax)*b;
A=x(1,:);
B=x(2,:);
Ez=x(3,:);

rx=41.7647;

%vykreslenie priebehu napätí
r=linspace (r1,r2,52);
f=figure(1);
sigmar=A-(B./(r.^2))+(lambda*Ez);
sigmat=A+(B./(r.^2))+(lambda*Ez);
sigmaz=(2*mi*A+(2*G+lambda)*Ez)*ones(1,length(r));
hold on;
grid on;
plot(r,sigmar,'-r','LineWidth',1.5);
plot(r,sigmat,'-b','LineWidth',1.5);
plot(r,sigmaz,'g','LineWidth',1.5);
legend('radiálne napätie','dotyčnicové napätie','napätie v smere z','FontSize',12,'FontName','Times New Roman');
set(f,'color',[1 1 1],'units','centimeters','position',[0 0 16 8]);
xlabel('r [mm]','FontSize',12,'FontName','Times New Roman');
ylabel('napätie [MPa]','FontSize',12,'FontName','Times New Roman');
title('Variant B','FontSize',12,'FontName','Times New Roman');

%výpočet rizika porušenia súčasti (prvý interval)
sigmar=@(r)  A-(B./(r.^2))+(lambda*Ez);
sigmat=@(r)  A+(B./(r.^2))+(lambda*Ez);
sigmaz=-p2;

sigma11=@(r) sigmar(r);
sigma21=sigmaz;
sigma31=@(r) sigmat(r);
sigman1=@(r,fi,theta) (sigma11(r).*(cos(fi)).^2+sigma21.*(sin(fi)).^2).*(sin(theta)).^2+sigma31(r).*(cos(theta)).^2;
sigmae1=@(r,fi,theta) (sigman1(r,fi,theta)+abs(sigman1(r,fi,theta)))./2;
f1=@(theta,fi,r) (2.*m+1).*l./(sigma0.^m).*r.*sin(theta).*sigmae1(r,fi,theta).^m;
B1=integral3(f1,0,pi/2,0,2*pi,r1,rx)

%výpočet rizika porušenia súčasti (druhý interval)
sigma12=sigmaz;
sigma22=@(r) sigmar(r);
sigma32=@(r) sigmat(r);
sigman2=@(r,fi,theta) (sigma12.*(cos(fi)).^2+sigma22(r).*(sin(fi)).^2).*(sin(theta)).^2+sigma32(r).*(cos(theta)).^2;
sigmae2=@(r,fi,theta) (sigman2(r,fi,theta)+abs(sigman2(r,fi,theta)))./2;
f2=@(theta,fi,r) (2.*m+1).*l./(sigma0.^m).*r.*sin(theta).*sigmae2(r,fi,theta).^m;
B2=integral3(f2,0,pi/2,0,2*pi,rx,r2)

%výpočet pravdepodobnosti porušenia súčasti
B=B1+B2;
Pf=(1-exp(-B))*100;



%zobrazenie výsledku
disp('Riziko porušenia súčasti B je')
B
disp('Pravdepodobnosť porušenia súčasti Pf v percentách je')
Pf