%výpoèet dráhového diagramu
%track diagraph calculation
    clear all; format compact; format shortg;
    
    %zdrojový soubor dráhy
    %source file of track
        filenameT='Data_draha.xlsx';
    %zdrojový soubor vozidla
    %source file of vehicle
        filenameV='Data_vozidlo.xlsx';
        
%název souboru výstupu simulace        
%!!!!zmenit název soouboru aby se nepøepsal vystup minulé simulace!!!!
%!!!!change the filename of otput file to avoid the previous output file
%rewrite!!!!

fileID=fopen('Vystup_simulace_Humpolec_test.txt','w');
        
    %parametry dráhy
    %track parameters
    [ndata,text,alldata]= xlsread(filenameT,1,'D1');
      trackname=cell2mat(text);
       
        %gravitaèní zrychlení
        %gravitation acceleration
        g=xlsread(filenameT,1,'N4');
        
        %adhezní soucinitel
        %coeficient of adhesion
        fi=xlsread(filenameT,1,'N5');
        
        %soucinitel odporu valeni
        %rolling resistance coeficient
        cR=xlsread(filenameT,1,'N6');
        
        %hustova vzduchu
        %air density
        rov=xlsread(filenameT,1,'N7');
        
        %profil dráhy
        
        vzdalenostz=xlsread(filenameT,1,'E:E');
        krivost=xlsread(filenameT,1,'F:F');
        vzdalenostu=xlsread(filenameT,1,'I:I');
        uhel_stoupani=xlsread(filenameT,1,'H:H');
            
        
    %parametry vozidla
    %vehicle parameters
        
        %hmotnost
        %weight
        m=xlsread(filenameV,1,'E30');
        
        %pomerna vzdalenost teziste od zadni napravy
        %relative position of COG by rear axle
        lr=xlsread(filenameV,1,'E32');
        
        %polomìr kola
        %wheel radius
        r=xlsread(filenameV,1,'N11');
       
        %prevodové pomìry
        %gear ratios
        i=[1 2 3 4 5 6 7 8];
        ii=xlsread(filenameV,1,'E21:L21');
        
        %soucinitele rotacnich hmot
        %coeficients of inertia mass
        rod=xlsread(filenameV,1,'E22:L22');
       
        %staly prevod
        %final drive ratio
        i_f=xlsread(filenameV,1,'E19');
        
        %vnejsi charakteristika motoru
        %
        n=xlsread(filenameV,1,'E3:U3');
        Mk=xlsread(filenameV,1,'E4:U4');
        
        %maximální otáèky motoru
        %engine rev limiter
        nmax=xlsread(filenameV,1,'E7');
        %soucinitel odporu vzduchu
        %aerodynamic ceficient
        cx=xlsread(filenameV,1,'E42');
        
        %celni plocha vozidla
        %front surface of vehicle
        Sv=xlsread(filenameV,1,'E40');
       
        %ucinnost pohonneho retezce
        %drive train eficiency
        ny=xlsread(filenameV,1,'E16');
        
        %rychlost vozidla pøi øazení
        %speed of the vehicle while shifting
        A6=xlsread(filenameV,1,'AA13');
        A5=xlsread(filenameV,1,'AB13');
        A4=xlsread(filenameV,1,'AC13');
        A3=xlsread(filenameV,1,'AD13');
        A2=xlsread(filenameV,1,'AE13');
        A1=xlsread(filenameV,1,'AF13');
        A0=xlsread(filenameV,1,'AG13');
        
      %vc1
        for vv=(n(4)/60)*2*pi*r/(ii(1)*i_f):0.0001:500
        
        vc1=vv*3.6;
        nc1=60*vv*ii(1)*i_f/(2*pi*r);
        
        if nc1>nmax
            Fkc1=0;
        else
        Fkc1=(A6*nc1^6+A5*nc1^5+A4*nc1^4+A3*nc1^3+A2*nc1^2+A1*nc1^1+A0*nc1^0)*ii(1)*i_f;
        end
        Fkcc=m*g*fi*lr*r;
        
        deltaFkc1=Fkc1-Fkcc;
        if  deltaFkc1<=0
             break      
        end
        end
        
        
        %v12
        for vv=(n(4)/60)*2*pi*r/(ii(2)*i_f):0.0001:500
        
        v12=vv*3.6;
        n11=60*vv*ii(1)*i_f/(2*pi*r);
        n12=60*vv*ii(2)*i_f/(2*pi*r);
        if n11>nmax
            Fk11=0;
        else
        Fk11=(A6*n11^6+A5*n11^5+A4*n11^4+A3*n11^3+A2*n11^2+A1*n11^1+A0*n11^0)*ii(1)*i_f;
        end
        Fk12=(A6*n12^6+A5*n12^5+A4*n12^4+A3*n12^3+A2*n12^2+A1*n12^1+A0*n12^0)*ii(2)*i_f;
        
        deltaFk11=Fk11-Fk12;
        if  deltaFk11<=0
             break      
        end
        end
        
             %v23
        for vv=(n(4)/60)*2*pi*r/(ii(3)*i_f):0.0001:500
        
        v23=vv*3.6;
        n22=60*vv*ii(2)*i_f/(2*pi*r);
        n23=60*vv*ii(3)*i_f/(2*pi*r);
        if n22>nmax
            Fk22=0;
        else
        Fk22=(A6*n22^6+A5*n22^5+A4*n22^4+A3*n22^3+A2*n22^2+A1*n22^1+A0*n22^0)*ii(2)*i_f;
        end
        Fk23=(A6*n23^6+A5*n23^5+A4*n23^4+A3*n23^3+A2*n23^2+A1*n23^1+A0*n23^0)*ii(3)*i_f;
        
        deltaFk22=Fk22-Fk23;
        if  deltaFk22<=0
             break      
        end
        end
        
             %v34
        for vv=(n(4)/60)*2*pi*r/(ii(4)*i_f):0.0001:500
        
        v34=vv*3.6;
        n33=60*vv*ii(3)*i_f/(2*pi*r);
        n34=60*vv*ii(4)*i_f/(2*pi*r);
        if n33>nmax
            Fk33=0;
        else
        Fk33=(A6*n33^6+A5*n33^5+A4*n33^4+A3*n33^3+A2*n33^2+A1*n33^1+A0*n33^0)*ii(3)*i_f;
        end
        Fk34=(A6*n34^6+A5*n34^5+A4*n34^4+A3*n34^3+A2*n34^2+A1*n34^1+A0*n34^0)*ii(4)*i_f;
        
        deltaFk33=Fk33-Fk34;
        if  deltaFk33<=0
             break      
        end
        end
        
             %v45
        for vv=(n(4)/60)*2*pi*r/(ii(5)*i_f):0.0001:500
        
        v45=vv*3.6;
        n44=60*vv*ii(4)*i_f/(2*pi*r);
        n45=60*vv*ii(5)*i_f/(2*pi*r);
        if n44>nmax
            Fk44=0;
        else
        Fk44=(A6*n44^6+A5*n44^5+A4*n44^4+A3*n44^3+A2*n44^2+A1*n44^1+A0*n44^0)*ii(4)*i_f;
        end
        Fk45=(A6*n45^6+A5*n45^5+A4*n45^4+A3*n45^3+A2*n45^2+A1*n45^1+A0*n45^0)*ii(5)*i_f;
        
        deltaFk44=Fk44-Fk45;
        if  deltaFk44<=0
             break      
        end
        end
        
             %v56
        for vv=(n(4)/60)*2*pi*r/(ii(6)*i_f):0.0001:500
        
        v56=vv*3.6;
        n55=60*vv*ii(5)*i_f/(2*pi*r);
        n56=60*vv*ii(6)*i_f/(2*pi*r);
        if n55>nmax
            Fk55=0;
        else
        Fk55=(A6*n55^6+A5*n55^5+A4*n55^4+A3*n55^3+A2*n55^2+A1*n55^1+A0*n55^0)*ii(5)*i_f;
        end
        Fk56=(A6*n56^6+A5*n56^5+A4*n56^4+A3*n56^3+A2*n56^2+A1*n56^1+A0*n56^0)*ii(6)*i_f;
        
        deltaFk55=Fk55-Fk56;
        if  deltaFk55<=0
             break      
        end
        end
        
             %v67
        for vv=(n(4)/60)*2*pi*r/(ii(7)*i_f):0.0001:500
        
        v67=vv*3.6;
        n66=60*vv*ii(6)*i_f/(2*pi*r);
        n67=60*vv*ii(7)*i_f/(2*pi*r);
        if n66>nmax
            Fk66=0;
        else
        Fk66=(A6*n66^6+A5*n66^5+A4*n66^4+A3*n66^3+A2*n66^2+A1*n66^1+A0*n66^0)*ii(6)*i_f;
        end
        Fk67=(A6*n67^6+A5*n67^5+A4*n67^4+A3*n67^3+A2*n67^2+A1*n67^1+A0*n67^0)*ii(7)*i_f;
        
        deltaFk66=Fk66-Fk67;
        if  deltaFk66<=0
             break      
        end
        end
        
             %v78
        for vv=(n(4)/60)*2*pi*r/(ii(8)*i_f):0.0001:500
        
        v78=vv*3.6;
        n77=60*vv*ii(7)*i_f/(2*pi*r);
        n78=60*vv*ii(8)*i_f/(2*pi*r);
        if n77>nmax
            Fk77=0;
        else
        Fk77=(A6*n77^6+A5*n77^5+A4*n77^4+A3*n77^3+A2*n77^2+A1*n77^1+A0*n77^0)*ii(7)*i_f;
        end
        Fk78=(A6*n78^6+A5*n78^5+A4*n78^4+A3*n78^3+A2*n78^2+A1*n78^1+A0*n78^0)*ii(8)*i_f;
        
        deltaFk77=Fk77-Fk78;
        if  deltaFk77<=0
             break      
        end
        end
        
        %èas pøeøazení
        %shift time
        tsh=xlsread(filenameV,1,'E24');
        
        %maximální skluz
        %maximal slip
        slpmax=xlsread(filenameV,1,'N10');
        
        %pohaena naprava
        %driving axle
        dax=xlsread(filenameV,1,'I15');
        
       if dax<=1
        xx=lr;
       elseif dax<=2
           xx=lr-1;
       else
           xx=1;
       end
       
       v00=sqrt(fi*g/krivost(1,:));
       
       if v00>140
           v00=0;
       else
           v00=v00;
       end
 
sim('drahovy_diagram_sim');
       
        
t=vystup_simulace(:,1);
s=vystup_simulace(:,2);
v=vystup_simulace(:,3);
a=vystup_simulace(:,4);
vlim=vystup_simulace(:,5);
sg=vystup_simulace(:,6);
engrevs=vystup_simulace(:,7);
slip=vystup_simulace(:,8);

sizer=size(engrevs);
engrevsl=zeros;
for blp=5:sizer(1)-5
engrevsl(blp-4)=engrevs(blp);
end

tlap=max(t);
vmax=max(v);
vmin=min(v);
slap=max(s);
vavg=(slap/tlap)*3.6;
amax=max(a);
amin=-min(a);
amaxg=amax/g;
aming=amin/g;
maxengrevs=max(engrevsl);
minengrevs=min(engrevsl);

sizek=size(krivost);
R=[];
for ksi=1:sizek(1)
R(ksi)=1/krivost(ksi);
end

disp('název dráhy / track name');
disp(trackname);
disp('èas kola / lap time =');
disp(tlap);
disp('[s]');
disp('maximální rychlost / top speed =');
disp(vmax);
disp('[km/h]');


figure (1)
plot(t,v,'b',t,vlim,':k',t,sg,'m');
title ('dráhový diagram / track diagraph');
xlabel('èas/time [s]');
xlim([0 t(end)]);
ylabel('rychlost/speed [km/h]     zaøazený stupeò/selected gear [1/10]');
ylim([0 vmax+15]);


figure (2)
plot(s,v,'b',s,vlim,':k',s,sg,'m');
title ('dráhový diagram / track diagraph');
xlabel('dráha/track [m]');
xlim([0 s(end)]);
ylabel('rychlost/speed [km/h]     zaøazený stupeò/selected gear [10]');
ylim([0 vmax+15]);


figure (3)
plot(s,a,'b');
title ('dráhový diagram / track diagraph');
xlabel('dráha/track [m]');
xlim([0 s(end)]);
ylabel('zrychlení/acceleration [m/s^2]');
ylim([min(a)-2 max(a)+2]);

figure (4)
plot(s,engrevs,'b');
title ('dráhový diagram / track diagraph');
xlabel('dráha/track [m]');
xlim([0 s(end)]);
ylabel('otáèky motoru/engine revs [1/min]');
ylim([0 maxengrevs+15]);

figure (5)
plot(s,slip,'b');
title ('dráhový diagram / track diagraph');
xlabel('dráha/track [m]');
xlim([0 s(end)]);
ylabel('skluz/slip [%]');
ylim([0 100]);

fprintf(fileID,'%6s\r\n','Vystup simulace / Simulation output');
fprintf(fileID,'%2s\r\n','--------------------------------------------------------------');
fprintf(fileID,'%6s\r\n','vstupní hodnoty /input values');
fprintf(fileID,'%2s\r\n','--------------------------------------------------------------');
fprintf(fileID,'%6s\r\n','data dráhy / track data');
fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%20s %45s\r\n','nazev drahy / track name',trackname);
fprintf(fileID,'%20s %20.2f %10s\r\n','gravitaèní zrychlení / gravitation acceleration',g,'[m/s^2]');
fprintf(fileID,'%20s %27.2f %10s\r\n','adhezní soucinitel / adhesion coeficient',fi,'[1]');
fprintf(fileID,'%20s %11.2f %10s\r\n','soucinitel odporu valeni / rolling resistance coeficient',cR,'[1]');
fprintf(fileID,'%20s %38.2f %10s\r\n','hustova vzduchu / air density',rov,'[kg/m^3]');
fprintf(fileID,'%20s\r\n','profil dráhy / track proile');
fprintf(fileID,'%6s','s       [m]');fprintf(fileID,'%12.2f',vzdalenostz);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s','R       [m]');fprintf(fileID,'%12.2f',R);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s','s       [m]');fprintf(fileID,'%12.2f',vzdalenostu);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s','alpha   [°]');fprintf(fileID,'%12.2f',uhel_stoupani);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s\r\n','--------------------------------------------------------------');
fprintf(fileID,'%6s\r\n','data vozidla / vehicle data');
fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%10s %74.2f %10s\r\n','hmotnost / weight',m,'[kg]');
fprintf(fileID,'%10s %8.2f %10s\r\n','pomìrná vzdalenost tìziste od zadní nápravy / relative position of COG by rear axle',lr,'[1]');
fprintf(fileID,'%10s %64.2f %10s\r\n','polomìr kola / wheel radius',r,'[m]');
fprintf(fileID,'%10s %41.2f %10s\r\n','soucinitel odporu vzduchu / aerodynamic coeficient',cx,'[1]');
fprintf(fileID,'%10s %47.2f %10s\r\n','celni plocha vozidla / front area of vehicle',Sv,'[m^2]');
fprintf(fileID,'%10s %41.2f %10s\r\n','uèinnost pohonneho øetìzce / drive train eficiency',ny,'[1]');
fprintf(fileID,'%10s %60.2f %10s\r\n','rychlost pøeøazení / shift time',tsh,'[s]');
fprintf(fileID,'%10s %61.0f %10s\r\n','maximální skluz / maximal slip',slpmax,'[%]');
if dax<2
   fprintf(fileID,'%10s %60s\r\n','pohánìná náprava / driving axle','FWD');
elseif dax<3
   fprintf(fileID,'%10s %60s\r\n','pohánìná náprava / driving axle','RWD');
else
    fprintf(fileID,'%10s %60s\r\n','pohánìná náprava / driving axle','AWD');
end
fprintf(fileID,'%10s %59.2f %10s\r\n','stalý pøevod / final drive ratio',i_f,'[1]');
fprintf(fileID,'%10s\r\n','prevodové pomìry / gear ratios                              [1]');
fprintf(fileID,'%12.0f',i);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%12.3f',ii);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%10s\r\n','soucinitele rotacnich hmot / coeficients of inertia mass     [1]');
fprintf(fileID,'%12.0f',i);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%12.2f',rod);fprintf(fileID,'%6s\r\n',' ');

fprintf(fileID,'%10s\r\n','vnejsi charakteristika motoru / outer engine characteristic');
fprintf(fileID,'%6s','n   [1/min]');
fprintf(fileID,'%12.0f',n);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s','Mk     [Nm]');
fprintf(fileID,'%12.2f',Mk);fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s\r\n','--------------------------------------------------------------');

fprintf(fileID,'%6s\r\n','vysledky simulace / simulation results');
fprintf(fileID,'%6s\r\n',' ');
fprintf(fileID,'%6s %37.3f %8s\r\n','èas kola / lap time',tlap,'[s]');
fprintf(fileID,'%6s %24.2f %9s\r\n','delka jednoho kola / lap lenght',slap,'[m]');
fprintf(fileID,'%6s %24.1f %10s\r\n','maximální rychlost / top speed',vmax,'[km/h]');
fprintf(fileID,'%6s %20.1f %10s\r\n','minimální rychlost / minimum speed',vmin,'[km/h]');
fprintf(fileID,'%6s %21.1f %10s\r\n','prùmìrná rychlost / average speed',vavg,'[km/h]');
fprintf(fileID,'%6s %7.0f %12s\r\n','maximální otáèky motoru / maximum engine revs',maxengrevs,'[1/min]');
fprintf(fileID,'%6s %7.0f %12s\r\n','minimální otáèky motoru / minimum engine revs',minengrevs,'[1/min]');
fprintf(fileID,'%6s %13.2f %9s\r\n','maximální zrychlení / maximum acceleration',amax,'[m/s^2]');
fprintf(fileID,'%6s %13.2f %9s\r\n','                                          ',amaxg,'[*G]');
fprintf(fileID,'%6s %12.2f %9s\r\n','maximální  zpomalení / maximum deceleration',amin,'[m/s^2]');
fprintf(fileID,'%6s %12.2f %9s\r\n','                                           ',aming,'[*G]');