Could someone please help me setup an interpolation script to find atmospheric conditions at various altitudes?

2 views (last 30 days)
Bilal Arshed
Bilal Arshed on 23 Nov 2019
Edited: dpb on 24 Nov 2019
Hi, i am trying to interpolate T, P, rho for any alt (altitude). Currently it grabs the closest values but is not good for number such as 85.6e3 (m) as this part of the code
k = find(satm(:,1) <= alt/1000, 1, 'last');
takes the last value. For which 86km is not included so 86km is not taken as the final value which is the value it should be taking.
Thanks for any help in advance!
function [T, P, rho] = standard_atm(alt)
%{
Function implementing the International Standard Atmospheric model, which
relies on look-up tables.
Currently, the function looks up the closest altitude in the table to the
one asked for in the input (alt), and returns that value. A better method
would be to assume a linear or polynomial interpolation for intermediate
values.
INPUTS
alt = altitude above Earth surface, m
OUTPUTS
T = atmospheric temperature (K)
P = atmospheric pressure (Pa)
rho = atmospheric density (kg/m3)
%}
% make sure you are working in the correct units! The table below uses km,
% not m for altitude
% alt sigma delta theta temp press dens a visc k.visc
% km ` K N/sq.m kg/cu.m m/s kg/m-s sq.m/s
satm = [-2 1.2067E+0 1.2611E+0 1.0451 301.2 1.278E+5 1.478E+0 347.9 18.51 1.25E-5;...
0 1.0000E+0 1.0000E+0 1.0000 288.1 1.013E+5 1.225E+0 340.3 17.89 1.46E-5;...
2 8.2168E-1 7.8462E-1 0.9549 275.2 7.950E+4 1.007E+0 332.5 17.26 1.71E-5;...
4 6.6885E-1 6.0854E-1 0.9098 262.2 6.166E+4 8.193E-1 324.6 16.61 2.03E-5;...
6 5.3887E-1 4.6600E-1 0.8648 249.2 4.722E+4 6.601E-1 316.5 15.95 2.42E-5;...
8 4.2921E-1 3.5185E-1 0.8198 236.2 3.565E+4 5.258E-1 308.1 15.27 2.90E-5;...
10 3.3756E-1 2.6153E-1 0.7748 223.3 2.650E+4 4.135E-1 299.5 14.58 3.53E-5;...
12 2.5464E-1 1.9146E-1 0.7519 216.6 1.940E+4 3.119E-1 295.1 14.22 4.56E-5;...
14 1.8600E-1 1.3985E-1 0.7519 216.6 1.417E+4 2.279E-1 295.1 14.22 6.24E-5;...
16 1.3589E-1 1.0217E-1 0.7519 216.6 1.035E+4 1.665E-1 295.1 14.22 8.54E-5;...
18 9.9302E-2 7.4662E-2 0.7519 216.6 7.565E+3 1.216E-1 295.1 14.22 1.17E-4;...
20 7.2578E-2 5.4569E-2 0.7519 216.6 5.529E+3 8.891E-2 295.1 14.22 1.60E-4;...
22 5.2660E-2 3.9945E-2 0.7585 218.6 4.047E+3 6.451E-2 296.4 14.32 2.22E-4;...
24 3.8316E-2 2.9328E-2 0.7654 220.6 2.972E+3 4.694E-2 297.7 14.43 3.07E-4;...
26 2.7964E-2 2.1597E-2 0.7723 222.5 2.188E+3 3.426E-2 299.1 14.54 4.24E-4;...
28 2.0470E-2 1.5950E-2 0.7792 224.5 1.616E+3 2.508E-2 300.4 14.65 5.84E-4;...
30 1.5028E-2 1.1813E-2 0.7861 226.5 1.197E+3 1.841E-2 301.7 14.75 8.01E-4;...
32 1.1065E-2 8.7740E-3 0.7930 228.5 8.890E+2 1.355E-2 303.0 14.86 1.10E-3;...
34 8.0709E-3 6.5470E-3 0.8112 233.7 6.634E+2 9.887E-3 306.5 15.14 1.53E-3;...
38 4.3806E-3 3.7218E-3 0.8496 244.8 3.771E+2 5.366E-3 313.7 15.72 2.93E-3;...
40 3.2615E-3 2.8337E-3 0.8688 250.4 2.871E+2 3.995E-3 317.2 16.01 4.01E-3;...
42 2.4445E-3 2.1708E-3 0.8880 255.9 2.200E+2 2.995E-3 320.7 16.29 5.44E-3;...
44 1.8438E-3 1.6727E-3 0.9072 261.4 1.695E+2 2.259E-3 324.1 16.57 7.34E-3;...
46 1.3992E-3 1.2961E-3 0.9263 266.9 1.313E+2 1.714E-3 327.5 16.85 9.83E-3;...
48 1.0748E-3 1.0095E-3 0.9393 270.6 1.023E+2 1.317E-3 329.8 17.04 1.29E-2;...
50 8.3819E-4 7.8728E-4 0.9393 270.6 7.977E+1 1.027E-3 329.8 17.04 1.66E-2;...
52 6.5759E-4 6.1395E-4 0.9336 269.0 6.221E+1 8.055E-4 328.8 16.96 2.10E-2;...
54 5.2158E-4 4.7700E-4 0.9145 263.5 4.833E+1 6.389E-4 325.4 16.68 2.61E-2;...
56 4.1175E-4 3.6869E-4 0.8954 258.0 3.736E+1 5.044E-4 322.0 16.40 3.25E-2;...
58 3.2344E-4 2.8344E-4 0.8763 252.5 2.872E+1 3.962E-4 318.6 16.12 4.07E-2;...
60 2.5276E-4 2.1668E-4 0.8573 247.0 2.196E+1 3.096E-4 315.1 15.84 5.11E-2;...
62 1.9647E-4 1.6468E-4 0.8382 241.5 1.669E+1 2.407E-4 311.5 15.55 6.46E-2;...
64 1.5185E-4 1.2439E-4 0.8191 236.0 1.260E+1 1.860E-4 308.0 15.26 8.20E-2;...
66 1.1668E-4 9.3354E-5 0.8001 230.5 9.459E+0 1.429E-4 304.4 14.97 1.05E-1;...
68 8.9101E-5 6.9593E-5 0.7811 225.1 7.051E+0 1.091E-4 300.7 14.67 1.34E-1;...
70 6.7601E-5 5.1515E-5 0.7620 219.6 5.220E+0 8.281E-5 297.1 14.38 1.74E-1;...
72 5.0905E-5 3.7852E-5 0.7436 214.3 3.835E+0 6.236E-5 293.4 14.08 2.26E-1;...
74 3.7856E-5 2.7635E-5 0.7300 210.3 2.800E+0 4.637E-5 290.7 13.87 2.99E-1;...
76 2.8001E-5 2.0061E-5 0.7164 206.4 2.033E+0 3.430E-5 288.0 13.65 3.98E-1;...
78 2.0597E-5 1.4477E-5 0.7029 202.5 1.467E+0 2.523E-5 285.3 13.43 5.32E-1;...
80 1.5063E-5 1.0384E-5 0.6893 198.6 1.052E+0 1.845E-5 282.5 13.21 7.16E-1;...
82 1.0950E-5 7.4002E-6 0.6758 194.7 7.498E-1 1.341E-5 279.7 12.98 9.68E-1;...
84 7.9106E-6 5.2391E-6 0.6623 190.8 5.308E-1 9.690E-6 276.9 12.76 1.32E+0;...
86 5.6777E-6 3.6835E-6 0.6488 186.9 3.732E-1 6.955E-6 274.1 12.53 1.80E+0];
k = find(satm(:,1) <= alt/1000, 1, 'last');
T = satm(k,5);
P = satm(k,6);
rho = satm(k,7);
end

Sign in to answer this question.

Accepted Answer

dpb
dpb on 23 Nov 2019
Edited: dpb on 24 Nov 2019
function [T, P, rho] = standard_atm(alt)
...
% alt sigma delta theta temp press dens a visc k.visc
% km ` K N/sq.m kg/cu.m m/s kg/m-s sq.m/s
satm = [-2 1.2067E+0 1.2611E+0 1.0451 301.2 1.278E+5 1.478E+0 347.9 18.51 1.25E-5;...
...
84 7.9106E-6 5.2391E-6 0.6623 190.8 5.308E-1 9.690E-6 276.9 12.76 1.32E+0;...
86 5.6777E-6 3.6835E-6 0.6488 186.9 3.732E-1 6.955E-6 274.1 12.53 1.80E+0];
% insert error checking code for range here if desired..
...
vq=interp1(satm(:,1),satm(:,[5:7]),alt); % linear interpolation by default
T = vq(1);
P = vq(2);
rho = vq(3);
end
Pick whatever level/form of interpolation wanted...default is linear, can use one of several alternatives. See doc interp1 for details.
  2 Comments
Bilal Arshed
Bilal Arshed on 23 Nov 2019
im trying this with my code but i get NAN..
function [T, P, rho] = standard_atm(alt)
%{
Function implementing the International Standard Atmospheric model, which
relies on look-up tables.
Currently, the function looks up the closest altitude in the table to the
one asked for in the input (alt), and returns that value. A better method
would be to assume a linear or polynomial interpolation for intermediate
values.
INPUTS
alt = altitude above Earth surface, m
OUTPUTS
T = atmospheric temperature (K)
P = atmospheric pressure (Pa)
rho = atmospheric density (kg/m3)
%}
% make sure you are working in the correct units! The table below uses km,
% not m for altitude
% alt sigma delta theta temp press dens a visc k.visc
% km ` K N/sq.m kg/cu.m m/s kg/m-s sq.m/s
satm = [-2 1.2067E+0 1.2611E+0 1.0451 301.2 1.278E+5 1.478E+0 347.9 18.51 1.25E-5;...
0 1.0000E+0 1.0000E+0 1.0000 288.1 1.013E+5 1.225E+0 340.3 17.89 1.46E-5;...
2 8.2168E-1 7.8462E-1 0.9549 275.2 7.950E+4 1.007E+0 332.5 17.26 1.71E-5;...
4 6.6885E-1 6.0854E-1 0.9098 262.2 6.166E+4 8.193E-1 324.6 16.61 2.03E-5;...
6 5.3887E-1 4.6600E-1 0.8648 249.2 4.722E+4 6.601E-1 316.5 15.95 2.42E-5;...
8 4.2921E-1 3.5185E-1 0.8198 236.2 3.565E+4 5.258E-1 308.1 15.27 2.90E-5;...
10 3.3756E-1 2.6153E-1 0.7748 223.3 2.650E+4 4.135E-1 299.5 14.58 3.53E-5;...
12 2.5464E-1 1.9146E-1 0.7519 216.6 1.940E+4 3.119E-1 295.1 14.22 4.56E-5;...
14 1.8600E-1 1.3985E-1 0.7519 216.6 1.417E+4 2.279E-1 295.1 14.22 6.24E-5;...
16 1.3589E-1 1.0217E-1 0.7519 216.6 1.035E+4 1.665E-1 295.1 14.22 8.54E-5;...
18 9.9302E-2 7.4662E-2 0.7519 216.6 7.565E+3 1.216E-1 295.1 14.22 1.17E-4;...
20 7.2578E-2 5.4569E-2 0.7519 216.6 5.529E+3 8.891E-2 295.1 14.22 1.60E-4;...
22 5.2660E-2 3.9945E-2 0.7585 218.6 4.047E+3 6.451E-2 296.4 14.32 2.22E-4;...
24 3.8316E-2 2.9328E-2 0.7654 220.6 2.972E+3 4.694E-2 297.7 14.43 3.07E-4;...
26 2.7964E-2 2.1597E-2 0.7723 222.5 2.188E+3 3.426E-2 299.1 14.54 4.24E-4;...
28 2.0470E-2 1.5950E-2 0.7792 224.5 1.616E+3 2.508E-2 300.4 14.65 5.84E-4;...
30 1.5028E-2 1.1813E-2 0.7861 226.5 1.197E+3 1.841E-2 301.7 14.75 8.01E-4;...
32 1.1065E-2 8.7740E-3 0.7930 228.5 8.890E+2 1.355E-2 303.0 14.86 1.10E-3;...
34 8.0709E-3 6.5470E-3 0.8112 233.7 6.634E+2 9.887E-3 306.5 15.14 1.53E-3;...
38 4.3806E-3 3.7218E-3 0.8496 244.8 3.771E+2 5.366E-3 313.7 15.72 2.93E-3;...
40 3.2615E-3 2.8337E-3 0.8688 250.4 2.871E+2 3.995E-3 317.2 16.01 4.01E-3;...
42 2.4445E-3 2.1708E-3 0.8880 255.9 2.200E+2 2.995E-3 320.7 16.29 5.44E-3;...
44 1.8438E-3 1.6727E-3 0.9072 261.4 1.695E+2 2.259E-3 324.1 16.57 7.34E-3;...
46 1.3992E-3 1.2961E-3 0.9263 266.9 1.313E+2 1.714E-3 327.5 16.85 9.83E-3;...
48 1.0748E-3 1.0095E-3 0.9393 270.6 1.023E+2 1.317E-3 329.8 17.04 1.29E-2;...
50 8.3819E-4 7.8728E-4 0.9393 270.6 7.977E+1 1.027E-3 329.8 17.04 1.66E-2;...
52 6.5759E-4 6.1395E-4 0.9336 269.0 6.221E+1 8.055E-4 328.8 16.96 2.10E-2;...
54 5.2158E-4 4.7700E-4 0.9145 263.5 4.833E+1 6.389E-4 325.4 16.68 2.61E-2;...
56 4.1175E-4 3.6869E-4 0.8954 258.0 3.736E+1 5.044E-4 322.0 16.40 3.25E-2;...
58 3.2344E-4 2.8344E-4 0.8763 252.5 2.872E+1 3.962E-4 318.6 16.12 4.07E-2;...
60 2.5276E-4 2.1668E-4 0.8573 247.0 2.196E+1 3.096E-4 315.1 15.84 5.11E-2;...
62 1.9647E-4 1.6468E-4 0.8382 241.5 1.669E+1 2.407E-4 311.5 15.55 6.46E-2;...
64 1.5185E-4 1.2439E-4 0.8191 236.0 1.260E+1 1.860E-4 308.0 15.26 8.20E-2;...
66 1.1668E-4 9.3354E-5 0.8001 230.5 9.459E+0 1.429E-4 304.4 14.97 1.05E-1;...
68 8.9101E-5 6.9593E-5 0.7811 225.1 7.051E+0 1.091E-4 300.7 14.67 1.34E-1;...
70 6.7601E-5 5.1515E-5 0.7620 219.6 5.220E+0 8.281E-5 297.1 14.38 1.74E-1;...
72 5.0905E-5 3.7852E-5 0.7436 214.3 3.835E+0 6.236E-5 293.4 14.08 2.26E-1;...
74 3.7856E-5 2.7635E-5 0.7300 210.3 2.800E+0 4.637E-5 290.7 13.87 2.99E-1;...
76 2.8001E-5 2.0061E-5 0.7164 206.4 2.033E+0 3.430E-5 288.0 13.65 3.98E-1;...
78 2.0597E-5 1.4477E-5 0.7029 202.5 1.467E+0 2.523E-5 285.3 13.43 5.32E-1;...
80 1.5063E-5 1.0384E-5 0.6893 198.6 1.052E+0 1.845E-5 282.5 13.21 7.16E-1;...
82 1.0950E-5 7.4002E-6 0.6758 194.7 7.498E-1 1.341E-5 279.7 12.98 9.68E-1;...
84 7.9106E-6 5.2391E-6 0.6623 190.8 5.308E-1 9.690E-6 276.9 12.76 1.32E+0;...
86 5.6777E-6 3.6835E-6 0.6488 186.9 3.732E-1 6.955E-6 274.1 12.53 1.80E+0;...
87 0 0 0 0 0 0 0 0 0];
k = find(satm(:,1) <= alt/1000, 1, 'last');
T = satm(k,5);
P = satm(k,6);
rho = satm(k,7);
vq=interp1(satm(:,1),satm(:,[5:7])); % linear interpolation by default
T = vq(1);
P = vq(2);
rho = vq(3);
end
dpb
dpb on 24 Nov 2019
Edited: dpb on 24 Nov 2019
You should get rid of the first code section using find no purpose in it.
You don't show the calling statement but must be passing in a value outside range of table then... interp1 won't extrapolate by default; returns NaN instead.
Why I suggested error-checking code in front.
Oh! I see I didn't add the lookup value to the call to interp1() when cut/pasted...
vq=interp1(satm(:,1),satm(:,[5:7]),alt(:)); % linear interpolation by default
NB: vq will be array if alt is a vector for multiple lookups at once. In that case use
T = vq(:,1);
P = vq(:,2);
rho = vq(:,3);
to return the vectors.
COMMENT: The row of zeros at the end of the table surely looks funky...that surely isn't from the actual model dataset is it???

Sign in to comment.

More Answers (0)

Categories

Find more on Interpolation in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!