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

error

%{ 
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);

Best Answer

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.
Related Question