function pirani_PID_odpor(Rcil,r,sumy,tmax)
%pirani_PID_odpor(Rcil,r,sumy,tmax)
%Simuluje casovy prubeh napeti, proudu, odporu a teploty jednoduche Piraniho merky 
% pri pouziti PID regulace. 
%"Rcil" je cilova hodnota odporu Piraniho merky (pocatecni hodnota je nastavena na 30 st.), 
%"r" je triprvkovy vektor parametru PID regulatoru (zesileni regulatoru, jednotky [V/Ohm V/Ohm*s s*V/Ohm]), 
%"sumy" je parametr ridici velikost pripadneho sumu, ktery se pricita k merene odchylce odporu od cilove hodnoty, 
%"tmax" je casovy interval [s], po ktery se chovani Piraniho merky bude pocitat. 

U0=0; %pocatecni napeti
T0=30; %pocatecni teplota vlakna [stupne Celsia]
mezeU=[0 30]; %meze, ve kterem napeti musi byt, tj. napeti nelze snizit pod mezeU(1) ani zvysit nad mezeU(2)
R0=24; %Odpor pri 0 stC 
alfaR=0.13; %teplotni soucinitel odporu
C=0.5; %Tepelna kapacita vlakna

%Koeficienty odvodu tepla od vlakna, vztazeny k abs. teplote [K]:
Kl=5e-3; %Koeficient ztraty tepla u linearniho clenu vuci T
Kz=5e-12; %Koeficient ztraty tepla u clenu s T^4
K0=-(Kl*300 + Kz*300^4); %Koeficient ztraty tepla u konstantniho clenu, pocita se s tim, ze vlakno je v rovnovaze s okolim pri teplote 300 K. 

dt=0.01; %Casovy krok vypoctu

if length(r)==2
 r=[r 0]; 
elseif length(r)==1
 r=[r 0 0]; 
end

takt=0; Takt=T0; Uakt=U0; 
rint=0; %Pocatecni integral odchylky odporu
R=[]; T=[]; U=[]; I=[]; t=[];

%Vlastni vypocet prubehu regulace: 
while takt<tmax, 
 Rakt=R0+alfaR*Takt; %Odpor pri aktualni teplote [stupne C]
 Iakt=Uakt/Rakt; %Aktualni proud
 dT=Uakt*Iakt - (K0 + Kl*(273.15+Takt) + Kz*(273.15+Takt)^4); %Tepelna bilance vlakna
 dT=dT/C; %Podeleni dodane energie tepelnou kapacitou
 dT=dT*dt; %Zmena teploty za cas dt
 R=[R; Rakt]; U=[U; Uakt]; I=[I; Iakt]; T=[T; Takt]; t=[t; takt];
 Takt=Takt+dT; 
 if length(r) ~=0
  if takt~=0
   odchmin=odch;
  end
  odch=Rcil-Rakt; %Odchylka ciloveho a aktualniho odporu
  if length(sumy)~=0
   odch=odch+sumy*randn; %Pricteni sumu ke ctene odchylce odporu 
  end
  rint=rint+odch*dt; %Integral odchylky
  if takt==0
   rder=0;
  else
   rder=(odch-odchmin)/dt; %Casova derivace odchylky
  end
  Uakt = U0 + r(1)*odch + r(2)*rint + r(3)*rder; %Regulace napeti
  %Omezeni mezi napeti: 
  if Uakt < mezeU(1)
   Uakt=mezeU(1);
  elseif Uakt > mezeU(2) 
   Uakt=mezeU(2);
  end
 end
 takt=takt+dt; 
end

%Grafeni: 
figure; plot(t,U,'b'); xlabel('Cas [s]'); ylabel('Napeti'); grid on; 
pos=get(gcf,'Position'); 
pos(1)=pos(1)-310; 
set(gcf,'Position',pos);

figure; plot(t,I,'k'); xlabel('Cas [s]'); ylabel('Proud'); grid on;
pos=get(gcf,'Position'); 
pos(1)=pos(1)+310; 
set(gcf,'Position',pos);

figure; plot(t,R,'g'); xlabel('Cas [s]'); ylabel('Odpor'); grid on;
pos=get(gcf,'Position'); 
pos(2)=pos(2)-500; pos(1)=pos(1)-310; 
set(gcf,'Position',pos);

figure; plot(t,T,'r'); xlabel('Cas [s]'); ylabel('Teplota'); grid on;
pos=get(gcf,'Position'); 
pos(2)=pos(2)-500; pos(1)=pos(1)+310; 
set(gcf,'Position',pos);
