' hohmann.bas March 16, 2017

' Hohmann two impulse orbit transfer between
' planar and non-coplanar circular orbits

' Micromite eXtreme version

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

option default float

dim hn1, hn2, hn3, dinc, rtol, maxiter, niter

dim alt1, alt2, inc1, inc2, rmin, rmax

' astrodynamic and utility constants

const dtr = pi / 180.0

const rtd = 180.0 / pi

' earth gravitational constant (km^3/sec^2)

const mu = 398600.436233

' earth equatorial radius (kilometers)

const req = 6378.1363

' root-finding tolerance

rtol = 1.0e-8

' maximum number of root-finding iterations allowed

maxiter = 100

' request inputs

print " "
print "Hohmann Orbit Transfer Analysis"
print "==============================="
print " "

do

print " "

print "please input the initial altitude (kilometers)"

input alt1

if (alt1 > 0.0) then

exit do

end if

loop

do

print " "

print "please input the final altitude (kilometers)"

input alt2

if (alt2 > 0.0) then

exit do

end if

loop

do

print " "

print "please input the initial orbital inclination (degrees)"
print "(0 <= inclination <= 180))"

input inc1

if (inc1 >= 0.0 and inc1 <= 180.0) then

exit do

end if

loop

do

print " "

print "please input the final orbital inclination (degrees)"
print "(0 <= inclination <= 180)"

input inc2

if (inc2 >= 0.0 and inc2 <= 180.0) then

exit do

end if

loop

' convert orbit inclinations to radians

inc1 = inc1 * dtr

inc2 = inc2 * dtr

''''''''''''''''''''''''''''''''''''
' solve the orbit transfer problem '
''''''''''''''''''''''''''''''''''''

' calculate total inclination change (radians)

dinc = abs(inc2 - inc1)

' compute geocentric radii of initial and final orbits (kilometers)

r1 = req + alt1

r2 = req + alt2

' compute "normalized" radii

hn1 = sqr(2.0 * r2 / (r2 + r1))

hn2 = sqr(r1 / r2)

hn3 = sqr(2.0 * r1 / (r2 + r1))

' compute "local circular velocity" of initial and final orbits (km/sec)

v1 = sqr(mu / r1)

v2 = sqr(mu / r2)

' compute transfer orbit semimajor axis (kilometers)

smat = 0.5 * (r1 + r2)

' compute transfer orbit eccentricity (non-dimensional)

if (r1 > r2) then

rmax = r1

else

rmax = r2

end if

if (r1 < r2) then

rmin = r1

else

rmin = r2

end if

ecct = (rmax - rmin) / (r1 + r2)

' compute transfer orbit perigee and apogee radii and velocities

rp = smat * (1.0 - ecct)

ra = smat * (1.0 + ecct)

vt1 = sqr(2.0 * mu * ra / (rp * (rp + ra)))

vt2 = sqr(2.0 * mu * rp / (ra * (rp + ra)))

' compute transfer orbit period (seconds)

taut = 2.0 * pi * sqr(smat^3 / mu)

tof = 0.5 * taut

if (abs(dinc) = 0.0) then

'''''''''''''''''''''''''
' coplanar orbit transfer
'''''''''''''''''''''''''

if (r2 > r1) then

' higher-to-lower transfer

dv1 = vt1 - v1

dv2 = v2 - vt2

else

''''''''''''''''''''''''''
' lower-to-higher transfer
''''''''''''''''''''''''''

dv1 = v1 - vt2

dv2 = vt1 - v2

end if

dinc1 = 0.0

dinc2 = 0.0

inct = inc1

else

'''''''''''''''''''''''''''''
' non-coplanar orbit transfer
'''''''''''''''''''''''''''''

realroot(0, dinc, rtol, maxiter, niter, xroot, froot)

' calculate delta-v's

dinc1 = xroot

dinc2 = dinc - dinc1

dv1 = v1 * sqr(1.0 + hn1 * hn1 - 2.0 * hn1 * cos(dinc1))

dv2 = v1 * sqr(hn2 * hn2 * hn3 * hn3 + hn2 * hn2 - 2.0 * hn2 * hn2 * hn3 * cos(dinc2))

if (inc2 > inc1) then

inct = inc1 + dinc1

else

inct = inc1 - dinc1

end if

end if

' print results

print " "
print "Hohmann Orbit Transfer Analysis"
print "==============================="
print " "

print "initial orbit altitude ", str$(alt1, 0, 4), " kilometers"
print " "
print "initial orbit radius ", str$(alt1 + req, 0, 4), " kilometers"
print " "
print "initial orbit inclination ", str$(inc1 * rtd, 0, 4), " degrees"
print " "
print "initial orbit velocity ", str$(1000.0 * v1, 0, 4), " meters/second"
print " "
print " "
print "final orbit altitude ", str$(alt2, 0, 4), " kilometers"
print " "
print "final orbit radius ", str$(alt2 + req, 0, 4), " kilometers"
print " "
print "final orbit inclination ", str$(inc2 * rtd, 0, 4), " degrees"
print " "
print "final orbit velocity ", str$(1000.0 * v2, 0, 4), " meters/second"
print " "
print " "
print "first inclination change ", str$(dinc1 * rtd, 0, 4), " degrees"
print " "
print "second inclination change ", str$(dinc2 * rtd, 0, 4), " degrees"
print " "
print "total inclination change ", str$(rtd * (dinc1 + dinc2), 0, 4), " degrees"
print " "
print " "
print "first delta-v ", str$(1000.0 * dv1, 0, 4), " meters/second"
print " "
print "second delta-v ", str$(1000.0 * dv2, 0, 4), " meters/second"
print " "
print "total delta-v ", str$(1000.0 * (dv1 + dv2), 0, 4), " meters/second"
print " "
print " "
print "transfer orbit semimajor axis ", str$(smat, 0, 4), " kilometers"
print " "
print "transfer orbit eccentricity ", str$(ecct, 0, 8)
print " "
print "transfer orbit inclination ", str$(rtd * inct, 0, 4), " degrees"
print " "
print "transfer orbit perigee velocity ", str$(1000.0 * vt1, 0, 4), " meters/second"
print " "
print "transfer orbit apogee velocity ", str$(1000.0 * vt2, 0, 4), " meters/second"
print " "
print " "
print "transfer orbit coast time ", str$(tof, 0, 4), " seconds"
print " ", str$(tof / 60.0, 0, 4), " minutes"
print " ", str$(tof / 3600.0, 0, 4), " hours"

end

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

sub hohmfunc (x, fx)

' inclination objective function

' input

' x = function argument

' output

' fx = function value

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

local hn, a, b

dinc1 = x

hn = hn2 * hn2

a = hn1 * sin(dinc1) / sqr(1.0 + hn1 * hn1 - 2.0 * hn1 * cos(dinc1))

b = hn * hn3 * (sin(dinc) * cos(dinc1) - cos(dinc) * sin(dinc1))

' calculate objective function value

fx = a - b / sqr(hn * hn3 * hn3 + hn - 2.0 * hn * hn3 * cos(dinc - dinc1))

end sub

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

sub realroot(xl, xu, tol, maxiter, niter, xroot, froot)

' real root of a single non-linear function subroutine

' input

' xl = lower bound of search interval
' xu = upper bound of search interval
' tol = convergence criteria
' maxiter = maximum number of iterations allowed

' output

' xroot = real root of f(x) = 0
' froot = function value
' niter = actual number of iterations

' note: requires sub hohmfunc

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

local beps, a, b

local c, d, e

local fa, fb, fcc

local tol1, xm, p

local q, r, s

local xmin, tmp

beps = 2.2e-16

e = 0.0

a = xl

b = xu

hohmfunc(a, fa)

hohmfunc(b, fb)

fc = fb

niter = 0

for iter = 1 to maxiter

' count number of iterations

niter = niter + 1

if (fb * fc > 0.0) then

c = a

fc = fa

d = b - a

e = d

end if

if (abs(fc) < abs(fb)) then

a = b

b = c

c = a

fa = fb

fb = fc

fc = fa

end if

tol1 = 2.0 * beps * abs(b) + 0.5 * tol

xm = 0.5 * (c - b)

if (abs(xm) <= tol1 or fb = 0) then exit for

if (abs(e) >= tol1 and abs(fa) > abs(fb)) then

s = fb / fa

if (a = c) then

p = 2.0 * xm * s

q = 1.0 - s

else

q = fa / fc

r = fb / fc

p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))

q = (q - 1.0) * (r - 1.0) * (s - 1.0)

end if

if (p > 0) then q = -q

p = abs(p)

xmin = abs(e * q)

tmp = 3.0 * xm * q - abs(tol1 * q)

if (xmin < tmp) then xmin = tmp

if (2.0 * p < xmin) then

e = d

d = p / q

else

d = xm

e = d

end if

else

d = xm

e = d

end if

a = b

fa = fb

if (abs(d) > tol1) then

b = b + d

else

b = b + sgn(xm) * tol1

end if

hohmfunc(b, fb)

next iter

froot = fb

xroot = b

end sub