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
A(:,1) = [k(1)*ones(6,1); k(2)*ones(6,1); k(3)*ones(4,1); k(4)*ones(7,1)];
A(:,2)=ones(23,1);

b = e1;

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

clear t0
t0 = a(1);
c0 = a(2);

%
clear A, b, a, y, e1
g = 10e2;
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]';
e1 = y-c0*ones(23,1);
t = [(0:30:150)'; (30:30:180)'; (0:60:180)'; (0:30:180)'];
T = [463.15*ones(6,1); 473.15*ones(6,1); 483.15*ones(11,1)];
x = [2.15*ones(6,1); 2.95*ones(6,1); 2.45*ones(4,1); 2.65*ones(7,1)];

A = [ones(23, 1) g./T log(x)];
b = log(-e1)-log(t+t0*ones(23, 1));

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

e = A*a-b;

k0 = -exp(a(1));
dE = -a(2)*8.314*g;
bt = a(3);


e=y-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1)
e'*e
clear A, a, b
%A(:,1) = k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1));
%A(:,2) = -k0./(8.314*T).*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))*10^3;
%A(:,3) = k0*exp(-dE./(8.314*T)).*log(x).*(x.^bt).*(t+t0*ones(23, 1));
%A(:,4) = k0*exp(-dE./(8.314*T)).*log(x).*(x.^bt)*10^2;
%A(:,4) = ones(23, 1);

%b = e;
%a = inv(A'*A)*A'*b
%e = A*a-b
af = 0;
clear de1, e1

de1 = 1;
e1 = e'*e;
while (de1>1e-7)

A = log(y).*y.^(1+af)/(1+af);
b = -e;

daf = inv(A'*A)*A'*b;
af = daf+af;

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

%correction of k0
A = k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1));
b = e;

dk0 = inv(A'*A)*A'*b;
k0 = k0*(1+dk0);

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

%correction of dE
A = -1./(8.314*T).*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))*k0;
b = e;

ddE = inv(A'*A)*A'*b;
dE = dE+ddE;

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

%correction of bt
A = k0*exp(-dE./(8.314*T)).*log(x).*(x.^bt).*(t+t0*ones(23, 1));
b = e;

dbt = inv(A'*A)*A'*b;
bt = bt+dbt;

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

%correction of t0
A = k0*exp(-dE./(8.314*T)).*log(x).*(x.^bt);
b =e;

dt0 = inv(A'*A)*A'*b;
t0 = t0+dt0;

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

%correction of c0
A = ones(23, 1);
b = e;

dc0 = inv(A'*A)*A'*b;
c0 = c0+dc0;

e=y.^(1+af)/(1+af)-k0*exp(-dE./(8.314*T)).*(x.^bt).*(t+t0*ones(23, 1))-c0*ones(23, 1);

de1 = e1-e'*e
e1 = e'*e;

end