Model rocket simulator

# Prevent Octave from thinking that this # is a function file:
1;
clear;

# Constants
    g   = 9.8; #m/s2
    rho = 1.2; #kg/m3

# Motor Data

    # Read thrust curve files
        load A8.m
        load B6.m
        load C6.m
        load D12.m

    # Motor variable format: Propellant mass - Mass before firing - Mass after firing

    # Original Estes Engines
        MotorDataA8   = [ .0033 ; .01635 ; .0102 ];
        MotorDataB6   = [ .0056 ; .01823 ; .0097 ];
        MotorDataC6   = [ .0108 ; .02400 ; .0094 ];
        MotorDataD12  = [ .0211 ; .04260 ; .0160 ];

    # Derived / Clustered
        MotorDataTwoDee   = 2*MotorDataD12;
        MotorDataThreeDee = 3*MotorDataD12;
        MotorDataFourDee  = 4*MotorDataD12;
        TwoDee   = [ D12(:,1) , 2*D12(:,2) ];
        ThreeDee = [ D12(:,1) , 3*D12(:,2) ];
        FourDee  = [ D12(:,1) , 4*D12(:,2) ];

# Rocket data
                                                                        # KURIOSUM: Et hønseæg vejer i gennemsnit 63 g. Kilde: Danæg

    # Rascal
    # m  = 0.0411  kg
    # d  = 0.0250  m
    # d2 = 0.3     m

    # Code Red
    # m  = 0.0666  kg
    # d  = 0.034   m
    # d2 = 0.3     m

    # Modular Model Rocket
    # m  = 0.095 kg
    # d  = 0.044 m
    # d2 = 0.3   m

    m         = 0.060; # Rocket mass in kg without engine
    d         = 0.050; # Rocket frontal diameter in m
    d2        = 0.300; # Frontal diameter after deployment (chute diameter)

    cd        = 0.65; # Rocket drag coefficient
    cd2       = 1.00; # Drag coefficient after deployment

    # Motor data
        Motor     = B6;
        MotorData = MotorDataB6;
        delay     = 4; # Delay charge burn time

# Set simulation timestep
    dt = 0.001;

# Interpolate thrust curve
    tf         = [ 0 : dt : max(Motor(:,1)) ];
    Thrust     = abs(interp1(Motor(:,1),Motor(:,2),tf,"cubic"));
    mDelay     = MotorData(2)-MotorData(3)-MotorData(1); # Mass of delay charge, ejection charge and clay cap
    mDelay     = mDelay/2; # Assuming ejection charge and cap together weighs about the same as the delay charge
    dmDelay    = dt/delay*mDelay;

### SIMULATION ###

m   = m + MotorData(2);
A   = pi*(d /2)^2;
A2  = pi*(d2/2)^2;

runtime = 300; # Maximum simulation runtime in seconds

i = 1;
t = 0;
h = 0;
v = 0;
a = 0;
I = 0;

done      = false;
ignition  = false;
liftoff   = false;
burnout   = false;
apogee    = false;
deploy    = false;
landing   = false;

printf("\n\n*** SIMULATION ***\n\n")

PropMass = MotorData(1)
Burntime = max(Motor(:,1))
TotalImp = sum(Thrust)*dt
SpeciImp = TotalImp/(PropMass*g)

while (!done)

    t = i * dt;

    if( t < Burntime )
        thrust = Thrust(i);
        TotalImpulseSoFar = thrust*dt;
    else
        thrust = 0;
    endif

    drag   = -0.5 * rho * cd * A * v^2 * sign(v) ;
    q      =  0.5 * rho *          v^2; # Dynamic Pressure https://en.wikipedia.org/wiki/Max_Q

    if( ignition && !burnout )
        dm = PropMass*TotalImpulseSoFar/TotalImp;
        m  = m-dm;
    endif

    weight = -m*g;

    if( burnout && !deploy )
        m  = m-dmDelay;
    endif

    F = thrust + drag + weight;

    a    = F/m;

    if( !liftoff && a<0 )
        a = 0;
    endif

    gees = abs(a/g+1);

    dv = a * dt;
    v  = v + dv;

    dh = v * dt;
    h  = h + dh;

    if( h < 0 )
        h=0;
    endif

    if( !ignition && thrust > 0 )
        ignition = true;
        printf("\nIgnition time  ")
        t
    endif

    if( !liftoff && h > 0 )
        liftoff = true;
        printf("\nLiftoff time   ")
        t
    endif

    if( ignition && !burnout && ( thrust == 0 ) )
        burnout = true;
        printf("\nBurnout at     ")
        t
        printf("Burnout height ")
        h
        printf("Burnout speed  ")
        v
    endif

    if( !apogee && liftoff && v<0 )
        apogee = true;
        printf("\nApogee at      ")
        t
        printf("Apogee height  ")
        h
    endif

    # Deploy condition - time or something else?

    if( t == ( Burntime + delay ) )
        deploy = true;
        cd = cd2;
        A  = A2;
        m=m-mDelay;
        printf("\nDeploy at      ")
        t
        printf("Deploy height  ")
        h
        printf("Deploy speed   ")
        v
    endif

    if( liftoff && ( h == 0 ) )
        landing = true;
        printf("\nLanding at     ")
        t
        printf("Landing speed  ")
        v
    endif

    if( deploy && ( h == 0 ) )
        done = true; # Rocket didn't leave the pad
    endif

    if( landing || ( t >= runtime ) )
        done = true;
        if( t>=runtime )
            printf("\nLanding not detected within preset runtime - stopping")
        endif
    endif

    t2(i)      = t;
    thrust2(i) = thrust;
    drag2(i)   = drag;
    F2(i)      = F;
    a2(i)      = a;
    v2(i)      = v;
    h2(i)      = h;
    m2(i)      = m;
    q2(i)      = q;
    gees2(i)   = gees;

    t = t + dt;
    i = i + 1;

endwhile

if( burnout && !liftoff )
    printf("\nRocket didn't leave the pad - too heavy!")
endif

graphics_toolkit ("gnuplot")

# Time Plots

plot(t2,h2);
grid on;
xlabel("time / [s]");
ylabel("altitude / [m]");
print h.png

plot(t2,a2);
grid on;
xlabel("time / [s]");
ylabel("acceleration / [m/s^2]");
print a.png

plot(t2,v2);
grid on;
xlabel("time / [s]");
ylabel("speed / [m/s]");
print v.png

plot(t2,drag2);
grid on;
xlabel("time / [s]");
ylabel("drag / [N]");
print drag.png

plot(t2,m2);
grid on;
xlabel("time / [s]");
ylabel("mass / [kg]");
print m.png

plot(tf,Thrust);
grid on;
xlabel("time / [s]");
ylabel("thrust / [N]");
print thrust.png

plot(t2,q2);
grid on;
xlabel("time / [s]");
ylabel("dynamic pressure / [Pa]");
print dynamicpressure.png

plot(t2,gees2);
grid on;
xlabel("time / [s]");
ylabel("gees / [G]");
print gees.png


printf("\n*** END ***\n\n")