% function [A0,A1,A2,A3,A4,A5] = numden2areas (num,den,Tdelay);
%
% Tp1..Tp8 process poles (time constants) 
% Tz1..Tz8 process zeros (time constants)
% Tdelay ... time delay
% Kpr ... process steady-state gain
% Kpmax ... max. proportional gain K*Kpr
%
% Kp ... controller gain, Ti .. int. time constant, Td ... dif. time const.
% alphad ... parameter
% A1 .. A5 areas
% a1 .. a6, b1 .. b6 calculated process parameters


function [A0,A1,A2,A3,A4,A5] = numden2areas (num,den,Tdelay);

a=max(size(den));
b=max(size(num));

den0 = den(a);
num0 = num(b);
Kpr = num0/den0;

for k = 1:(a-1),
   av(k) = den(a-k);
end;
for k = a:8,
   av(k) = 0;
end;

for k = 1:(b-1),
   bv(k) = num(b-k);
end;
for k = b:8,
   bv(k) = 0;
end;

av=av/den0;
bv=bv/num0;

Td2 = Tdelay*Tdelay;
Td3 = Td2*Tdelay;
Td4 = Td3*Tdelay;
Td5 = Td4*Tdelay;

A0 = 1;
A1 = av(1)-bv(1)+Tdelay;
A2 = bv(2)-av(2)+A1*av(1)-Tdelay*bv(1)+Td2/2;
A3 = av(3)-bv(3)+A2*av(1)-A1*av(2)+Tdelay*bv(2)-Td2*bv(1)/2+Td3/6;
A4 = bv(4)-av(4)+A3*av(1)-A2*av(2)+A1*av(3)-Tdelay*bv(3)+Td2*bv(2)/2-Td3*bv(1)/6+Td4/24;
A5 = av(5)-bv(5)+A4*av(1)-A3*av(2)+A2*av(3)-A1*av(4)+Tdelay*bv(4)-Td2*bv(3)/2+Td3*bv(2)/6-Td4*bv(1)/24+Td5/120;

A0=Kpr;
A1=A1*Kpr;
A2=A2*Kpr;
A3=A3*Kpr;
A4=A4*Kpr;
A5=A5*Kpr;