% M-plik part2;
% pkt.4 
% Identyfikacja modelu na podstawie cech odpowiedzi skokowej
% przy charakterze oscylacyjnym odpowiedzi
% autor: Marcin Maślanka

load data u y t Fs
C= 0.6*10^(-6); % zadana pojemność w [F]

clc; 
disp(['Wyniki obliczeń dla zadanej odpowiedzi skokowej:']); 
disp(' ');
A1= max(y)-1; % wartość maksymalna -1 --> A1
disp(['* Amplituda A1= ',num2str(A1)]);
ind1= find(y==(A1+1)); % numer próbki oodpowiadający A1
disp(['* Numer próbki odpowiadający amplitudzie A1: ',num2str(ind1)]);
t1= t(ind1); % czas t1 oodpowiadający amplitudzie A1
disp(['* Czas t1 odpowiadający amplitudzie A1: t1= ',num2str(t1),'[s]']);

bufor= y(ind1:length(y)); % chwilowe obcięcie odpowiedzi w celu wyznaczenia A2
A2= 1- min(bufor); % 1- wartość minimalna z tak obciętego sygnału --> A1
disp(['* Amplituda A2= ',num2str(A2)]);
ind2= find(y==(1-A2)); % numer próbki oodpowiadający A2
disp(['* Numer próbki odpowiadający amplitudzie A2: ',num2str(ind2)]);
t2= t(ind2);
disp(['* Czas t2 odpowiadający amplitudzie A2: t2= ',num2str(t2),'[s]']);
clear ind1 ind2 bufor 
disp(' ');

% Okres drgań tłumionych:
T= 2*(t2-t1);
disp(['* Okres drgań tłumionych T= ',num2str(T),'[s]']);
% Pulsacja drgań tłumionych:  
omega= 2*pi/T;
disp(['* Pulsacja drgań tłumionych omega= ',num2str(omega),'[rad/s]']);
% Tłumienność:
delta= log(A1/A2)*omega/pi;
disp(['* Tłumienność delta= ',num2str(delta),'[rad/s]']);
% Pulsacja drgań własnych:
omega0= sqrt(omega^2+delta^2);
disp(['* Pulsacja drgań własnych omega0= ',num2str(omega0),'[rad/s]']);
% Okres drgań nietłumionych:
T0= 1/omega0;
disp(['* Okres drgań nietłumionych T0= ',num2str(T0),'[s]']);
% Współczynnik tłumienia:
dzeta= delta/omega0;
disp(['* Współczynnik tłumienia dzeta= ',num2str(dzeta)]);
disp(' ');

% Poszukiwane parametry R,L:
disp(['Poszukiwane parametry wynoszą:']); 
R= 2*T0*dzeta/C;
disp(['* Rezystancja R= ',num2str(R),'[Om]']);
L= T0^2/C;
disp(['* Impedancja L= ',num2str(L),'[Henr]']);
disp(' ');

% Transmitancja operatorowa modelu:
disp('Otrzymana transmitancja operatorowa modelu:');
num= [1];	den= [T0^2 2*T0*dzeta 1];
gc= tf(num,den)

yh = lsim(gc,u,t); % symulowana odpowiedź modelu na wymuszenie u

fit = norm(yh - y)/sqrt(length(y)); % obliczenie błędu średniokwadratowego

a= figure; 
set(a,'Numbertitle','off','name', ...
	'Wynik weryfikacji otrzymanego modelu funkcją compare');

% Odpowiedź układu o wyznaczonej transmitancji na wymuszenie u:
a= figure;
set(a,'Numbertitle','off','name', ...
'Różnica między odpowiedzią otrzymanego modelu i badanego układu na skok jedn.');

plot((y-y1)); xlabel('numer próbki'); ylabel('amplituda błędu weryfikacji'); grid on;