% Plots the residuals for two possible altitude models to see which one is better --------------
d = [0 1 2 3 4 5 6 7 8 9 10; % Altitude (km)
1 .90748 .82162 .74214 .66868 .60091 .53853 .48123 .42871 .38069 .33690]; % density ratio
% The data above is from the 1976 US Standard Atmosphere, US Government Printing Office:
% https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf
% starting on page 67 of the pdf document (last column of the table)
ktf = 1000/(25.4 * 12); % convert km to thousands of feet
a = d(1,:) * ktf; % altitude in thousands of feet
s = d(2,:); % density ratio
alt = (0:.01:10) * ktf; % x axis for extrapolated values (every 10 meters)
k = 1:100:1001; % indices into alt for every 1000 meters
tid = {'Calc' 'True' 'Resid'}; xlbl = 'Thousands of feet';
% The density ratio must decay to 0 at infinite altitude and is defined to be 1.0 at sea level.
% The simplest model satisfying these conditions is the single parameter model: e ^ (-alt/p)
p = fminsearch(@(p,a,s) norm(exp(-a/p) - s),1,[],a,s);
t = prin('sigma = e ^ (-alt/%6w)',p); % show the model in the figure name
S = exp(-alt/p); % extrapolate density ratio from model
plt(alt,S,a,s,a,(S(k)-s)*100,'Ylabel',{'Density ratio' 'residual (%)'},...
'FigName',t,'Xlabel',xlbl,'Styles','-ox','MarkerSize',12,'TraceID',tid);
% This two-parameter model also satisfies those conditions: (1-alt/b) ^ c --------------------
p = fminsearch(@(p,a,s) norm((1-a/p(1)).^p(2) - s),[100 1],[],a,s);
t = prin('sigma = (1 - alt/%6w) ^ %6w',p); % show the model in the figure name
S = (1-alt/p(1)) .^ p(2); % extrapolate density ratio from model
plt(alt,S,a,s,a,(S(k)-s)*1e6,'Ylabel',{'Density ratio' 'residual (ppm)'},...
'FigName',t,'Xlabel',xlbl,'Styles','-ox','MarkerSize',12,'TraceID',tid);