I am getting this error below, i do not understand how to fix it.
Show older comments
%{
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 Comments
Walter Roberson
on 24 Nov 2019
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
on 24 Nov 2019
Edited: Bilal Arshed
on 24 Nov 2019
Accepted Answer
Walter Roberson
on 24 Nov 2019
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 Comments
Bilal Arshed
on 24 Nov 2019
Edited: Bilal Arshed
on 24 Nov 2019
Walter Roberson
on 24 Nov 2019
T = interp1(satm(:,1),satm(:,5), k,'makima', 0);
would tell it to use 0 for anything outside the data range in the first variable.
Bilal Arshed
on 24 Nov 2019
Edited: Bilal Arshed
on 24 Nov 2019
Walter Roberson
on 24 Nov 2019
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.
More Answers (0)
Categories
Find more on Aerospace Applications in Help Center and File Exchange
on 24 Nov 2019
on 24 Nov 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
