clear all; close all; clc;

global alfa L v_f v_max n m r_kol mu_a g yT xTA xTB ksi_val M_prmax T_s1 T_sum

load('alfa.mat')
load('L.mat')
load('v_f.mat')
load('v_max.mat')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input

n=31;                   % pocet useku trati
m=1936+75+461+128.3;    % [kg] LDR operation mass
r_kol=0.323;            % [m] polomer kola LDR
mu_a=0.6;               % soucinitel adheze povrchu pryz-mokra hlina
mu=0.6;                 % soucinitel smykoveho treni pryz-mokra hlina
g=9.81;                 % gravitacni zrychleni
r1=0.0822/2;            % polomer sprocket clusteru vystupni na hrideli prevodovky
r2=0.25913/2;           % polomer sprocketu na hrideli kola
i_pr=12;                % prevodovy pom?r p?evodovky
eta_pr=0.93;            % ú?innost p?evodovky

% Poloha teziste stroje
yT=0.316;           % [m] y souradnice teziste
xTA=0.221;          % [m] x vzdalenost teziste od osy kola A
xTB=0.597;          % [m] x vzdalenost teziste od osy kola B
% A-B = 0.818 m

% Ramena valivych odporu
ksi_val=0.04*r_kol;       % [m] rameno valiveho odporu penumatik
ksi_zl=0.0000;            % [m] rameno valiveho odporu lozisek hrideli kol - zanedbano

% Prevodove pomery
i_rA=r2/r1;             % >1 vystup prevodovky-vstup kolo A
i_rB=r2/r1;             % >1 vystup prevodovky-vstup kolo B
i_r=r2/r1;              % >1

% Ucinnosti
eta_rA=0.945;        % retezovy prevod na kolo A
eta_rB=0.945;        % retezovy prevod na kolo B
eta_r=0.945;

% maximalni vykon mozny prevest na vystupni hridel prevodovky a brzdny
M_pr_uT=(ksi_val+ksi_zl)*0.5*m*g*cos(alfa*pi/180)/(eta_r*i_r);                                      % udržovací moment pohonného mechanismu stroje

Pm_max=        13.3 *ones(n,1); %[kW]                                                               % maximalni vykon 1 elektomotor
Pm_brmax=      13.3 *ones(n,1); %[kW]                                                               % maximalni brzdny vykon 1 elektromotor
P_max_pr=      Pm_max*eta_pr;                                                                       % maximalni pohonny vykon stroje na vystupni hrideli prevodovky
P_brmax_pr=    Pm_brmax/eta_pr;                                                                     % maximalni brzdny vykon stroje na vystupni hrideli prevodovky

% Setrvacne ucinky soustavy
    I0=6.47e-4;     % prevodovka
    I1=1.65e-3;     % SPROCKET CLUSTER 1 - NA HRIDELI PREVODOVKY
    I2=4.548e-2;    % SPORCKET 2 - NA HRIDELI KOL
    IH=1.774e-3;    % HRIDEL KOL
    IK=9.38e-1;     % KOLO STROJE
    mr=3.93;        % HMOTNOST RETEZU

v_s1=10.3; % [km/h] pocatecni rychlost
s_s1=0; % [m] pocatecni draha
T_s1=0; % [s] pocatecni cas
t_inc=0.02; % Casovy inkrement vypoctu

M_prmax=[];
M_pr_u=[];
for i=1:n
    % Normalove sily v kontaktnich bodech kol
    NB=0.5*(m*g*(cos(alfa(i)*pi/180)*xTA-sin(alfa(i)*pi/180)*(r_kol+yT)))/(xTA+xTB);
    NA=0.5*m*g*cos(alfa(i)*pi/180)-NB;
    
    % Tecne sily s predpokladanym prokluzem v kole A
    TA=mu_a*NA;
    TB=0.5*m*g*sin(alfa(i)*pi/180)-TA;
    if TB<=mu_a*NB
        TB=mu*NB;
    end
    
        % Vypocet krouticiho momentu
        M_kolA=TA*r_kol+NA*ksi_val;
        M_kolB=TB*r_kol+NB*ksi_val;
        M_rA=M_kolA+NA*ksi_zl;
        M_rB=M_kolB+NB*ksi_zl;
        M_prmax_A=M_rA/i_rA/eta_rA+M_rB/i_rB/eta_rB;

    % Tecne sily s predpokladanym prokluzem v kole B
    TB=mu_a*NB;
    TA=0.5*m*g*sin(alfa(i)*pi/180)-TB;
    if TA<=mu_a*NA
        TA=mu*NA;
    end
    
        % Vypocet krouticiho momentu
        M_kolB=TB*r_kol+NB*ksi_val;
        M_kolA=TA*r_kol+NA*ksi_val;
        M_rB=M_kolB;
        M_rA=M_kolA;   
        M_prmax_B=M_rA/i_rA/eta_rA+M_rB/i_rB/eta_rB;
    
    % Moment prokluzu je mensi ze dvou vypoctenych momentu
    M_prmax_i=min(M_prmax_A, M_prmax_B);
   
    % Vypocet udrzovaciho momentu
    B=r_kol/i_rA/eta_rA;
    C=r_kol/i_rB/eta_rB;
    D=0.5*m*g*sin(alfa(i)*pi/180);
    E=(ksi_val+ksi_zl)*(NA/i_rA/eta_rA+NB/i_rB/eta_rB);
    M_pr_u_i=C*D+E;
    
    if M_prmax_i<M_pr_u_i;
        disp(i)
        disp('Neni mozne projet usek.')
    end
    
    M_prmax=[ M_prmax; M_prmax_i ];
    M_pr_u=[M_pr_u; M_pr_u_i];
end

M_prmax_celk=min(M_prmax)
M_prmin_celk=max(M_pr_u)
