dt=1e-5;%[s]
tkonec=2.5;
ikonec=floor(tkonec/dt);

%parametry regulatoru
kr=0.01;

Taur=10;
%Taur=1/(10^0.75);
%Taur=0.13;
%Taur=0.027;
%Taur=0.3;

fp=5000;%frekvence nosne pily, Hz
urmax=1;%rozsah ridiciho napeti mustku
U=10;%napajeci napeti [V]
R=1/0.3;%odpor vinuti [Ohm]
L=1e-2*R;%indukcnost vinuti [H]
TauI=1e-3;

%pocatecni podminky
t=0;i=0;sum=0;ur=0;otw=0;ot=0;
u=0;Ui=0;Mz=0;kfi=0.1;
vys=zeros(7,ikonec);

%vlastni vypocetni smycka
for icas=1:ikonec
    t=(icas-1)*dt;
    %definovani pozadovaneho prubehu otacek
    if t>0.02 otw=100; end;
    %definovani prubehu poruchy
%if t>1.6 Mz=0.1; end;
if t>1 Mz=0.1; end;
%if t>0.8 Ui=1; end;
%Ui=kfi*ot;
    
    e=otw-ot; %regulacni odchylka

    %PI (P) regulator bez omezovace
    sum=sum+1/Taur*e*dt; %integrace reg. odchylky
    ur=kr*(e+sum);
%ur=ur+Ui/U*urmax;
    %PI (P) regulator s omezovacem
    %if abs(ur)<urmax sum=sum+1/Taur*e*dt; end;%integrace reg. odchylky, pokud neni vystup omezen
    %ur=kr*(e+sum);
    %if ur>urmax ur=urmax; end;
    %if ur<-urmax ur=-urmax; end;
        
    % vypocet diference proudu (z rovnice u=R*i+L*di/dt)
    di_dt=1/L*(u-R*i-Ui);
    % vypocet diference napeti (z prenosu u/ur=U/urmax/(1+p/(2*fp)) )
    du_dt=2*fp*(U/urmax*ur-u);
    dot=i-Mz/kfi;
    % numericka integrace Eulerovou metodou
    i=i+di_dt*dt;
    u=u+du_dt*dt;
    ot=ot+1/TauI*dot*dt;
    % zapis vysledku do pole
    vys(:,icas)=[t,otw,ur,u,Ui,ot,Mz]';
end;


if isempty(findobj('type','figure','name','RT03 prechodove'))
    fh=figure('name','RT03 prechodove');
    set(gcf,'color',[1,1,1]);
    pismo=13;

    MzGraf=subplot(3,1,3);
    set(gca,'nextplot','add')
    set(gca,'outerposition',[0,0.005,1.0000,0.15])
    set(gca,'FontName','Helvetica','FontSize',pismo);

    UiGraf=subplot(3,1,2);
    set(gca,'nextplot','add')
    set(gca,'outerposition',[0,0.15,1.0000,0.15])
    set(gca,'FontName','Helvetica','FontSize',pismo);

    Prechodova=subplot(3,1,1);
    set(gca,'nextplot','add')    
    set(gca,'outerposition',[0,0.3,1.0000,0.65])
    set(gca,'FontName','Helvetica','FontSize',pismo);
    xlabel('time [s]');
end

set(fh, 'currentaxes', Prechodova);
plot(vys(1,:),vys(2,:),'k','LineWidth',2);
plot(vys(1,:),vys(6,:),'LineWidth',2);
legend('w','y');

set(fh, 'currentaxes', UiGraf);
plot(vys(1,:),vys(5,:),'LineWidth',2);
legend('Ui');

set(fh, 'currentaxes', MzGraf);
plot(vys(1,:),vys(7,:),'LineWidth',2);
legend('Mz');