clear all;
close all;
clc;

% valec lopaty + geometria
D1=57.2; %priemer valca
p1=24; %tlak vo valci v MPa
S1=(pi()*D1^2)/4;
Fv1=p1*S1;
a=??;
b=??;
c=??;
d=??;
f=??;
h=??;
m=??; %roztiahnuty valec
k=m-??; %stiahnuty valec
l=??;
delta=??;
i=??;
e=0;
g=k+e;

% valec nasady 29 Mpa
r=??; %rameno
D2=??; %priemer valca
p2=29; %tlak vo valci v MPa
S2=(pi()*D2^2)/4;
Fv2=p2*S2; %maximalna sila pred prepustenim
Mo2=Fv2*r; %Nmm - moment od prepustacieho valca

%pomocne premenne
z=0;
Fk_max=0;
Fk_max_p=0;
Mo_max=0;

%vytvorenie poli buniek
x=[];
M_ohyb_b=[];
M_ohyb_p=[];
Fk_b=[];
Fk_p=[];

while g<=m

%pocitadlo
z=z+1;
    
% vypocet geometrie
u=sqrt(h^2+f^2);
alfa1=atand(f/h);
alfa2=acosd((u^2+g^2-d^2)/(2*u*g));
alfa3=alfa1+alfa2;
alfa=90-alfa3;
beta1=90-alfa1;

%podmienka pre uhol beta2
if (d^2+u^2>g^2)
    beta2=asind(g*(sind(alfa2))/d);
else 
    beta2p=asind(g*(sind(alfa2))/d);
    rozdiel=90-beta2p;
    beta2=90+rozdiel;
end

beta3=beta1+beta2;
beta4=180-beta3;
v=sqrt(a^2+d^2-2*a*d*cosd(beta4));
gama1=acosd((v^2+c^2-b^2)/(2*v*c));
gama2=asind(a*(sind(beta4))/v);
gama3=gama1+gama2;

%podmienka pre uhol gama4
if (b^2+c^2<v^2)
    gama4p=asind(v*(sind(gama1))/b);
    rozdiel=90-gama4p;
    gama4=90+rozdiel;
else
    gama4=asind(v*(sind(gama1))/b);
end

gama5=gama3-beta3;
delta1=gama4+gama5;
delta2=delta1-delta;
delta3=90-delta2;
delta4=180-delta1;

%vysledna matica reakcii
A=[1 0 -1 0 0 0 0 0 0
    0 1 0 -1 0 0 0 0 0
    d*sind(beta3) d*cosd(beta3) 0 0 0 0 0 0 0
    0 0 1 0 -1 0 0 0 0
    0 0 0 1 0 -1 0 0 0
    0 0 0 0 c*sind(gama5) -c*cosd(gama5) 0 0 0
    0 0 0 0 1 0 -1 0 -cosd(delta3)
    0 0 0 0 0 1 0 -1 sind(delta3)
    0 0 0 0 -b*sind(delta4) -b*cosd(delta4) 0 0 l];

%vektor pravych stran
y=[-Fv1*cosd(alfa);Fv1*sind(alfa);0;0;0;0;0;0;0];

%vypocet reakcii
X=inv(A)*y;
Fk=X(9); %rypna sila pred prepustenim

%aktualny zdvih valca
x(z)=g;

%vypocet ohyboveho momentu pred prepustenim v Nmm
M_ohyb_b(z)=Fk*(l+i*cosd(delta2));

%podmienka pre moment zmenseny pri prepusteni ventilu
if M_ohyb_b(z)>Mo2
    M_ohyb_p(z)=Mo2;
else
    M_ohyb_p(z)=M_ohyb_b(z);
end

%rypna sila bez prepustenia
Fk_b(z)=Fk;

%rypna sila po prepusteni
if M_ohyb_b(z)>Mo2
    Fk_p(z)=Mo2/(l+i*cosd(delta2));
else
    Fk_p(z)=Fk_b(z);
end

%poloha pre maximalny ohybovy moment bez prepustenia
if Mo_max<M_ohyb_b(z)
Mo_max=M_ohyb_b(z);
poloha_ohyb=g;
Fk_ohyb=Fk;
end

%poloha maximalnej rypnej sily bez prepustenia
if Fk_max<Fk
Fk_max=Fk;
poloha=g;
end

%maximalna rypna sila s uvazovanim prepustenia
if Fk_max_p<Fk_p(z)
Fk_max_p=Fk_p(z);
poloha_p=g;
Ax=X(1);
Ay=X(2);
Bx=X(3);
By=X(4);
Cx=X(5);
Cy=X(6);
Dx=X(7);
Dy=X(8);
alfa_max=alfa;
delta3_max=delta3;
end

e=e+0.01; %krok

g=k+e;

end

%vypis Fk_max pred a po prepusteni, Mo_max a ich polohy
fprintf('Najväèšia rypná sila je %.2f N pri zdvihu valca %.2f mm.\n\n', Fk_max,poloha);

fprintf('Najväèší ohybový moment je %.2f Nm pri zdvihu valca %.2f mm.\n Rypná sila pri najväèšom ohybovom momente je %.2f N.\n', Mo_max/1000,poloha_ohyb,Fk_ohyb);

fprintf('\n\nNajväèšia rypná sila s uvažovaním prepustenia je %.2f N pri zdvihu valca %.2f mm.\n\n', Fk_max_p,poloha_p);

subplot(2,1,1); 
 y1=M_ohyb_b/1000; y2=M_ohyb_p/1000;
 plot(x,y1,'r-')
 hold on
 plot(x,y2,'b-')
 hold off
 grid on
 xlabel('Zdvih valca [mm]') %popis x-ovej osy
 ylabel('M ohyb[Nm]')
 title('Závislos� ohybového momentu na zdvihu valca') %nadpis grafu
 legend('bez prepustenia','s prepustením') %popis krivek nachazejicich se v grafu
 
 subplot(2,1,2); 
 y1=Fk_b; y2=Fk_p;
 plot(x,y1,'r-')
 hold on
 plot(x,y2,'b-')
 hold off
 grid on
 xlabel('Zdvih valca [mm]') %popis x-ovej osy
 ylabel('Rypná sila[N]')
 title('Závislos� rypnej sily na zdvihu valca') %nadpis grafu
 legend('bez prepustenia','s prepustením') %popis krivek nachazejicich se v grafu