%% 2. järku mudeli identifitseerimine

syms T_1 T_2 Ko u % Mudel sümbolkujul
syms y(t)
Mudel=(T_1*T_2)*diff(diff(y(t),t))+(T_1+T_2)*diff(y(t),t)+y(t)==Ko*u

% Hüppelise sisendi +500 W KATSEANDMED: ajavektor tt [h] ja nendele ajahetkedele vastavad väljundi väärtused y1(t)

tt=[0 0.5 1   2    3    4    5    6  10]
y1=[8 9.0 11 13.6 14.8 15.2 15.4 15.45 15.5]
plot(tt,y1)% katseandmete graafiku väljatrükk

%
% 1.Staatilise ülekandeteguri Ko=\frac{y\left(\infty \right)-y\left(0\right)}{u\left(\infty \right)-u\left(0\right)} määramine:

deltaU=500; % sisendi muutus
deltaY=15.5 - 8 % väljundi muutus
K0=deltaY/deltaU

%
% 2. Väljundi karakteristiku  normeerimine:
%

yy=y1-8;
yy=yy/7.5; % normeeritud karakteristik yy
plot(tt,yy), grid on %väljatrükk

%
% 3. Abimuutujate leidmine (G.Ormanns'i meetod):

t7=spline(yy,tt,0.7) % normeerit. y väärtusel 0.7
t3=t7/3
y3=spline(tt,yy,t3)

%
% 4. Vastavalt y3-le leiame tabelist \chi^2

y3n=[0.189 0.20 0.21 0.22 0.24 0.260 0.28 0.30 0.32];
hi2=[0.000 0.11 0.20 0.27 0.41 0.525 0.63 0.75 0.87];
hii_ruut=spline(y3n,hi2,y3); 
plot(hi2,y3n), hold on, plot(hii_ruut,y3,'o'), xlabel('χ 2'), ylabel('y3'), grid on
hii_ruut, hii=sqrt(hii_ruut), hold off

%
% 5. Ajakonstantide T1 T2 leidmine

T= t7/2.4; T1= T*(1+hii), T2= T*(1-hii)

%
% 6. Mudeli dif võrrandi lahend sümbolkujul ja parameetritega

Output=dsolve(Mudel,y(0)==8) % poolikute algtingimustega
Mudel_p = subs(Mudel,[T_1,T_2,Ko,u],[T1, T2, K0, 500+533]);
Dy=diff(y,t); %algtingimuse esitamiseks
valjund=dsolve(Mudel_p,[y(0)==8, Dy(0)==0])

%
% 7. Võrdleme mudeli väljundit ja katseandmeid samal graafikul

fplot(valjund,[0 10]), hold on, plot(tt,y1), hold off