%Реализация графико-аналитического способа идентификации по переходному
%процессу
%TODO:Заменить на данные с платы
Kzad=0.0794;
Tzad=0.6010;
Ezad=0.4082;
wzad=1;
w = tf([Kzad],[power(Tzad,2) 2*Ezad*Tzad 1]);
%Строим график
ltiview

%по переходному процессу определяем время максимума (tmax), значение в
%максимуме (hmax), установившееся значение (k)
tmax=2.24;
hmax=0.569;
K=0.455/wzad;
%Зная что производная функции w  в максимуме = 0 => находим
E=abs(log(hmax/K-1)/sqrt(power(pi,2)+power(log(hmax/K-1),2)));
T=sqrt(1-power(E,2))*tmax/pi;
%Анализируем точность идентификации
procenterrorofE =(Ezad-E)*100/Ezad;
procenterrorofT =(Tzad-T)*100/Tzad;
procenterrorofK =(Kzad-K)*100/Kzad;
%Из E,T,K найдем Kпр,A,H,Kду,D


%Зададим априорно начальные приближения для алгоритма Качмажа
clear C;
clear d;
clear a;
C(1,1)=0; C(1,2)=0; C(1,3)=0;
U=   [5,164  5,149  5,150  5,148  5,174  5,180  5,180  5,180  5,180  5,180  5,175  5,175  5,175  5,178  5,179  5,180  5,174  5,160  5,154  5,165  5,169  5,160  5,141  5,157  5,162  5,176  5,179  5,180  5,181  5,181  5,178  5,175  5,176  5,179  5,180  5,181  5,181  5,170  5,166  5,153  5,173  5,174  5,168  5,158  5,161  5,166  5,179  5,181  5,180  5,180  5,179  5,177 5,164  5,149  5,150  5,148  5,174  5,180  5,180  5,180  5,180  5,180  5,175  5,175  5,175  5,178  5,179  5,180  5,174  5,160  5,154  5,165  5,169  5,160  5,141  5,157  5,162  5,176  5,179  5,180  5,181  5,181  5,178  5,175  5,176  5,179  5,180  5,181  5,181  5,170  5,166  5,153  5,173  5,174  5,168  5,158  5,161  5,166  5,179  5,181  5,180  5,180  5,179  5,177];
Wzad=[32,582 33,600 34,080 34,188 33,978 34,080 34,080 34,560 33,120 34,080 34,080 34,080 33,600 34,188 34,560 34,231 33,792 34,334 34,080 33,600 34,560 33,600 34,080 34,080 34,080 33,600 34,080 34,188 34,080 33,882 34,080 34,188 33,600 33,978 34,080 33,693 33,513 33,280 33,513 33,109 33,513 33,047 33,600 34,816 33,120 33,431 33,600 33,280 33,120 33,431 33,120 33,513 32,582 33,600 34,080 34,188 33,978 34,080 34,080 34,560 33,120 34,080 34,080 34,080 33,600 34,188 34,560 34,231 33,792 34,334 34,080 33,600 34,560 33,600 34,080 34,080 34,080 33,600 34,080 34,188 34,080 33,882 34,080 34,188 33,600 33,978 34,080 33,693 33,513 33,280 33,513 33,109 33,513 33,047 33,600 34,816 33,120 33,431 33,600 33,280 33,120 33,431 33,120 33,513 ];
d(1,1)=0;
d(1,2)=0;
d(1,3)=0;
for k=2:length(U)-1
    x(k,1)=U(k);
    x(k,2)=U(k-1);
    x(k,3)=Wzad(k-1);
    
    y(k)=U(k+1);
    z(k)=0;
    for j=1:3
    z(k)=z(k)+power(x(k,j),2);   
    end
    
    for j=1:3
    d(k,j)=d(k-1,j)+(y(k)-d(k-1,j)*x(k,j))/z(k);   
        if  (((y(k)-C(k-1,j)*x(k,j))/(z(k)))>(y(k-1)-C(k-1,j)*x(k-1,j))/z(k))
        a(k,j) = 1;
        else
        a(k,j) = 0;
        end;
    C(k,j)=d(k,j)+a(k,j)*(C(k-1,j)-d(k,j));    
    end
end
C
Tint=0.512;
%Задаем H,A
r=30;
diag=35;
h=sqrt(power(diag,2)-power(r,2))
material=7700;
m=material*pi*power(r/1000,2)*h/1000
H=0.5*m*power(r/1000,2)
A=0.25*m*power(r/1000,2)+m*power(h/1000,2)/12
Kpr=(2-C(k,1))*A/T
D=-A*(C(k,2)+C(k,1)-1)/power(T,2)
Kdy=C(k,3)*A/H

K=H*Kdy/Kpr
T=sqrt(A/Kpr)
e=D/(2*sqrt(A*Kpr))