% function [A0,A1,A2,A3,A4,A5,A6,A7] = zp2areas (z,p,K,Tdelay);
%
% Function zp2areas calculates characteristic areas (A0 to A7) from the process TF parameters:
%
% Gpr = K*(s+z(1))*(s+z(2))*...*exp(-Tdelay*s)/((s+p(1))*(s+p(2))*...)

function [A0,A1,A2,A3,A4,A5,A6,A7] = zp2areas (z,p,K,Tdelay);

[num,den] = zp2tf(z,p,K);

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

if (b<9),
   den=[zeros(size(1:(9-b))) den];
   b = 9;
end

if (a<b),
   num=[zeros(size(1:(b-a))) num];
end

a0 = den(b);
den = den/a0;
a1 = den(b-1);
a2 = den(b-2);
a3 = den(b-3);
a4 = den(b-4);
a5 = den(b-5);
a6 = den(b-6);
a7 = den(b-7);

b0 = num(b);
num = num/b0;
b1 = num(b-1);
b2 = num(b-2);
b3 = num(b-3);
b4 = num(b-4);
b5 = num(b-5);
b6 = num(b-6);
b7 = num(b-7);

Kpr = b0/a0;

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

A0 = Kpr;
A1 = A0*(a1-b1+Tdelay);
A2 = A0*(b2-a2-Tdelay*b1+Td2/2)+A1*a1;
A3 = A0*(a3-b3+Tdelay*b2-Td2*b1/2+Td3/6)+A2*a1-A1*a2;
A4 = A0*(b4-a4-Tdelay*b3+Td2*b2/2-Td3*b1/6+Td4/24)+A3*a1-A2*a2+A1*a3;
A5 = A0*(a5-b5+Tdelay*b4-Td2*b3/2+Td3*b2/6-Td4*b1/24+Td5/120)+A4*a1-A3*a2+A2*a3-A1*a4;
A6 = A0*(b6-a6-Tdelay*b5+Td2*b4/2-Td3*b3/6+Td4*b2/24-Td5*b1/120+Td6/720)+A5*a1-A4*a2+A3*a3-A2*a4+A1*a5;
A7 = A0*(a7-b7+Tdelay*b6-Td2*b5/2+Td3*b4/6-Td4*b3/24+Td5*b2/120-Td6*b1/720+Td7/5040)+A6*a1-A5*a2+A4*a3-A3*a4+A2*a5-A1*a6;