I am getting this error below, i do not understand how to fix it.

2 次查看（过去 30 天）
Bilal Arshed
Bilal Arshed 2019-11-24
评论： Walter Roberson 2019-11-24
%{
Error using griddata (line 85)
Input coordinates cannot be complex.
Error in get_cl (line 50) (scroll to bottom, you will find this)
cl = griddata(alpha_list(:), Ma_list(:), cl_matrix, alpha, Ma);
Error in Lift_Drag (line 3)
cl = get_cl(h, v, alpha);
Error in myeqn (line 18)
[L, D] = Lift_Drag(h, alpha, v, rho);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode113 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in run_myeqn (line 8)
[t,v] = ode113(@myeqn, t0:tend, v0,options);
%}
I get this error from running, "run_myeqn"
%% first script - run_myeqn
clc
clear
v0=[-610.1, 6973.4, 6871e3, 0];
t0=0;
tend=2000; %if you try tend=620 for example i find an error.
options=odeset('Events',@myEvent,'reltol',1e-7,'abstol',1e-7);
[t,v] = ode113(@myeqn, t0:tend, v0,options);
r = v(:,3);
theta = v(:,4);
v_t = v(:,2);
v_r = v(:,1);
%Shuttle Orbit and Landing
figure();
polarplot(theta,r) ; hold on
th = linspace (0,2*pi,50);
r = 6371e3;
polarplot(th, r+zeros(size(th))) ; hold on
polarplot(v(1,4),v(1,3),'r*');
polarplot(v(end,4),v(end,3),'b*');
title('Shuttle Orbit and Landing');
%%second script - myeqn
function dv = myeqn(t, x)
alpha = 0; % starting AoA
%State vector, intial conditions deifined in run_myeqn
theta=x(4); %starting theta angle
r=x(3); %starting postion above earth
v_t=x(2); % tanget velocity
v_r=x(1); % radial velocity
gamma = tan(x(1)./x(2));
v=sqrt(v_t^2+v_r^2); % velocity vector
mu_e = 3.986e14;
Re=6371e3;
h = r-Re;
m=5000;
[T, P, rho] = standard_atm(h);
[L, D] = Lift_Drag(h, alpha, v, rho);
dv_r = (-(mu_e)/(r^2)) + ((v_t^2)/r) + ((1/m)*(-D*sin(gamma)+L*cos(gamma))); % radial accel equation
dv_t= -((v_r*v_t)/r) + 1/m*(-D*cos(gamma) - L*sin(gamma)); % tangent accel equation
dr = v_r;
dtheta = v_t/r;
dv=[dv_r dv_t dr dtheta]';
return
%%script which is giving error
function cl = get_cl(h, v, alpha)
%{
INPUTS
h -- altitude (m)
v -- velocity (m/s)
alpha -- angle of attack (rad)
OUTPUT
cl -- Interpolated coefficient of lift
%}
% Computation of Mach number
k=1.4;
R=287.058; % J/(kg*K)
[T, P, rho] = standard_atm(h); % atmospheric parameters, including T temperature
c=(k*R*T)^(1/2); % speed of sound, m/s
Ma=v/c;
% You can load in a text file, or copy & paste the data as a matrix in here
clfile = [-4 -2 0 2 4 6 8 10 12 14 16 18 20 22.5 25 30 35 40;...
0.3 -0.0982 -0.0222 0.0539 0.1299 0.2059 0.2854 0.3719 0.4649 0.564 0.6691 0.7787 0.8905 1.0336 1.1843 1.4933 1.8025 2.1016;...
0.6 -0.0971 -0.0209 0.0552 0.1313 0.2074 0.287 0.3735 0.4666 0.5658 0.6708 0.7803 0.8928 1.0366 1.188 1.4993 1.8119 2.1143;...
0.9 -0.0922 -0.017 0.0582 0.1333 0.2085 0.2872 0.3729 0.4652 0.5638 0.6684 0.7779 0.8912 1.0367 1.1898 1.5058 1.825 2.1465;...
1.2 -0.0711 -0.0093 0.0504 0.114 0.1856 0.2566 0.3014 0.3403 0.3787 0.4177 0.4578 0.4943 0.5351 0.5746 0.6437 0.693 0.7352;...
1.5 -0.0715 -0.024 0.0219 0.0698 0.1226 0.1773 0.2318 0.2835 0.3332 0.3779 0.4178 0.4541 0.4981 0.5413 0.6178 0.6712 0.7166;...
2 -0.0618 -0.0217 0.0176 0.0581 0.1016 0.1463 0.191 0.2368 0.282 0.3235 0.3619 0.3976 0.4411 0.4841 0.5616 0.6193 0.6702;...
3 -0.0549 -0.0201 0.014 0.0496 0.0864 0.1244 0.1627 0.2018 0.2408 0.2784 0.3148 0.3498 0.3928 0.4357 0.5141 0.575 0.6299;...
5 -0.0498 -0.0192 0.0116 0.0432 0.0761 0.1096 0.1435 0.1777 0.212 0.2466 0.2816 0.316 0.3584 0.4009 0.4793 0.5418 0.5984;...
7.5 -0.0471 -0.0187 0.0104 0.0403 0.0712 0.1026 0.135 0.1683 0.2017 0.2358 0.2706 0.3047 0.3469 0.3892 0.4677 0.5308 0.5883;...
10 -0.0457 -0.0185 0.0098 0.039 0.0686 0.0994 0.1314 0.1645 0.1981 0.2322 0.267 0.3011 0.3433 0.3856 0.4641 0.5273 0.5848;...
15 -0.0448 -0.0184 0.0092 0.0379 0.0666 0.0973 0.1293 0.1623 0.1958 0.2298 0.2644 0.2984 0.3405 0.3829 0.4615 0.5248 0.5825;...
20 -0.0441 -0.0183 0.0087 0.037 0.0651 0.0956 0.1279 0.1608 0.1941 0.2278 0.2622 0.2961 0.3383 0.3806 0.4593 0.5228 0.5807;...
40 -0.0441 -0.0183 0.0087 0.037 0.0651 0.0956 0.1279 0.1608 0.1941 0.2278 0.2622 0.2961 0.3383 0.3806 0.4593 0.5228 0.5807];
% the list of alpha (angle of attack) values are in the first row
% convert from deg -> rad
alpha_list = clfile(1,2:end).*(pi/180);
% the list of Mach values are in the first column, in both cases, you need to skip
% the first entry (NaN)
Ma_list = clfile(2:end,1);
cl_matrix = clfile(2:end, 2:end);
% look up griddata in the help file of matlab to learn how it works
cl = griddata(alpha_list(:), Ma_list(:), cl_matrix, alpha, Ma);
  2 个评论
Walter Roberson
Walter Roberson 2019-11-24
standard_atm is not a defined function. We can grab it from one of your other postings, but there appears to be a lot of overlap between that and the clfile code, so it is not clear that it is intended to be used with that funciton.
Lift_Drag is not a defined function, and you do not appear to have posted code for it.
Bilal Arshed
Bilal Arshed 2019-11-24
编辑：Bilal Arshed 2019-11-24
%i have put the code below for both the
%functions that i have not posted code for
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 = alt/1000;
T=interp1(satm(:,1),satm(:,5), k,'makima');
P=interp1(satm(:,1),satm(:,6), k,'makima');
rho=interp1(satm(:,1),satm(:,6), k,'makima');
end
function [L, D] = Lift_Drag(h, alpha, v, rho)
S=297.15;
cl = get_cl(h, v, alpha);
cd = get_cd(h, v, alpha);
L=0.5*rho*v^2*S*cl;
D=0.5*rho*v^2*S*cd;
end
function cd = get_cd(h, v, alpha)
%{
INTERPOL_cd loads the Drag coefficient cd dataset, and interpolates a
single value of cd for given discrete values of altitude(h), velocity (v)
and angle of attack (alpha)
INPUTS
h -- altitude (m)
v -- velocity (m/s)
alpha -- angle of attack (rad)
OUTPUT
cd -- Interpolated coefficient of Drag
%}
% Computation of Mach number
k=1.4;
R=287.058; % J/(kg*K)
[T, P, rho] = standard_atm(h); % atmospheric parameters, incduding T temperature
c=(k*R*T)^(1/2); % speed of sound, m/s
Ma=v/c;
% You can load in a text file, or copy & paste the data as a matrix in here
cdfile =[-4 -2 0 2 4 6 8 10 12 14 16 18 20 22.5 25 30 35 40;...
0.3 0.0385 0.0325 0.0315 0.0365 0.0481 0.0665 0.0932 0.1293 0.1776 0.2422 0.4052 0.7744 1.3828 1.7971 2.8139 4.0548 5.468;...
0.6 0.0376 0.0325 0.0319 0.0365 0.0475 0.0656 0.0911 0.1256 0.173 0.2363 0.39 0.7176 1.3207 1.7197 2.7047 3.9168 5.3035;...
0.9 0.0401 0.0364 0.0364 0.0404 0.0503 0.0668 0.0894 0.1203 0.1636 0.2213 0.3524 0.6182 1.1713 1.5304 2.4271 3.5466 4.8951;...
1.2 0.1192 0.1173 0.1192 0.1262 0.1379 0.1555 0.1696 0.1872 0.2068 0.2296 0.2551 0.2862 0.3306 0.3775 0.4809 0.6001 0.7246;...
1.5 0.1212 0.119 0.1186 0.1222 0.129 0.14 0.1552 0.1743 0.1949 0.2179 0.2431 0.2741 0.3189 0.3666 0.4722 0.5936 0.7201;...
2 0.1104 0.1064 0.105 0.1073 0.1131 0.1225 0.1357 0.1533 0.1729 0.1956 0.2212 0.2523 0.2968 0.3438 0.4476 0.5671 0.6914;...
3 0.0946 0.0893 0.0877 0.09 0.0961 0.1059 0.1193 0.1367 0.156 0.1787 0.2041 0.2348 0.2783 0.3241 0.4253 0.5423 0.6641;...
5 0.0814 0.0754 0.0741 0.0773 0.0849 0.0958 0.1096 0.1269 0.146 0.1683 0.1931 0.2233 0.2659 0.3107 0.4093 0.5225 0.6398;...
7.5 0.0761 0.07 0.069 0.0732 0.0813 0.0919 0.1055 0.1225 0.1412 0.1634 0.1879 0.2173 0.2589 0.3026 0.3993 0.5117 0.6288;...
10 0.074 0.0679 0.0671 0.0718 0.0798 0.0904 0.1038 0.1206 0.1394 0.1614 0.1857 0.2149 0.256 0.2992 0.3949 0.5068 0.6236;...
15 0.0726 0.0664 0.0657 0.0708 0.0786 0.0892 0.1025 0.1192 0.1379 0.1597 0.1838 0.2127 0.2535 0.2964 0.3917 0.5034 0.6201;...
20 0.0719 0.0655 0.0651 0.0704 0.078 0.0885 0.1019 0.1185 0.137 0.1586 0.1825 0.2112 0.2518 0.2945 0.3894 0.501 0.6178;...
40 0.0719 0.0655 0.0651 0.0704 0.078 0.0885 0.1019 0.1185 0.137 0.1586 0.1825 0.2112 0.2518 0.2945 0.3894 0.501 0.6178];...
% the list of alpha (angle of attack) values are in the first row
% convert from deg -> rad
alpha_list = cdfile(1,2:end).*(pi/180);
% the list of Mach values are in the first column, in both cases, you need to skip
% the first entry (NaN)
Ma_list = cdfile(2:end,1);
cd_matrix = cdfile(2:end, 2:end);
% look up griddata in the help file of matlab to learn how it works
cd = griddata(alpha_list(:), Ma_list(:), cd_matrix, alpha, Ma);

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2019-11-24
Your starting altitude works out to 500 km. Your pressure tables only extend to 86 km. A temperature of about -620 is predicted. Then in
c=(k*R*T)^(1/2);
with k and R being positive and T being negative, k*R*T is negative, and square root of that is complex valued.
You cannot go much past 85 km to avoid this problem.
Also you have
P=interp1(satm(:,1),satm(:,6), k,'makima');
rho=interp1(satm(:,1),satm(:,6), k,'makima');
Notice that the right hand sides are the same, so P and rho will be the same. To be consistent with the comments, you should be using column 7 for rho not column 6.
  4 个评论
显示 2更早的评论
Bilal Arshed
Bilal Arshed 2019-11-24
编辑：Bilal Arshed 2019-11-24
Thanks, i think that did the job. My script seems to be running now.
Now that c=0, my mach number (ma) seems to be inf.
Walter Roberson
Walter Roberson 2019-11-24
space.com :
The exact temperature of the thermosphere can vary substantially, but the average temperature above 180 miles (300 km) is about 800 degrees Fahrenheit (427 degrees Celsius) at solar minimum and 1,700 degrees Fahrenheit (927 degrees Celsius) at solar maximum.
... not 0.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by