airplane01

% model parameters:
g = 9.81;      % gravity in m s^{-2}
v_t = 5.5;   % trim velocity in m s^{-1}   
L_D = 5.2;  % drag coefficient L/D


%%% set initial conditions %%%
v0 = v_t+100;   % start at the trim velocity (or add a delta)
theta0 = -pi/3; % initial angle of trajectory
x0 = 0;     % horizotal position is arbitrary
y0 = 1000;  % initial altitude



% Establish differential equations (theta, v are unknowns)
dv_dt = @(v, theta) -g*sin(theta) - 1/L_D * g/v_t^2 .* v^2;
dtheta_dt = @(v, theta) -g/v*cos(theta) + g/v_t^2 .* v;


% Break these next equations up into tiny time steps starting at t=0 (numerical intergation)
x = @(v, theta) v*cos(theta);
y = @(v, theta) v*sin(theta);


% u = matrix of vectors [v, theta, x, y]
% use Euler method: solve at initial coditions + time step, set that as new initial conditions, repeat


tstep = 1; % set time step for Euler method
t = [0:tstep:100]; % time vector (seconds)



[TOUT,YOUT] = ode23(dv_dt,[0 100],y0);
plot(TOUT, YOUT);
%axis([0 100 0 1200]) % set axis limits
    

% plan: optimize in v, then theta, then v, then theta...