% Control3_2.m
% Solution to problem 2 of control 3
% JMS, May.2024

clear

% Get experimental data
load('DatosExp_C3.mat');
h = Altura';
nh = numel(h);
T = zeros(nh,1);
for ih = 1:nh
  T(ih) = meanValue(T_medidas(:,ih));
  DT(ih) = standDev(T_medidas(:,ih));
end

% Fit data
% h0 = 48.8;                                 % temptative mean height (km)
% b = [ 268; 0.91; -0.18; -1e-2; -2.5e-4 ];  % temptative parameters
h0 = sum(h.*T)/sum(T);
b = regress( T, [ones(size(T)), h-h0, (h-h0).^2, (h-h0).^3, (h-h0).^4]);
Tfit = @(h)( b(1) + b(2)*(h-h0) + b(3)*(h-h0).^2 ...
           + b(4)*(h-h0).^3 + b(5)*(h-h0).^4 );

% Find heights with T>260K
hmesh = min(h):1e-3:max(h);  % a fine mesh of heights (interval of 1 m)
Tmesh = Tfit(hmesh);         % fit temperatures at mesh heights
ih = find(Tmesh>260);        % index of heights with T>260K
hmin = min(hmesh(ih));       % min height with T>260K (with 1m accuracy)
hmax = max(hmesh(ih));       % max height with T>260K

% Plot data and fit
figure(1)
errorbar(h,T,DT,'o'), hold on
plot( hmesh,Tmesh,'b', hmesh(ih),Tmesh(ih),'r', 'LineWidth',2), hold off
xlabel('height (km)')
ylabel('temperature (K)')
grid on
set(gca,'FontSize',15)
