% Control3_2.m
% Solution to problem 2 of control 3, May.2018
% J.M.Soler, Apr.2018

% Set data
clear
l0 = [0.175;0.180;0.185];  % equil. spring lengths (m)
k0 = [290;250;275];        % spring const. at T=0 (N/m)
D = [20;16;100];           % temperature correction spring const (N/m)
T0 = [100;200;350];        % spring temperature parameters (K)
m = [0.10;0.15];           % hanging masses (kg)
L = 0.65;                  % sum of spring lengths (m)

% Find temperature-dependent spring lengths
T = 0:1:1000;
nT = numel(T);
for iT = 1:nT
    d(:,iT) = equilibrium_3_springs(k0, D, T0, l0, m, L, T(iT));
end
l = repmat(l0,1,nT)+d;     % spring lengths

% Find T for which la=lb
dl = l(1,:)-l(2,:);        % difference of spring lengths
iT = find(dl>0,1)-1;       % last temperature index at which dl<0
Tab = (T(iT)*dl(iT+1)-T(iT+1)*dl(iT)) / (dl(iT+1)-dl(iT))  % lever rule

% Plot spring lengths, and Tab
figure(1)
lmin = min(min(l));
lmax = max(max(l));
plot(T,l(1,:),T,l(2,:),T,l(3,:),[Tab,Tab],[lmin,lmax])
legend('a','b','c')
xlabel('Temperature (K)')
ylabel('spring lengths (m)')
