function [a,P]=aP(Q,n)

% aP.m Function  [a,P]=aP(Q,n)
% a and P are process identification parameters for
% system with two (n=2) or three (n=3) time constants.
% Q = t_0.66 / t_0.33  from S-shaped step response.
%
%   100% |---------____  
%        |        /
%    66% |-------/
%        |      /|
%    33% |-----/ |    Q = t_0.66 / t_0.33
%        |    /| |
%     0% |___/ | |
%        |-----|-|-----   
%         <---> t_0.33
%         <------> t_0.66
%
% Method is described in "Bertil Thomas - Modern reglerteknik"
% LIBER 1992 ISBN 91-47-05085-3
% Table values are taken from fig 7.1 and 7.2
%
% T is the leading time constant. Usage of a and P:
%  n = 2   1.92 < Q < 2.71  T = t_0.66/P(1+a)
%                 K
%    G(s) = -------------  Two time constants
%           (1+Ts)(1+aTs)
%
%  n = 3   1.68 < Q < 2.71  T = t_0.66/P(1+a+a^2)
%                    K
%    G(s) = ----------------------  Three time constants
%           (1+Ts)(1+aTs)(1+a^2Ts)
%

warning('off', 'all') % don't display warnings about NaN

q=1.5:0.02:2.7;

a3=[ NaN  NaN  NaN  NaN NaN  ...
NaN  NaN  NaN  NaN  1.00 ...
0.83 0.72 0.63 0.57 0.52 ...
0.485 0.45 0.425 0.40 0.37 ...
0.35 0.335 0.315 0.30 0.28 ...
0.265 0.25 0.235 0.225 0.21 ...
0.20 0.19 0.17 0.16 0.15 ...
0.145 0.14 0.13 0.12 0.11 ...
0.105 0.10 0.09 0.085 0.08 ...
0.0725 0.065 0.06 0.055 0.05 ...
0.045 0.04 0.035 0.03 0.025 ...
0.02 0.015 0.0125 0.01 0.005  0 ];

a2= [ NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN   1   0.8 0.55 0.475 ...
0.4125 0.375 0.34 0.31 0.29 ...
0.27 0.25 0.23 0.215 0.20 ...
0.18 0.165 0.155 0.14 0.13 ...
0.11875 0.1075 0.10 0.09 0.08 ...
0.0725 0.065 0.06 0.055 0.05 ...
0.045 0.04 0.035 0.03 0.025 ...
0.02 0.015 0.015 0.01 0.005  0 ];

P3=[ NaN  NaN  NaN  NaN  NaN  ...
NaN  NaN  NaN  NaN  1.145 ...
1.14 1.135 1.13075 1.1265 1.12325 ...
1.12 1.118 1.116 1.1135 1.11175 ...
1.11 1.1085 1.1075 1.106 1.10475 ...
1.1035 1.1025 1.1015 1.1005 1.0995 ...
1.0985 1.0975  1.097 1.0965  1.096 ...
1.0955 1.095 1.0945 1.0942 1.094 ...
1.0935 1.0935 1.0935 1.0935 1.0935 ...
1.0935 1.0937  1.0937 1.0937 1.094 ...
1.094 1.094 1.0945 1.0947 1.095 ...
1.095 1.0955 1.0955 1.096 1.0962 1.0965 ];

P2=[ NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN  NaN  NaN  NaN  NaN ...
NaN   1.145   1.13975 1.1345 1.127 ...
1.1225 1.1175 1.1135 1.11 1.1075 ...
1.105 1.1035 1.102 1.10 1.099 ...
1.098 1.097 1.0965 1.096 1.0952 ...
1.095 1.0945 1.094 1.0937 1.0935 ...
1.0935 1.0935 1.0935 1.0937 1.094 ...
1.094 1.094 1.0945 1.0947 1.095 ...
1.095 1.0955 1.0955 1.096 1.0962 1.0965 ];

if (n==2) && (Q>1.92) && (Q<=2.71)
    a=interp1(q,a2,Q);
    P=interp1(q,P2,Q); 
    subplot(1,2,1)
    plot(q,a2,q,a3,Q,a,'o')
    title('a(Q), n=2 success')
    xlabel('Q')
    grid on
    subplot(1,2,2)
    plot(q,P2,q,P3,Q,P,'o')
    title('P(Q), n=2 success')
    xlabel('Q')
    grid on
elseif (n==2)
    disp('Observe! n = 2  1.92 < Q < 2.71')
    a=NaN;
    P=NaN;
    subplot(1,2,1)
    plot(q,a2,q,a3)
    title('a(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
    subplot(1,2,2)
    plot(q,P2,q,P3)
    title('P(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
elseif (n==3) && (Q > 1.68) && (Q < 2.71)
    a=interp1(q,a3,Q);
    P=interp1(q,P3,Q);
    subplot(1,2,1)
    plot(q,a2,q,a3,Q,a,'o')
    title('a(Q), n=3 success')
    xlabel('Q')
    grid on
    subplot(1,2,2)
    plot(q,P2,q,P3,Q,P,'o')
    title('P(Q), n=3 success')
    xlabel('Q')
    grid on
elseif (n==3)
    disp('Observe! n = 3  1.68 < Q < 2.71')
    a=NaN;
    P=NaN;
    subplot(1,2,1)
    plot(q,a2,q,a3)
    title('a(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
    subplot(1,2,2)
    plot(q,P2,q,P3)
    title('P(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
else    
    disp('Observe! n must be 2 or 3')
    a=NaN;
    P=NaN;
    subplot(1,2,1)
    plot(q,a2,q,a3)
    title('a(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
    subplot(1,2,2)
    plot(q,P2,q,P3)
    title('P(Q), n=3, n=2 failure')
    xlabel('Q')
    grid on
end