% function [A0,A1,A2,A3,A4,A5,A6,A7] = tc2areas (Tpv,Tzv,Tdelay,Kpr);
%
% Function zp2areas calculates characteristic areas (A0 to A7) from the process 
% time constants:
%
% Gpr = Kpr*(1+s*Tzv(1))*(1+s*Tzv(2))*...*exp(-Tdelay*s)/((1+s*Tpv(1))*(1+s*Tpv(2))*...)

function [A0,A1,A2,A3,A4,A5,A6,A7] = tc2areas (Tpv,Tzv,Tdelay,Kpr);

Len_num=length(Tzv);
Len_den=length(Tpv);

num = 1;
for k = 1:Len_num
   num = conv(num,[Tzv(k) 1]);
end;
	   
den = 1;
for k = 1:Len_den
   den = conv(den,[Tpv(k) 1]);
end;

num = real(num);		% Just in case if complex numerator or denominator are calculated
den = real(den);

Len_num=length(num);
Len_den=length(den);

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

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

Len_num=length(num);
Len_den=length(den);

a1 = den(Len_den-1);
a2 = den(Len_den-2);
a3 = den(Len_den-3);
a4 = den(Len_den-4);
a5 = den(Len_den-5);
a6 = den(Len_den-6);
a7 = den(Len_den-7);
a8 = den(Len_den-8);

b1 = num(Len_num-1);
b2 = num(Len_num-2);
b3 = num(Len_num-3);
b4 = num(Len_num-4);
b5 = num(Len_num-5);
b6 = num(Len_num-6);
b7 = num(Len_num-7);
b8 = num(Len_num-8);

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;