% function [Ki,K,Kd] = area2PID (A0,A1,A2,A3,A4,A5,Kmax);
%
% Function araa2PID calculates parameters of the PID controller:
%
% u = (Ki/s + K + Kd*s) * e
%
% from the measured areas of the process A0 to A5 (A0 is the process steady-state gain)
% The parameter Kmax represents the highest allowed open-loop gain K*A0 


function [Ki,K,Kd] = area2PID (A0,A1,A2,A3,A4,A5,Kmax);

Num = A3*A3 - A1*A5;												% Numerator
Den = 2*(A1*A2*A3+A0*A1*A5-A1*A1*A4-A0*A3*A3);			% Denominator

if (Num == 0)			
   K = 0;
elseif (Den == 0)			% Division with zero
   if (A0 ~= 0)
      K = Kmax/A0;
   else
      K = Kmax;
   end;
else							% Everything is OK
   K = Num/Den;
end;

Tmp = K*A0;							% Nominal gain

if (Tmp > Kmax) | (Tmp < 0)	% If the calculated gain is too high or negative	
   K = Kmax/A0;
end;

if (A1 ~= 0)
   Kd = (2*K*(A1*A2-A0*A3)-A3)/(2*A1*A1);
else
   Kd = 0;
end;

Tmp = Kd*A0;							% Nominal gain of D-term

if (Tmp < 0)							% If derivative time constant is negative
   Kd = 0;
end;

if (Kd == 0)							% If negative, calculate PI controller parameters
   [Ki,K] = area2PI (A0,A1,A2,A3,Kmax);
else
   Ki = (2*K*A0 + 1)/(2*A1);
end;