' program rocket1.bas March 19, 2017

' MicroMite eXtreme version

' determines flight performance
' of single stage model rockets

' altitude at burnout, meters
' velocity at burnout, meters per second
' coast time and total flight time, seconds
' maximum altitude, meters

'''''''''''''''''''''''''''

option default float

dim altsite, tempsite, tduration, impulse, massi, mprop, fd, cd

' surface gravity (meters/second^2)

const gravity = 9.80665

' sea level atmospheric density (kilograms/cubic meter)

const rhosl = 1.22557

print " "
print "program rocket1"
print "---------------"

do

print " "

print "please input the launch site altitude (meters)"

print "(sites above sea level are positive, sites below sea level are negative)"

input altsite

loop until (abs(altsite) >= 0.0)

print " "

print "please input the launch site temperature (degrees fahrenheit)"

input tempsite

do
print " "

print "please input the thrust duration (seconds)"

input tduration

loop until (tduration > 0.0)

do

print " "

print "please input the total impulse (newton-seconds)"

input impulse

loop until (impulse > 0.0)

do

print " "

print "please input the initial mass (grams)"

input massi

loop until (massi > 0.0)

do

print " "

print "please input the propellant mass (grams)"

input mprop

loop until (mprop > 0.0)

do

print " "

print "please input the frontal diameter (millimeters)"

input fd

loop until (fd > 0.0)

do

print " "

print "please input the drag coefficient"

input cd

loop until (cd > 0.0)

' convert mass to kilograms and diameter to square meters

massi = 0.001 * massi

mprop = 0.001 * mprop

fd = pi * fd * fd / 4000000.0

' compensate for launch site altitude and temperature

rho = rhosl * density(altsite) / (1.0 + (tempsite - 59.0) / 518.67)

' determine altitude performance

thrust = impulse / tduration

mass = (massi - 0.5 * mprop)

k2 = 0.5 * rho * fd * cd

weight = mass * gravity

b = tduration * sqr(k2 * (thrust - weight)) / mass

c = exp(b)

d = exp(-b)

e = 0.5 * (c + d)

f = (c - d) / (c + d)

' calculate burnout conditions

altbo = (mass / k2) * log(e)

velbo = f * sqr((thrust - weight) / k2)

massbo = massi - mprop

weightbo = massbo * gravity

' compute coast conditions

tcoast = sqr(massbo / (k2 * gravity)) * atn(velbo * sqr(k2 / weightbo))

altcoast = (massbo / (2.0 * k2)) * log(k2 * velbo * velbo / weightbo + 1.0)

' calculate total flight time and maximum altitude

tflight = tduration + tcoast

altmax = altbo + altcoast

' print results

print " "
print "program rocket1"
print "---------------"

print
print "burnout altitude ", str$(altbo, 0, 4), " meters"

print " "
print "burnout velocity ", STR$(velbo, 0, 4), " meters per second"

print " "
print "burnout mass ", STR$(1000.0 * massbo, 0, 4), " grams"

print " "
print "coast time ", str$(tcoast, 0, 4), " seconds"

print " "
print "total flight time ", str$(tflight, 0, 4), " seconds"

print " "
print "maximum altitude ", str$(altmax, 0, 4), " meters"

print " "
print "average thrust ", str$(thrust, 0, 4), " newtons"
print " "

end

''''''''''''''''''''''''''
''''''''''''''''''''''''''

function density(x) as float

' atmospheric density function

' input

' x = altitude (meters)

' output

' density = atmospheric density (kilograms/cubic meter)

''''''''''''''''''''''''''''''''''''''''''''''''''''''''

density = (1.0 - 0.000022556913 * x) ^ 4.256116

end function

' program rocket2.bas April 2, 2017

' determines flight performance of model rockets

' altitude at burnout (meters)
' velocity at burnout (meters per second)
' coast time and total flight time (seconds)
' maximum altitude (meters)

' numerical integration method

' MicroMite eXtreme version

'''''''''''''''''''''''''''

option default float

option base 1

const gravity = 9.80665

const rhosl = 1.22557

const dtr = pi / 180.0

const rtd = 180.0 / pi

' integrator coefficients

const a0 = 0.045
const b0 = 0.3
const c0 = 13.0 / 126.0
const d0 = 5.0 / 18.0
const e0 = 5.0 / 42.0
const f0 = 7.0 / 600.0
const g0 = 7.0 / 30.0
const h0 = 7.0 / 15.0
const i0 = 7.0 / 6.0
const j0 = 25.0 / 63.0
const k0 = 0.7
const l0 = 19.0 / 78.0
const m0 = 35.0 / 312.0
const n0 = 15.0 / 104.0
const o0 = 64.0 / 39.0
const p0 = 70.0 / 39.0
const q0 = 15.0 / 13.0

' global arrays and variables

dim d1(2), d2(2), d3(2), d4(2)

dim acc(2), xp(2), xs(2), vp(2), vs(2)

dim langle, fltflg%, f5, maxthr, susthr, vex, vex2

dim massi, mprop, tt, tp, dt, tsus, tmax, tdur, t6, t7

dim thrust, mass, drag, ffactor

' begin simulation

print " "
print "program rocket2"
print "---------------"

do
print " "
print "please input the launch site altitude (meters)"
input altsite
loop until (altsite >= 0.0)

do
print " "
print "please input the launch site temperature (degrees f)"
input tempsite
loop until (tempsite >= 0.0)

do
print " "
print "please input the maximum thrust (newtons)"
input maxthr
loop until (maxthr > 0.0)

do
print " "
print "please input the sustainer thrust (newtons)"
input susthr
loop until (susthr > 0.0)

do
print " "
print "please input the time of maximum thrust (seconds)"
input tmax
loop until (tmax > 0.0)

do
print " "
print "please input the time of sustainer thrust (seconds)"
input tsus
loop until (tsus > 0.0)

do
print " "
print "please input the thrust duration (seconds)"
input tdur
loop until (tdur > 0.0)

do
print " "
print "please input the initial mass (grams)"
input massi
loop until (massi > 0.0)

massi = 0.001 * massi

do
print " "
print "please input the propellant mass (grams)"
input mprop
loop until (mprop > 0.0)

mprop = 0.001 * mprop

do
print " "
print "please input the frontal diameter (millimeters)"
input fdia
loop until (fdia > 0.0)

fdia = pi * fdia * fdia / 4000000.0

do
print " "
print "please input the drag coefficient (non-dimensional)"
input cd
loop until (cd > 0.0)

do
print " "
print "please input the launch angle"
print "(degrees [ 0 to 90 ] - measured from the ground)"
print "(for example, a vertical launch is an angle of 90 degrees)"
input langle
loop until (langle >= 0.0 and langle <= 90.0)

langle = langle * dtr

do
print " "
print "please input the boost step size (seconds)"
print "(a value of 0.01 seconds is recommended)"
input dtboost
loop until (dtboost > 0.0)

do
print " "
print "please input the coast step size (seconds)"
print "(a value of 0.1 seconds is recommended)"
input dtcoast
loop until (dtcoast > 0.0)

print " "
print "integrating equations of motion ..."
print " "

' compensate for launch site altitude and temperature

rhosite = rhosl * fna(altsite) / (1.0 + (tempsite - 59.0) / 518.67)

f5 = maxthr - susthr

ffactor = 0.5 * rhosite * fdia * cd

t6 = tsus - tmax

' calculate exhaust velocity

vex = (maxthr * 0.5 * tsus + susthr * (tdur - 0.5 * tmax - 0.5 * tsus)) / mprop

vex2 = 2.0 * vex

tp = 0.0

xp(1) = 0.0

xp(2) = 0.0

vp(1) = 0.0

vp(2) = 0.0

fltflg% = 0

' powered flight

np% = tdur / dtboost

dt = dtboost

for ii% = 1 to np%

integrate

tt = tp

plmass

plthrust

' un-comment the following line to display the
' flight conditions during the thrusting flight

' pdata

next ii%

' burnout conditions

xbo = xp(2)

vbo = sqr(vp(1) ^ 2 + vp(2) ^ 2)

' coasting flight

dt = dtcoast

fltflg% = 1

do

if (vp(2) > 0.0) then

integrate

tt = tp

plmass

plthrust

' un-comment the following line to display the
' flight conditions during the coasting flight

' pdata

else

exit do

end if

loop

tcoast = tp - tdur

' print final conditions

print " "
print "program rocket2"
print "---------------"
print " "
print "burnout altitude ", str$(xbo, 0, 4), " meters"
print " "
print "burnout velocity ", str$(vbo, 0, 4), " meters per second"
print " "
print " "
print "coast time ", str$(tcoast, 0, 4), " seconds"
print " "
print "total flight time ", str$(tp, 0, 4), " seconds"
print " "
print " "
print "maximum altitude ", str$(xp(2), 0, 4), " meters"
print " "

end

''''''''''''''''''''''
''''''''''''''''''''''

function fna(x) as float

' atmospheric density function

''''''''''''''''''''''''''''''

fna = (1.0 - 0.000022556913 * x) ^ 4.256116

end function

''''''''
''''''''

sub diffeq

' two-dimensional equations of motion subroutine

''''''''''''''''''''''''''''''''''''''''''''''''

vmag = sqr(vs(1) ^ 2 + vs(2) ^ 2)

if (xp(2) < 1.0) then

uvx = cos(langle)

uvy = sin(langle)

else

uvx = vs(1) / vmag

uvy = vs(2) / vmag

end if

drag = vmag ^ 2 * ffactor * fna(xs(2))

mass = massi - mprop

thrust = 0.0

if (fltflg% = 0) then

t7 = tt - tmax

plmass

plthrust

end if

tacc = (thrust - drag) / mass

acc(1) = tacc * uvx

acc(2) = tacc * uvy - gravity

end sub

'''''''''''
'''''''''''

sub integrate

' numerical integration subroutine

''''''''''''''''''''''''''''''''''

tt = tp

for i% = 1 to 2

xs(i%) = xp(i%)

vs(i%) = vp(i%)

next i%

diffeq

tt = tp + b0 * dt

for i% = 1 to 2

d1(i%) = dt * acc(i%)

xs(i%) = xp(i%) + dt * (b0 * vp(i%) + a0 * d1(i%))

vs(i%) = vp(i%) + b0 * d1(i%)

next i%

diffeq

tt = tp + k0 * dt

for i% = 1 to 2

d2(i%) = dt * acc(i%)

xs(i%) = xp(i%) + dt * (k0 * vp(i%) + f0 * d1(i%) + g0 * d2(i%))

vs(i%) = vp(i%) - h0 * d1(i%) + i0 * d2(i%)

next i%

diffeq

tt = tp + dt

for i% = 1 to 2

d3(i%) = dt * acc(i%)

xs(i%) = xp(i%) + dt * (vp(i%) + l0 * d1(i%) + m0 * d2(i%) + n0 * d3(i%))

vs(i%) = vp(i%) + o0 * d1(i%) - p0 * d2(i%) + q0 * d3(i%)

next i%

diffeq

tp = tp + dt

for i% = 1 to 2

d4(i%) = dt * acc(i%)

xp(i%) = xp(i%) + dt * (vp(i%) + c0 * d1(i%) + d0 * d2(i%) + e0 * d3(i%))

vp(i%) = vp(i%) + c0 * (d1(i%) + d4(i%)) + j0 * (d2(i%) + d3(i%))

next i%

end sub

'''''''
'''''''

sub pdata

' print flight conditions subroutine

''''''''''''''''''''''''''''''''''''

vmag = sqr(vp(1) ^ 2 + vp(2) ^ 2)

rmag = sqr(xp(1) ^ 2 + xp(2) ^ 2)

fpa = rtd * atn(vp(2) / vp(1))

accm = sgn(acc(2)) * sqr(acc(1) ^ 2 + acc(2) ^ 2)

print " "
print "flight time (seconds)"; tp
print " "
print " "
print "altitude (meters)"; xp(2)
print " "
print "horizontal range (meters)"; xp(1)
print " "
print "velocity (meters per second)"; vmag
print " "
print "flight path angle (degrees)"; fpa
print " "
print "acceleration (meters/second/second)"; accm
print " "
print " "
print "mass (grams)"; 1000.0 * mass
print " "
print "thrust (newtons)"; thrust
print " "
print "drag (newtons)"; drag
print " "

end sub

''''''''
''''''''

sub plmass

' piecewise-linear mass subroutine

''''''''''''''''''''''''''''''''''

if (tp <= tmax) then

mass = massi - maxthr * tt ^ 2 / (vex2 * tmax)

exit sub

elseif (tp > tmax and tp <= tsus) then

mass = massi - (maxthr / vex) * (tt - 0.5 * tmax) + f5 * t7 ^ 2 / (vex2 * t6)

exit sub

elseif (tp > tsus and tp < tdur) then

mass = massi - (maxthr * tsus + susthr * t6) / vex2 - (susthr / vex) * (tt - tsus)

exit sub

elseif (tp >= tdur) then

mass = massi - mprop

exit sub

end if

end sub

''''''''''
''''''''''

sub plthrust

' piecewise-linear thrust subroutine

''''''''''''''''''''''''''''''''''''

if (tp <= tmax) then

thrust = maxthr * tt / tmax

exit sub

elseif (tp > tmax and tp <= tsus) then

thrust = maxthr - f5 * (tt - tmax) / t6

exit sub

elseif (tp > tsus and tp < tdur) then

thrust = susthr

exit sub

elseif (tp >= tdur) then

thrust = 0.0

exit sub

end if

end sub

' program rocket3.bas April 3, 2017

' for a given model rocket engine and aerodynamic characteristics,
' this program determines the optimal launch mass of a single-stage
' model rocket which maximizes total altitude (Bengen's maxima)

' MicroMite eXtreme version

'''''''''''''''''''''''''''

option default float

option base 1

' global constants

const gravity = 9.80665

const rhosl = 1.22557

' global variables

dim rho, tduration, impulse, mprop, fd, cd

dim altbo, velbo, massbo, weightbo

dim tcoast, altcoast, tflight

' begin simulation

print " "
print "program rocket3"
print "---------------"

do
print " "
print "launch site altitude (meters)"
input altsite
loop until (altsite >= 0.0)

do
print " "
print "launch site temperature (degrees fahrenheit)"
input tempsite
loop until (tempsite >= 0.0)

do
print
print "thrust duration (seconds)"
input tduration
loop until (tduration > 0.0)

do
print
print "total impulse (newton-seconds)"
input impulse
loop until (impulse > 0.0)

do
print
print "propellant mass (grams)"
input mprop
loop until (mprop > 0.0)

do
print
print "frontal diameter (millimeters)"
input fd
loop until (fd > 0.0)

do
print
print "drag coefficient (non-dimensional)"
input cd
loop until (cd > 0.0)

' convert mass to kilograms and diameter to square meters

mprop = 0.001 * mprop

fd = pi * fd ^ 2 / 4000000.0

' compensate for launch site altitude and temperature

rho = rhosl * fna(altsite) / (1.0 + (tempsite - 59.0) / 518.67)

' determine "best" initial launch mass

m1 = 1.001 * mprop

m2 = 10.0 * mprop

tolm = 1.0e-8

minima(m1, m2, tolm, massi, altmax)

' print results

print " "
print "program rocket3"
print "---------------"
print " "
print "initial mass ", str$(1000.0 * massi, 0, 4), " grams"
print " "
print "burnout altitude ", str$(altbo, 0, 4), " meters"
print " "
print "burnout velocity ", str$(velbo, 0, 4), " meters per second"
print " "
print "burnout mass ", str$(1000.0 * (massi - mprop), 0, 4), " grams"
print " "
print " "
print "coast time ", str$(tcoast, 0, 4), " seconds"
print " "
print "total flight time ", str$(tflight, 0, 4), " seconds"
print " "
print " "
print "maximum altitude ", str$(-altmax, 0, 4), " meters"
print " "

end

''''''''''''''''''''''
''''''''''''''''''''''

function fna(x) as float

' atmospheric density function

''''''''''''''''''''''''''''''

fna = (1.0 - 0.000022556913 * x) ^ 4.256116

end function

'''''''''''''''''''''''''
'''''''''''''''''''''''''

sub alt_func(massi, altmax)

' calculate maximum altitude subroutine

'''''''''''''''''''''''''''''''''''''''

thrust = impulse / tduration

mass = massi - 0.5 * mprop

k2 = 0.5 * rho * fd * cd

weight = mass * gravity

b = tduration * sqr(k2 * (thrust - weight)) / mass

c = exp(b)

d = exp(-b)

e = 0.5 * (c + d)

f = (c - d) / (c + d)

' burnout conditions

altbo = (mass / k2) * log(e)

velbo = f * sqr((thrust - weight) / k2)

massbo = massi - mprop

weightbo = massbo * gravity

' coast conditions

tcoast = sqr(massbo / (k2 * gravity)) * atn(velbo * sqr(k2 / weightbo))

altcoast = (massbo / (2.0 * k2)) * log(k2 * velbo * velbo / weightbo + 1.0)

' total flight time and maximum altitude

tflight = tduration + tcoast

altmax = -(altbo + altcoast)

end sub

''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''

sub minima(a, b, tolm, xmin, fmin)

' one-dimensional minimization

' Brent's method

' input

' a = initial x search value
' b = final x search value
' tolm = convergence criterion

' output

' xmin = minimum x value

' note

' user-defined objective subroutine
' coded as alt_func(x, fx)

' remember: a maximum is simply a minimum
' with a negative attitude!

'''''''''''''''''''''''''''''''''''''

' machine epsilon

LOCAL epsm = 2.23e-16

' golden number

LOCAL c = 0.38196601125

LOCAL d, e, t2, p, q

local r, u, fu

LOCAL x, xm, w

local v, fx, fw

LOCAL tol1, fv

x = a + c * (b - a)

w = x

v = w

e = 0.0
p = 0.0
q = 0.0
r = 0.0

alt_func(x, fx)

fw = fx

fv = fw

for iter% = 1 to 100

if (iter% > 50) then

print ("error in subroutine minima!")
print ("(more than 50 iterations)")

end if

xm = 0.5 * (a + b)

tol1 = tolm * abs(x) + epsm

t2 = 2.0 * tol1

if (abs(x - xm) <= (t2 - 0.5 * (b - a))) then

xmin = x

fmin = fx

exit sub

end if

if (abs(e) > tol1) then

r = (x - w) * (fx - fv)

q = (x - v) * (fx - fw)

p = (x - v) * q - (x - w) * r

q = 2.0 * (q - r)

if (q > 0.0) then

p = -p

end if

q = abs(q)

r = e

e = d

end if

if ((abs(p) >= abs(0.5 * q * r)) or (p <= q * (a - x)) or (p >= q * (b - x))) then

if (x >= xm) then

e = a - x

else

e = b - x

end if

d = c * e

else

d = p / q

u = x + d

if ((u - a) < t2) or ((b - u) < t2) then

d = sgn(xm - x) * tol1

end if

end if

if (abs(d) >= tol1) then

u = x + d

else

u = x + sgn(d) * tol1

end if

alt_func(u, fu)

if (fu <= fx) then

if (u >= x) then

a = x

else

b = x

end if

v = w

fv = fw

w = x

fw = fx

x = u

fx = fu

else

if (u < x) then

a = u

else

b = u

end if

if ((fu <= fw) or (w = x)) then

v = w

fv = fw

w = u

fw = fu

elseif ((fu <= fv) Or (v = x) Or (v = w)) then

v = u

fv = fu

end if

end if

next iter%

end sub