How to debug this Matlab code to model glider descent?

Question

I have this Matlab project but for some reason I just cannot stop thinking about it because I could not get it to work.

Objective:

This is a MATLAB script that would calculate the change of pressure, temperature and density of a glider that is being dropped from 10000 feet. As it falls, we want to use those new values calculated and then plugged in a function that has 4 equations that need to be differentiated at every point using ode45 as well as the new values of P T and Rho.

Here is the main code:

% HouseKeeping:

clc
clear all
close all

% Constants:
S = 232;                                         % ft^2
Cd0 = 0.02;
K = 0.07;
W = 11000;                                       % lbf
Cl_max = sqrt(Cd0/K);
Cd_max = 2*K*Cl_max^2;
Rho_10000 = .001756;                             % slugs/ ft^3

%Initial conditions:
t = 0;                                           % Sec
x = 0;                                           % ft
h = 10000;                                       % ft
v = sqrt((2*W)/(Rho_10000*S*Cl_max));            %ft/s
gamma = -Cd_max/Cl_max; 

% Find Endurance:
V_RD= sqrt((2*W)/(S* Rho_10000* sqrt(3*Cd0/K)));
RD= V_RD/((sqrt(3*Cd0/K))/(2*Cd0))  ;            % ft/s
Endurcance= h/RD;                                % 958.3515 sec

% Sea Level values:
TSL = 518.69;                                    % Rankine
PSL = 2116.199414;                               % lbs/ft^2
RhoSL = 0.0023769;                          % slugs/ft^3

while h > 0

    tspan = [t t+1];
    i=1;
    X = [x;h;v;gamma;Rho_10000];
    Time(i)= t;

    % Calling ODE45:
    [F] = ode45(@ D,tspan, X)

    % Hight Varying Parameters:
    T = TSL - 0.00356616*h;
    P = (1.137193514E-11)*(T)^5.2560613;
    Rho = (RhoSL * TSL / PASL)*(P / T);
    a = 49.0214 * (T)^.5;

    H_Del(i) = (-Cd_max/Cl_max)*(plotted_x(i))+10000;
    V_Del(i) = sqrt((2*W)/(Rho*S*Cl_max));
    Gamma_Del(i) = -Cd_max/Cl_max;
    i= i+1;

    X = [ x ; H_Del(i); V_Del(i); Gamma_Del(i); Rho];

end



% Plots:

%1
figure (1)
plot(F(:,1),'-r',F(:,3),'-b')
title('Velocity vs Distance');
xlabel('x (ft)');
ylabel('v (ft/s)');


%2 
Figure (2)
plot(F(:,1),'-r',F(:,2),'-b')
title('Altitude vs Distance ');
xlabel('x (ft)');
ylabel('h (ft)');


%3
figure (3)
plot(F(:,1),'-r',F(:,4),'-b')
title('Gamma vs Distance');
xlabel('x (ft)');
ylabel('Gamma (rad)');

%4
figure (4)
plot(t,'-r',F(:,3),'-b')
title('Velocity vs Time');
xlabel(' t (s)');
ylabel('v  (ft/s)');

%5
figure (5)
plot(t,'-r',F(:,1),'-b')
title(' Distance vs Time ');
xlabel('t (s)');
ylabel('x (ft)');


%6
figure (6)
plot(t,'-r',F(:,4),'-b')
title('Gamma vs Time ');
xlabel('t (s)');
ylabel('Gamma (rad)');

%7 
figure (7)
plot(t,'-r',F(:,3),'-b')
title('Velocity vs Time');
xlabel('t (s)');
ylabel('V (ft/s)');

Here is the Function

function [F] = D(X)

% Constants:
S = 232;                                         % ft^2
Cd0 = 0.02;
K = 0.07;
W = 11000;                                       % lbf
Cl_max = sqrt(Cd0/K);
Cd_max = 2*K*Cl_max^2;
Rho_10000 = .001756;                             % slugs/ ft^3

% Initial conditions:
t = 0;                                           % Sec
x = 0;                                           % ft
h = 10000;                                       % ft
v = sqrt((2*W)/(Rho_10000*S*Cl_max));            % ft/s
gamma = -Cd_max/Cl_max; 
g= 32.2;                                         % ft/s^2

% Equations:
X_Pr = X(3)*cos(X(4));
H_Pr = X(3)*sin(X(4));
V_Pr = (-32.2./W)*(0.5*X(5)*Cd_max*S*X(3)^2 + W*sin(X(4)));
G_Pr = (32.2./(W*X(3)))*(0.5*X(5)*Cl_max*S*X(3)^2 - W*cos(X(4)));

F = [X_Pr;H_Pr;V_Pr;G_Pr];

I am not very good with Matlab but I did my very best! I went to my professors for help but they said they were too busy. I even stalked every senior I knew and they all said they did not know how to do it. My next project will be assigned soon and I think that if I cannot do this then I would not be able to do the next one.

Answer 1

Your code produces the following error:

Error using main>D

Too many input arguments.

This means that ode45 tries to call your provided function D with too many input arguments. You should check the desired odefun format in the ode45 documentation: dydt = odefun(t,y)

So, you should change your function declaration of D to

function [F] = D(t, X)

This should solve your first problem, but a following error pops up:

D returns a vector of length 4, but the length of initial conditions vector is 5. The vector returned by D and the initial conditions vector must have the same number of elements.

Again, you should check the ode45 documentation. Your function should return the derivatives of all your input variables X: F = dX/dt. You should also return the derivate of the fifth element Rho_10000.

Next, I got some error about undefined variables such as PASL. Probably because you did not post your full code.

Besides of the errors, you should really check your code again. You have written an infinite while loop while h > 0. You never change h in the loop, nor you use the output of ode45 in your loop. Furthermore you always overwrite your i and X value at the beginning of the loop, which is probably not what you want.

This is not a full answer to your question, but I hope you will be able to continue and post smaller, well defined problems, instead of one big problem which is very difficult to answer completely.