clc
clear all

%% VSTUPNÍ PARAMETRY
%a = 0:(2*pi)/359:2*pi;                              % úhel natočení klikového hřídele 1x360
%a = 0:(2*pi)/719:2*pi;                              % úhel natočení klikového hřídele 1x720
a = 0:(4*pi)/1079:4*pi;                               % úhel natočení klikového hřídele 1x1080
dt0= 233*pi/180;                                          % úhel odklonu 2. vedlejší ojniční kliky
dt = 240*pi/180;                                           % úhel odklonu osy 2. vedlejšího válce
rk = 43.45;                                                      % rk - poloměr zalomení klikové hřídele
rls = 51.2;                                                        % délka vedlejší ojniční kliky
rlt = rls;                                                            % délka vedlejší ojniční kliky
ls =86.42;                                                        % délka 2. vedlejší ojnice 
lt = ls;                                                                % délka 1. vedlejší ojnice 
D=76.51;                                                         % vrtání
Zpist = 29.27;                                                 % kompresní výška pístu
lm = 138;                                                         % délka hlavní ojnice
Em = 10.3                                                        % kompresní poměr hlavního válce
%c = 360;                                                         % krok
%c = 720;                                                         % krok
c =1080;                                                           % krok
ds = 120*pi/180;                                           % úhel odklonu osy 1. vedlejšího válce
ds0 = 127*pi/180;                                         % úhel odklonu 1. vedlejší ojniční kliky od osy hlavní ojnice
fis = ds0-ds;                                                     % úhel udávající rozdíl mezi ds0 a ds
n = 3000;                                                         % otáčky za minutu
w=2*pi*(n/60);                                             % úhlová rychlost
fit = dt0 - dt;                                                    % úhel udávající rozdíl mezi dt0 a dt
patm = 0.1;                                                      % atmosferický tlak [Mpa]
patm = patm*10^6;                                      % atmosferický tlak [Pa]

%% VÝPOČET POLOHY KLIKOVÉ HŘÍDELE

for i = 1:length(a)
Bm(i) = asin((rk*sin(a(i)))/lm);                  % úhel odklonu hlavní ojnice
sm(i)=rk*cos(a(i))+lm*cos(Bm(i));           % dráha hlavního pístu
end

% Výpočet velikosti spalovacího prostoru hlavního válce
Zm = max(sm) - min(sm);                           % Zdvih hlavního válce
Vzm = Zm*pi*(D^2/4);                                % Zdvihový objem hlavního válce
Vkm = Vzm/(Em-1);                                      % Kompresní objem hlavního válce
Zkm = (4*Vkm)/((D^2)*pi);                        % Výška kompresního objemu
Z = Zkm + Zpist + lm + rk;                             % Vzdálenost osy otáčení KH od stropu spalovacího prostoru

% První vedlejší válec
for i = 1:length(a)
    B2s(i) = Bm(i) - fis;                                     % Úhel odklonu 1. vedlejší kliky od osy 1. vedlejšího válce
    Bcs(i) = ds - a(i);                                         % úhel natočení klikové hřídele vůči ose prvního vedlejšího válce
    Bs(i) = asin((rk*sin(Bcs(i))+(rls*sin(B2s(i))))/ls);  % Úhel odklonu první vedlejší ojnice vůči ose prvního vedlejšího válce
    ss(i) = rk*cos(Bcs(i))+rls*cos(B2s(i))+ls*cos(Bs(i));    % Dráha pístu 1. vedlejšího válce
end

% Výpočet velikosti spalovacího prostoru prvního vedlejšího válce
Zs = max(ss)-min(ss);                                   % Zdvih pístu 1. vedlejšího válce
Zks = Z - max(ss) - Zpist;                              % Výška kompresního objemu 1. vedlejšího válce
Vks = Zks*pi*(D^2)/4;                                 % Kompresní objem 1. vedlejšího válce
Vzs = Zs*pi*(D^2)/4;                                    % Zdvihový objem 1. vedlejšího válce
Es = (Vzs+Vks)/Vks                                       % Kompresní poměr 1. vedlejšího válce

% Druhý vedlejší válec
for i = 1:length(a)
    B2t(i) = Bm(i) - fit;                                      % Úhel odklonu 2. vedlejší kliky od osy 2. vedlejšího válce
    Bct(i) = dt - a(i);                                          % úhel natočení KH vůči ose 2. válce
    Bt(i) = asin((rk*sin(Bct(i))+(rlt*sin(B2t(i))))/lt);    % Úhel odklonu 2. vedlejší ojnice vůči ose 2. vedlejší ojnice
    st(i) = rk*cos(Bct(i))+rlt*cos(B2t(i))+lt*cos(Bt(i)); % Dráha pístu 2. vedlejšího válce
end

% Výpočet velikosti spalovacího prostoru druhého vedlejšího válce
Zt = max(st)-min(st);                                    % Zdvih pístu 2. vedlejšího válce
Zkt = Z - max(st) - Zpist;                               % Výška kompresního objemu 2. vedlejšího válce
Vkt = Zkt*pi*(D^2)/4;                                  % Kompresní objem 2. vedlejšího válce
Vzt = Zt*pi*(D^2)/4;                                     % Zdvihový objem 2. vedlejšího válce
Et = (Vzt+Vkt)/Vkt                                         % Kompresní poměr 2. vedlejšího válce

%plot(a,sm,a,st,a,ss);

%% VÝPOČTY REDUKOVANÝCH HMOTNOSTÍ
lvo = 86.42;                                                     %délka vedlejší ojnice [mm] 
lrvo = lvo - 39.3;                                            %vzdálenost těžiště od počátku u vedlejší ojnice [mm]
mvo = 0.155;                                                  % hmotnost vedlejší ojnice [kg]
mpist = 285/1000;                                        % hmotnost pístu [Kg]
lvo = lvo/1000;                                               % délka vedlejší ojnice [m]
lrvo = lrvo/1000;                                           % vzdálenost těžiště od počátku u vedlejší ojnice [m]

mA = mvo*((lvo-lrvo)/lvo);                        %redukovana hmotnost posuvnych casti vedlejsi ojnice [Kg]
mB = mvo*(lrvo/lvo);                                   %redukovana hmotnost rotacnich casti vedlejsi ojnice [Kg]

mA = mA*1000;                                            %redukovana hmotnost posuvnych casti vedlejsi ojnice [g]
mB = mB * 1000;                                           %redukovana hmotnost rotacnich casti vedlejsi ojnice [g]


lho = 138;                                                         % délka hlavní ojnice [mm]
lrho = 5.7;                                                         % vzdálenost těžiště od počátku u hlavní ojnice [mm]
mho = 1.025;                                                   % hmotnost sestavy hlavní ojnice s dalšími díly [kg]

lho = lho/1000;                                               % délka hlavní ojnice [m]
lrho = lrho/1000;                                            % vzdálenost těžiště od počátku u hlavní ojnice [m]

mC = mho*((lho-lrho)/lho);                        %redukovana hmotnost rotacnich casti hlavni ojnice [Kg]
mD = mho*(lrho/lho);                                 %redukovana hmotnost posuvnych casti hlavni ojnice [Kg]

mC = mC*1000;                                             %redukovana hmotnost rotacnich casti hlavni ojnice [g]
mD = mD * 1000;                                          %redukovana hmotnost posuvnych casti hlavni ojnice [g]

mzal = 0.763;                                                 % hmotnost zalomení klikové hřídele [Kg]
lr = 42.1 ;                                                         %poloha težiště zalomení klikové hřídele [mm]
R = 43.45 ;                                                       % poloměr zalomení klikové hřídele [mm]

lr = lr/1000;                                                    %poloha težiště zalomení klikové hřídele [m]
R=R/1000;                                                       % poloměr zalomení klikové hřídele [m]

mE = mzal*(lr/R);                                          % redukovaná hmotnost zalomení klikové hřídele [Kg]

%% SÍLA OD TLAKU PLYNŮ
x= 0:720/1079:720;                                     % rozložení dvou otáček klikové hřídele na 1080 kroků ve stupních
atlak = 0:(2*pi)/1079:2*pi;                        % rozložení dvou otáček klikové hřídele na 1080 kroků v radiánech
p = importdata('otto3000rpm.dat');       % importování tlaků
pm = importdata('Tlak_m_valec.txt');     % import tlaku v hlavním válci [Pa]
ps = importdata('Tlak_s_valec.txt');         % import tlaku v prvním vedlejším válci [Pa]
pt = importdata('Tlak_t_valec.txt');          % import tlaku ve druhém vedlejším válci [Pa]

for i = 1:length(atlak)
    
Fpm(i) = ((pm(i) - patm)*((pi*(D/1000)^2)/4));         % síla od tlaku plynů v hlavním válci [N]
Fps(i) = ((ps(i) - patm)*((pi*(D/1000)^2)/4));             % síla od tlaku plynů v prvním vedlejším válci [N]
Fpt(i) = ((pt(i) - patm)*((pi*(D/1000)^2)/4));              % síla od tlaku plynů ve druhém vedlejším válci [N]

end

%plot(x,pm,x,ps,x,pt)
%plot(x,Fpm,x,Fps,x,Fpt)

%% RYCHLOST A ZRYCHLENÍ PÍSTU

%Fourierovi členy pro rychlost a zrychlení pístu
for i = 1:length(a)
Ax1(i) = sm(i)*sin(1*a(i));
Ax2(i) = sm(i)*sin(2*a(i));
Ax3(i) = sm(i)*sin(3*a(i));
Ax4(i) = sm(i)*sin(4*a(i));
Ay1(i) = st(i)*sin(1*a(i));
Ay2(i) = st(i)*sin(2*a(i));
Ay3(i) = st(i)*sin(3*a(i));
Ay4(i) = st(i)*sin(4*a(i));
Az1(i) = ss(i)*sin(1*a(i));
Az2(i) = ss(i)*sin(2*a(i));
Az3(i) = ss(i)*sin(3*a(i));
Az4(i) = ss(i)*sin(4*a(i));
Bx1(i) = sm(i)*cos(1*a(i));
Bx2(i) = sm(i)*cos(2*a(i));
Bx3(i) = sm(i)*cos(3*a(i));
Bx4(i) = sm(i)*cos(4*a(i));
By1(i) = st(i)*cos(1*a(i));
By2(i) = st(i)*cos(2*a(i));
By3(i) = st(i)*cos(3*a(i));
By4(i) = st(i)*cos(4*a(i));
Bz1(i) = ss(i)*cos(1*a(i));
Bz2(i) = ss(i)*cos(2*a(i));
Bz3(i) = ss(i)*cos(3*a(i));
Bz4(i) = ss(i)*cos(4*a(i));
end
Am1 = sum(Ax1) * (2/c);
Am2 = sum(Ax2) * (2/c);
Am3 = sum(Ax3) * (2/c);
Am4 = sum(Ax4) * (2/c);
At1 = sum(Ay1) * (2/c);
At2 = sum(Ay2) * (2/c);
At3 = sum(Ay3) * (2/c);
At4 = sum(Ay4) * (2/c);
As1 = sum(Az1) * (2/c);
As2 = sum(Az2) * (2/c);
As3 = sum(Az3) * (2/c);
As4 = sum(Az4) * (2/c);
Bm1 = sum(Bx1) * (2/c);
Bm2 = sum(Bx2) * (2/c);
Bm3 = sum(Bx3) * (2/c);
Bm4 = sum(Bx4) * (2/c);
Bt1 = sum(By1) * (2/c);
Bt2 = sum(By2) * (2/c);
Bt3 = sum(By3) * (2/c);
Bt4 = sum(By4) * (2/c);
Bs1 = sum(Bz1) * (2/c);
Bs2 = sum(Bz2) * (2/c);
Bs3 = sum(Bz3) * (2/c);
Bs4 = sum(Bz4) * (2/c);

%Rychlost a zrychlení pístů
for i= 1:length(a)
    vm(i)=(Am1*cos(a(i))-Bm1*sin(a(i))+2*Am2*cos(2*a(i))-2*Bm2*sin(2*a(i)))*w;%+3*Am3*cos(3*a(i))-3*Bm3*sin(3*a(i))+4*Am4*cos(4*a(i))-4*Bm4*sin(4*a(i)))*(w);
     vs(i)=(As1*cos(a(i))-Bs1*sin(a(i))+2*As2*cos(2*a(i))-2*Bs2*sin(2*a(i)))*w;%+3*As3*cos(3*a(i))-3*Bs3*sin(3*a(i))+4*As4*cos(4*a(i))-4*Bs4*sin(4*a(i)))*(w);
     vt(i)=(At1*cos(a(i))-Bt1*sin(a(i))+2*At2*cos(2*a(i))-2*Bt2*sin(2*a(i)))*w;%+3*At3*cos(3*a(i))-3*Bt3*sin(3*a(i))+4*At4*cos(4*a(i))-4*Bt4*sin(4*a(i)))*(w);
     am(i)=(-1)*(-Am1*sin(a(i))-Bm1*cos(a(i))-4*Am2*sin(2*a(i))-4*Bm2*cos(2*a(i)))*(w^2);%-9*Am3*sin(3*a(i))-9*Bm3*cos(3*a(i))-16*Am4*sin(4*a(i))-16*Bm4*cos(4*a(i)))*(w^2);
      as(i)=(-As1*sin(a(i))-Bs1*cos(a(i))-4*As2*sin(2*a(i))-4*Bs2*cos(2*a(i)))*(w^2);%-9*As3*sin(3*a(i))-9*Bs3*cos(3*a(i))-16*As4*sin(4*a(i))-16*Bs4*cos(4*a(i)))*(w^2);
      at(i)=(-At1*sin(a(i))-Bt1*cos(a(i))-4*At2*sin(2*a(i))-4*Bt2*cos(2*a(i)))*(w^2);%-9*At3*sin(3*a(i))-9*Bt3*cos(3*a(i))-16*At4*sin(4*a(i))-16*Bt4*cos(4*a(i)))*(w^2);
end

vm=vm/1000;                                                % Rychlost hlavního pístu [m/s]
vs=vs/1000;                                                    % Rychlost prvního vedlejšího pístu [m/s]
vt=vt/1000;                                                     % Rychlost prvního vedlejšího pístu [m/s]
am=am/1000;                                                % Zrychlení hlavního pístu [m/s^2]
as=as/1000;                                                    % Zrychlení prvního vedlejšího pístu [m/s^2]
at=at/1000;                                                    % Zrychlení druhého vedlejšího pístu [m/s^2]
%plot(a,vm,a,vs,a,vt)   
%plot(a,am,a,as,a,at)

%% ÚHLOVÁ RYCHLOST A ZRYCHLENÍ OJNIC

% Fourierovi členy pro úhlovou rychlost a zrychlení ojnic
for i= 1:length(a)
Aa1(i) = Bm(i)*sin(1*a(i));                          % sinusový fourierův člen 
Aa2(i) = Bm(i)*sin(2*a(i));
Aa3(i) = Bm(i)*sin(3*a(i));
Aa4(i) = Bm(i)*sin(4*a(i));
Ab1(i) = Bt(i)*sin(1*a(i));
Ab2(i) = Bt(i)*sin(2*a(i));
Ab3(i) = Bt(i)*sin(3*a(i));
Ab4(i) = Bt(i)*sin(4*a(i));
Ac1(i) = Bs(i)*sin(1*a(i));
Ac2(i) = Bs(i)*sin(2*a(i));
Ac3(i) = Bs(i)*sin(3*a(i));
Ac4(i) = Bs(i)*sin(4*a(i));
Ba1(i) = Bm(i)*cos(1*a(i));
Ba2(i) = Bm(i)*cos(2*a(i));
Ba3(i) = Bm(i)*cos(3*a(i));
Ba4(i) = Bm(i)*cos(4*a(i));
Bb1(i) = Bt(i)*cos(1*a(i));
Bb2(i) = Bt(i)*cos(2*a(i));
Bb3(i) = Bt(i)*cos(3*a(i));
Bb4(i) = Bt(i)*cos(4*a(i));
Bc1(i) = Bs(i)*cos(1*a(i));
Bc2(i) = Bs(i)*cos(2*a(i));
Bc3(i) = Bs(i)*cos(3*a(i));
Bc4(i) = Bs(i)*cos(4*a(i));
end

Amm1 = sum(Aa1) * (2/c);
Amm2 = sum(Aa2) * (2/c);
Amm3 = sum(Aa3) * (2/c);
Amm4 = sum(Aa4) * (2/c);
Att1 = sum(Ab1) * (2/c);
Att2 = sum(Ab2) * (2/c);
Att3 = sum(Ab3) * (2/c);
Att4 = sum(Ab4) * (2/c);
Ass1 = sum(Ac1) * (2/c);
Ass2 = sum(Ac2) * (2/c);
Ass3 = sum(Ac3) * (2/c);
Ass4 = sum(Ac4) * (2/c);
Bmm1 = sum(Ba1) * (2/c);
Bmm2 = sum(Ba2) * (2/c);
Bmm3 = sum(Ba3) * (2/c);
Bmm4 = sum(Ba4) * (2/c);
Btt1 = sum(Bb1) * (2/c);
Btt2 = sum(Bb2) * (2/c);
Btt3 = sum(Bb3) * (2/c);
Btt4 = sum(Bb4) * (2/c);
Bss1 = sum(Bc1) * (2/c);
Bss2 = sum(Bc2) * (2/c);
Bss3 = sum(Bc3) * (2/c);
Bss4 = sum(Bc4) * (2/c);

%Úhlová rychlost a zrychlení ojnic
for i= 1:length(a)
    wm(i)=(Amm1*cos(a(i))-Bmm1*sin(a(i))+2*Amm2*cos(2*a(i))-2*Bmm2*sin(2*a(i))+3*Amm3*cos(3*a(i))-3*Bmm3*sin(3*a(i))+4*Amm4*cos(4*a(i))-4*Bmm4*sin(4*a(i)))*(w);
     ws(i)=(Ass1*cos(a(i))-Bss1*sin(a(i))+2*Ass2*cos(2*a(i))-2*Bss2*sin(2*a(i))+3*Ass3*cos(3*a(i))-3*Bss3*sin(3*a(i))+4*Ass4*cos(4*a(i))-4*Bss4*sin(4*a(i)))*(w);
     wt(i)=(Att1*cos(a(i))-Btt1*sin(a(i))+2*Att2*cos(2*a(i))-2*Btt2*sin(2*a(i))+3*Att3*cos(3*a(i))-3*Btt3*sin(3*a(i))+4*Att4*cos(4*a(i))-4*Btt4*sin(4*a(i)))*(w);
     Eum(i)=(-Amm1*sin(a(i))-Bmm1*cos(a(i))-4*Amm2*sin(2*a(i))-4*Bmm2*cos(2*a(i))-9*Amm3*sin(3*a(i))-9*Bmm3*cos(3*a(i))-16*Amm4*sin(4*a(i))-16*Bmm4*cos(4*a(i)))*(w^2);
      Eus(i)=(-Ass1*sin(a(i))-Bss1*cos(a(i))-4*Ass2*sin(2*a(i))-4*Bss2*cos(2*a(i))-9*Ass3*sin(3*a(i))-9*Bss3*cos(3*a(i))-16*Ass4*sin(4*a(i))-16*Bss4*cos(4*a(i)))*(w^2);
      Eut(i)=(-Att1*sin(a(i))-Btt1*cos(a(i))-4*Att2*sin(2*a(i))-4*Btt2*cos(2*a(i))-9*Att3*sin(3*a(i))-9*Btt3*cos(3*a(i))-16*Att4*sin(4*a(i))-16*Btt4*cos(4*a(i)))*(w^2);
end

%plot(a,wm,a,ws,a,wt)
%plot(a,Em,a,Es,a,Et)

%% SETRVAČNÉ SÍLY POSUVNÝCH ČÁSTÍ
mA = mA/1000;                                         %redukovana hmotnost posuvnych casti vedlejsi ojnice [Kg] 
mD = mD/1000;                                         %redukovana hmotnost posuvnych casti hlavni ojnice [Kg]

mph = mpist + mD;                                       % hmotnost posuvnych casti hlavniho valce
mpv = mpist + mA;                                        % hmotnost posuvnych casti vedlejsiho valce


for i=1:length(a)
    am1(i)=((-Am1*sin(a(i))-Bm1*cos(a(i)))*(w^2))/1000;     % první harmonická složka zrychlení hlavního pístu
    Fsm1(i) = (-mph)*am1(i);                        % Setrvačná síla posuvných částí prvního řádu hlavního válce
    as1(i)=((-As1*sin(a(i))-Bs1*cos(a(i)))*(w^2))/1000;    % první harmonická složka zrychlení prvního vedlejšího pístu 
    Fss1(i) = (-mpv)*as1(i);                            % Setrvačná síla posuvných částí prvního řádu prvního vedlejšího válce
    at1(i)=((-At1*sin(a(i))-Bt1*cos(a(i)))*(w^2))/1000;    % první harmonická složka zrychlení druhého vedlejšího pístu
    Fst1(i) = (-mpv)*at1(i);                            % Setrvačná síla posuvných částí prvního řádu druhého vedlejšího válce
    am2(i)=((-4*Am2*sin(2*a(i))-4*Bm2*cos(2*a(i)))*(w^2))/1000;   % druhá harmonická složka zrychlení hlavního pístu
    Fsm2(i) = (-mph)*am2(i);                       % Setrvačná síla posuvných částí druhého řádu hlavního válce
    as2(i)=((-4*As2*sin(2*a(i))-4*Bs2*cos(2*a(i)))*(w^2))/1000;   % druhá harmonická složka zrychlení prvního vedlejšího pístu
    Fss2(i) = (-mpv)*as2(i);                           % Setrvačná síla posuvných částí druhého řádu prvního vedlejšího válce
    at2(i)=((-4*At2*sin(2*a(i))-4*Bt2*cos(2*a(i)))*(w^2))/1000;   % druhá harmonická složka zrychlení druhého vedlejšího pístu
    Fst2(i) = (-mpv)*at2(i);                            % Setrvačná síla posuvných částí druhého řádu druhého vedlejšího válce
    Fsm(i) = (-mph)*am(i);                            % Setrvačná síla posuvých částí hlavního válce celková
    Fss(i) = (-mpv)*as(i);                                % Setrvačná síla posuvých částí prvního vedlejšího válce celková
    Fst(i) = (-mpv)*at(i);                                 % Setrvačná síla posuvých částí druhého vedlejšího válce celková
end

Fsx1m = Fsm1*sin(0);                                  % X-ová složka setrvačné síly posuvných částí prvního řádu hlavního válce
Fsx2m = Fsm2*sin(0);                                  % X-ová složka setrvačné síly posuvných částí druhého řádu hlavního válce
Fsy1m = Fsm1*cos(0);                                 % Y-ová složka setrvačné síly posuvných částí prvního řádu hlavního válce
Fsy2m = Fsm2*cos(0);                                 % Y-ová složka setrvačné síly posuvných částí druhého řádu hlavního válce
Fsx1s = Fss1*sin(ds);                                   % X-ová složka setrvačné síly posuvných částí prvního řádu prvního vedlejšího válce
Fsx2s = Fss2*sin(ds);                                   % X-ová složka setrvačné síly posuvných částí druhého řádu prvního vedlejšího válce
Fsy1s = Fss1*cos(ds);                                  % Y-ová složka setrvačné síly posuvných částí prvního řádu prvního vedlejšího válce
Fsy2s = Fss2*cos(ds);                                  % Y-ová složka setrvačné síly posuvných částí druhého řádu prvního vedlejšího válce
Fsx1t = Fst1*sin(dt);                                     % X-ová složka setrvačné síly posuvných částí prvního řádu druhého vedlejšího válce
Fsx2t= Fst2*sin(dt);                                      % X-ová složka setrvačné síly posuvných částí druhého řádu druhého vedlejšího válce
Fsy1t = Fst1*cos(dt);                                    % Y-ová složka setrvačné síly posuvných částí prvního řádu druhého vedlejšího válce
Fsy2t = Fst2*cos(dt);                                    % Y-ová složka setrvačné síly posuvných částí druhého řádu druhého vedlejšího válce

Fsx1 = Fsx1m + Fsx1s + Fsx1t;                   % X-ová složka výsledné setrvačné síly prvního řádu
Fsx2 = Fsx2m + Fsx2s + Fsx2t;                   % X-ová složka výsledné setrvačné síly druhého řádu
Fsy1 = Fsy1m + Fsy1s + Fsy1t;                  % Y-ová složka výsledné setrvačné síly prvního řádu
Fsy2 = Fsy2m + Fsy2s + Fsy2t;                  % Y-ová složka výsledné setrvačné síly druhého řádu

%plot(Fsx1,Fsy1)
%plot(Fsx2,Fsy2)

%% ODSTŘEDIVÁ SÍLA ROTAČNÍCH ČÁSTÍ
mkloz = 0.033;                                               %hmotnost hlavního kluzného ložiska [Kg]
mC = mC/1000;                                     %redukovana hmotnost rotacnich casti hlavni ojnice [Kg]

Fod = (mE + mkloz + mC)*R*(w^2); %odstředivá síla rotačních částí

%% Vyvažování motoru

Fvp = (max(Fsx1)+max(Fsy1))/2;              %síla působící od vývažku [N]

%Vyvažování setrvačných sil posuvných částí prvního řádu
for i = 1:length(a)
Fvpx(i) = Fvp * sin(a(i));                               % x-ová složka síly od vývažku [N]
Fvpy(i) = Fvp * cos(a(i));                              % y-ová složka síly od vývažku [N]
Fnpx(i) = Fsx1(i) - Fvpx(i);                            % zbytková nevývaha v ose x 
Fnpy(i) = Fsy1(i) - Fvpy(i);                           % zbytková nevývaha v ose y
end

%plot(Fnpx,Fnpy)

mv = 1.190;                                                    %hmotnost vývažku
lrv = 41.9;                                                        %vzdálenost těžiště vývažku od osy otáčení KH [mm]
lrv = lrv/1000;                                                %vzdálenost těžiště vývažku od osy otáčení KH [m]

Fv = Fvp+Fod;                                                 %celková síla působící od vývažku     
xp=(Fvp/(2*(w^2)))+(Fod/(2*(w^2)))    %součin hmotnosti vývažku a polohy jeho těžiště požadovaný
xr = mv * lrv                                                    %součin hmotnosti vývažku a polohy jeho těžiště reálný

%% CELKOVÉ SÍLY NA PÍST

for i = 1:length(a)
    Fcm(i) = Fpm(i) - Fsm(i);                          %celková síla působící na hlavní píst [N]
    Fcs(i) = Fps(i) - Fss(i);                               %celková síla působící na první vedlejší píst [N]
    Fct(i) = Fpt(i) - Fst(i);                                %celková síla působící na druhý vedlejší píst [N]
end

%plot(a,Fcm,a,Fcs,a,Fct)

%% SÍLY PŘENÁŠENÉ VEDLEJŠÍMI OJNICEMI + ZRYCHLENÍ DO SOUŘADNÉHO SYSTÉMU

for i = 1:length(a)
    Fos(i) = (Fcs(i)/cos(Bs(i)));
    Fot(i) = (Fct(i)/cos(Bt(i)));
    asx(i) = as(i)*cos(Bs(i));                           % Zrychlení prvního vedlejšího pístu působící v ose ojnice [m/s^2]
    asy(i) = as(i)*sin(Bs(i));                            % Zrychlení prvního vedlejšího pístu působící kolmo k ose ojnice [m/s^2]
    atx(i) = at(i)*cos(Bt(i));                            % Zrychlení druhého vedlejšího pístu působící v ose ojnice [m/s^2]
    aty(i) = at(i)*sin(Bt(i));                             % Zrychlení druhého vedlejšího pístu působící kolmo k ose ojnice [m/s^2]
end

% Zobrazení sil od tlaku plynů s posuvnými setrvačnými silami
subplot(2,4,1)
plot(a,Fps)
subplot(2,4,2)
plot(a,Fss)
subplot(2,4,3)
plot(a,Fcs)
subplot(2,4,4)
plot(a,Fos)
subplot(2,4,5)
plot(a,Fpt)
subplot(2,4,6)
plot(a,Fst)
subplot(2,4,7)
plot(a,Fct)
subplot(2,4,8)
plot(a,Fot)
