% Script to plot and analyize Mechatronics 1st Order lab data

clear all, close all;

% Initialize knowns
Cp = 4218; % J/(kg.K)
m1 = 0.090;     % insert your value here

% load 1st data set (known mass)
load thermodata
t1 = thermodata(:,1);
T1 = thermodata(:,2)*1000;
T1i = min(T1);

% use log regression to find Tau1 and iterate to find T1f
T1f = max(T1)+10; % initial guess
dT = 0.05;
sumSqrOld = 1e10;
sumSqr = 1e9;
while(sumSqrOld-sumSqr>=0),
    sumSqrOld = sumSqr;
    y1 = log(T1f-T1);
    [P] = polyfit(t1,y1,1);
    Tau1 = -1/P(1);
    mod1 = T1i+(T1f-T1i)*(1-exp(-t1/Tau1));
    resid = mod1-T1;
    sumSqr = sum(resid.*resid);
    T1f = T1f-dT;
end

% plot data & model to check fit
plot(t1,T1,t1,mod1);

% output qs, Tau1, and T1f
Tau1
T1f
qs = T1f*m1*Cp/Tau1

% load 1st data set (known mass)
load thermodata2
t2 = thermodata2(:,1);
T2 = thermodata2(:,2);
T2i = min(T2);

% use log regression to find Tau2 and iterate to find T2f
T2f = max(T2)+10; % initial guess
sumSqrOld = 1e10;
sumSqr = 1e9;
while(sumSqrOld-sumSqr>=0),
    sumSqrOld = sumSqr;
    y2 = log(T2f-T2);
    [P] = polyfit(t2,y2,1);
    Tau2 = -1/P(1);
    mod2 = T2i+(T2f-T2i)*(1-exp(-t2/Tau2));
    resid = mod2-T2;
    sumSqr = sum(resid.*resid);
    T2f = T2f-dT;
end

% plot data & model to check fit
figure(2)
plot(t2,T2,t2,mod2);

% output qs, Tau1, and T1
Tau2
T2f
m2 = Tau2*qs/(T2f*Cp)