MATLAB: Trying to understand a MATLAB .m file

[borjomi]

I'm trying to create a paper plane simulation in C and I'm trying to convert MATLAB code to C. Admittedly, I don't know MATLAB (and it's honestly been a while since I've done a lot of math), but I've been able to decode some parts of the following MATLAB code that I found online:

S		=	0.017;			% Reference Area, m^2
AR		=	0.86;			% Wing Aspect Ratio
e		=	0.9;			% Oswald Efficiency Factor;
m		=	0.003;			% Mass, kg
g		=	9.8;			% Gravitational acceleration, m/s^2
rho		=	1.225;			% Air density at Sea Level, kg/m^3	
CLa		=	3.141592 * AR/(1 + sqrt(1 + (AR / 2)^2));
							% Lift-Coefficient Slope, per rad
CDo		=	0.02;			% Zero-Lift Drag Coefficient
epsilon	=	1 / (3.141592 * e * AR);% Induced Drag Factor	
CL		=	sqrt(CDo / epsilon);	% CL for Maximum Lift/Drag Ratio
CD		=	CDo + epsilon * CL^2;	% Corresponding CD
LDmax	=	CL / CD;			% Maximum Lift/Drag Ratio
Gam		=	-atan(1 / LDmax);	% Corresponding Flight Path Angle, rad
V		=	sqrt(2 * m * g /(rho * S * (CL * cos(Gam) - CD * sin(Gam))));
							% Corresponding Velocity, m/s
Alpha	=	CL / CLa;			% Corresponding Angle of Attack, rad
	
%	Oscillating Glide due to Zero Initial Flight Path Angle
	xo		=	[V;0;H;R];
	[tb,xb]	=	ode23('EqMotion',tspan,xo);

function xdot = EqMotion(t,x)
%	Fourth-Order Equations of Aircraft Motion

	global CL CD S m g rho
	
	V 	=	x(1);
	Gam	=	x(2);
	q	=	0.5 * rho * V^2;	% Dynamic Pressure, N/m^2
	
	xdot	=	[(-CD * q * S - m * g * sin(Gam)) / m
				 (CL * q * S - m * g * cos(Gam)) / (m * V)
				 V * sin(Gam)
				 V * cos(Gam)];

I've figured out that the 'xdot' function is the key to all of this and I think the first line of xdot (i.e. -CD *q * S, etc) represents the velocity, the second line is the gamma (or flight path angle), the third line is the x position and the last line is the y position.

I've also learned how to solve differential equations using the Runge-Kutta Fourth Order Method. However, when I try to solve for the next timestep using RK4, the values do not correspond with the graphs in MATLAB, which means I'm doing something wrong.

So my questions are, first, am I correct in assuming that the four lines of xdot correspond to velocity, gamma, x position, and y position respectively? If this is true, can I use the velocity and gamma equations to solve for v(n+1) and g(n+1) using the Runge-Kutta method?

I'm sorry for my lack of knowledge, but I hope you can help in answering my question. Please feel free to ask any questions and I'll do my best to clarify.

 

[jedishrfu]

The ode23 function is defined as:

https://www.mathworks.com/help/matlab/ref/ode23.html?s_tid=srchtitle
From MATLAB the ode23 function takes the following inputs:
- the time derivative of your function
- the range of time values as an array of values
- the initial starting condition as a tuple of values

and returns the following tuple of values:
- the time array of values
- the x array of values

In your case, I don't see time ie the 't' variable being used in your equations at all as shown in your EqMotion() function. It's hard to proceed from here.

I think we need to see your four equations of motion so we can compare to how they were coded into the EqMotion() function that will be called repeatedly by the ode23 function. Basically the x input tuple must match the xdot output tuple in meaning of tuple elements.

x0 = (x0, x0', x0'') initial values of state

then the EqMotion will work with these state values to create a new state xdot tuple with the same
meaning of elements namely (x, x', x''). The new state will be applied as input to EqMotion on the next ode23 iteration until all t values are exhausted and ode23 ends returning the final answer as a time array and an associated x array.

Here's some more info on how to use the ode23 function in matlab:

http://www.eng.auburn.edu/~tplacek/courses/3600/ode45waterloo.pdf
then the output