rocket.m

1; %Need this to make sure Matlab/Octave reads it as a script file.

% Some initial conditions.
h0 = 0.0; % Assumes sea level. Adjust to correct pressure altitude.
v0 = 0.0; % Assuming start with zero velocity.
m0 = 0.7; % Initial mass of the rocket. Adjust as necessary.
dt = 0.01; % Time step;

% Function definitions.

% Standard quadratic drag.
function f = drag(v)
    % Drag with A = 30cm^2, Cd = 1.
    % Density will be computed separately.
    k = 0.5 * 1 * 0.003;
    f = k*v^2;
end

% Returns constant thrust of 80N for 10s.
% You can replace it with a more realistic thrust profile.
function f = thrust(t)
    if t<=10
        f = 80.0;
    else
        f = 0.0;
    end
end

% Mass as funciton fo time.
% Since thrust is constant above, this is going to be a simple function.
% You might want to adjust it, or even integrate over it for general case.
function m = mass(m0, t)
    % Assuming composite rocket Isp of 2,000m/s.
    isp = 2000.0;

    if t<=10
        m = m0 - 80.0*t/isp;
    else
        m = m0 - 80.0*10.0/isp;
    end
end

% Air density as function of altitude.
function r = density(h)
    % Using standard atmosphere and 20 deg C for density.
    r = 1.2041;
    % Using 8.5km for scale height.
    sh = 8500.0;

    r = r * exp(-h/sh);
end

% Gravity currently returns a constant value.
% Adjust to function of altitude if needed.
function g = gravity(h)
    g = 9.807;
end

printf("Initializing...\n");

% Actual code starts here.

v = v0;
h = h0;
aold = 0.0;
t = 0.0;
heng = 0.0;
veng = 0.0;
teng = 0.0;

while v >=0.0
%while t<10
    a = -gravity(h) + (thrust(t) - density(h)*drag(v))/mass(m0, t);

    % Do a Velocity Verlet step.
    v = v + 0.5*(a + aold)*dt;
    h = h + v*dt + 0.5*a*dt*dt;

    if thrust(t)>0.0
        heng = h;
        veng = v;
        teng = t;
    end

    aold = a;
    t = t+dt;

    %printf("<%f\t%f\t%f\t%f\t%f>\n", t, h, v, a, mass(m0,t));
end

printf("Maximum altitude achieved: %f meters at %f seconds.\n", h, t);
printf("Engine cut out at %f seconds, cruising at %f meters and %f m/s.\n", teng, heng, veng);