clear
%first group t0 estimation
y=[32.02	31.44	30.75	30.17	29.69	29.11]';
t = (0:30:150)';

A = [y ones(length(y), 1)];

b = t;
a = inv(A'*A)*A'*b;

k(1) = 1/a(1);
t0(1) = a(2);

e1 = y-k(1)*t;

%second group t0 estimation
y=[ 31.43	30.96	30.35	29.52	28.87	28.11]';
t = (30:30:180)';

A = [y ones(length(y), 1)];

b = t;
a = inv(A'*A)*A'*b;

k(2) = 1/a(1);
t0(2) = a(2);

e1 = [e1; y-k(2)*t];
%third group t0 estimation
y=[ 31.34	28.83	26.32	24.26]';
t = (0:60:180)';

A = [y ones(length(y), 1)];

b = t;
a = inv(A'*A)*A'*b;

k(3) = 1/a(1);
t0(3) = a(2);

e1 = [e1; y-k(3)*t];
%fourth group t0 estimation
y=[32.13	31.20	30.23	29.35	27.98	27.07	26.40]';
t = (0:30:180)';

A = [y ones(length(y), 1)];

b = t;
a = inv(A'*A)*A'*b;

k(4) = 1/a(1)
t0(4) = a(2)

e1 = [e1; y-k(4)*t];

clear A, b, a, y, e1
y = [32.02	31.44	30.75	30.17	29.69	29.11 31.43	30.96	30.35	29.52	28.87	28.11 31.34	28.83	26.32	24.26  32.13	31.20	30.23	29.35	27.98	27.07	26.40]';
t = [(0:30:150)'; (30:30:180)'; (0:60:180)'; (0:30:180)'];
K = [k(1)*ones(6,1); k(2)*ones(6,1); k(3)*ones(4,1); k(4)*ones(7,1)];
b = ones(23,1);

a = zeros(100,1);
for i=1:100
    af = (i-50)*0.02;
    A = (1/(1+af))./K.*y.^(1+af)-t;
    a(i) = A'*b/(A'*A);
end