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

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

[Walter Roberson]

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]

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

Answer 1

[Walter Roberson]

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.

[Bilal Arshed]

i see that the interp1 is pulling -620 for temp, would adding a row of zeros to the bottom of the table help?

Since the space vehicle does not enter the atmosphere until altitude is 86km

or....

Could i try

real(c=(k*R*T)^(1/2);)

[Walter Roberson]

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]

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]

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.