' hyper1.bas July 20, 2017

' Battin's method for single impulse coplanar hyperbolic
' injection from a circular Earth park orbit

' MMBASIC for MMX, Raspberry Pi and DOS

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

option default float

option base 1

dim pv(3), shat(3), rhat1(3), rhat2(3)

dim rpo1(3), vpo1(3), oev_po1(6)

dim rpo2(3), vpo2(3), oev_po2(6)

dim rhyper1(3), vhyper1(3), oev_hyper1(6)

dim rhyper2(3), vhyper2(3), oev_hyper2(6)

dim deltav1(3), deltav2(3)

' angular conversion factors

const pi2 = 2.0 * pi, rtd = 180.0 / pi, dtr = pi / 180.0

' gravitational constant of the Earth (kilometers^3/second^2)

const mu = 398600.4415

' equatorial radius of the Earth (kilometers)

const req = 6378.14

' request user inputs

print " "
print "Interplanetary Injection from a Circular Earth Park Orbit"
print "---------------------------------------------------------"
print " "

do
print " "
print "please input the altitude of the circular Earth park orbit (kilometers)"
input alt

if (alt > 0.0) then exit do

loop

do
print " "
print "please input the orbital inclination of the Earth park orbit (degrees)"
print "(0 <= inclination <= 180)"

input orbinc

if (orbinc >= 0.0 and orbinc <= 180.0) then exit do

loop

do
print " "
print "please input the C3 of the departure hyperbola (kilometers^2/second^2)"
print "(C3 > 0)"

input c3

if (c3 > 0.0) then exit do

loop

do
print " "
print "please input the right ascension of the outgoing asymptote (degrees)"
print "(0 <= right ascension <= 360)"

input rasc

if (rasc >= 0.0 and rasc <= 360.0) then exit do

loop

do
print " "
print "please input the declination of the outgoing asymptote (degrees)"

input decl

if (abs(decl) <= 90.0) then exit do

loop

if (abs(decl) >= orbinc) then

print " "
print "*********************************************************************"
print "Please note: this program is valid for |DLA| < park orbit inclination"
print "*********************************************************************"

end

end if

' geocentric radius of park orbit (kilometers)

rpmag = req + alt

' convert angular elements to radians

xinc = orbinc * dtr

decl_asy = decl * dtr

rasc_asy = rasc * dtr

' v-infinity (km/sec)

vinf = sqr(c3)

' unit vector along the Earth's spin axis

pv(1) = 0.0
pv(2) = 0.0
pv(3) = 1.0

cdecl_asy = cos(decl_asy)

sdecl_asy = sin(decl_asy)

crasc_asy = cos(rasc_asy)

srasc_asy = sin(rasc_asy)

' compute unit asymptote vector

shat(1) = cdecl_asy * crasc_asy

shat(2) = cdecl_asy * srasc_asy

shat(3) = sdecl_asy

' beta = angle between unit asymptote vector and spin axis

beta = acos(vdot(pv(), shat()))

sum1 = acos(cos(beta) / sin(xinc))

' alpha = angle between unit asymptote vector and position vector

alpha = asin(1.0 / (1.0 + rpmag * vinf * vinf / mu))

''''''''''''''''''''''
' tangential injection
''''''''''''''''''''''

' injection argument of latitude #1

arglat1 = modulo(sum1 - alpha)

' injection argument of latitude #2

arglat2 = modulo(-sum1 - alpha)

' injection raan #1

raan1 = modulo(rasc_asy + pi + asin(cotan(beta) / tan(xinc)))

' injection raan #2

raan2 = modulo(rasc_asy + 2.0 * pi - asin(cotan(beta) / tan(xinc)))

' opportunity #1 - park orbit state vector and orbital elements

oev_po1(1) = rpmag

oev_po1(2) = 0.0

oev_po1(3) = xinc

oev_po1(4) = 0.0

oev_po1(5) = raan1

oev_po1(6) = arglat1

orb2eci(oev_po1(), rpo1(), vpo1())

uvector(rpo1(), rhat1())

ctheta = vdot(shat(), rhat1())

d = sqr(mu / ((1.0 + ctheta) * rpmag) + vinf * vinf / 4.0)

' velocity vector of hyperbola at injection

for i% = 1 to 3

vhyper1(i%) = (d + 0.5 * vinf) * shat(i%) + (d - 0.5 * vinf) * rhat1(i%)

next i%

' injection delta-v vector

for i% = 1 to 3

deltav1(i%) = vhyper1(i%) - vpo1(i%)

rhyper1(i%) = rpo1(i%)

next i%

eci2orb(rhyper1(), vhyper1(), oev_hyper1())

print " "
print "---------------------------------------------------------"
print "Interplanetary Injection from a Circular Earth Park Orbit"
print "---------------------------------------------------------"
print " "

print "departure hyperbola characteristics"
print "-----------------------------------"
print " "

print "c3 ", str$(c3, 0, 4) + " kilometers^2/second^2"
print " "
print "asymptote right ascension ", str$(rtd * rasc_asy, 0, 6) + " degrees"
print " "
print "asymptote declination ", str$(rtd * decl_asy, 0, 6) + " degrees"

print " "
print "orbital elements and state vector of park orbit at injection - opportunity #1"
print "-----------------------------------------------------------------------------"
print " "

oeprint(oev_po1())

svprint(rpo1(), vpo1())

print " "
print "orbital elements and state vector of hyperbola at injection - opportunity #1"
print "----------------------------------------------------------------------------"
print " "

oeprint(oev_hyper1())

svprint(rhyper1(), vhyper1())

print " "
print "injection delta-v vector and magnitude - opportunity #1"
print "-------------------------------------------------------"
print " "

print "x-component of delta-v ", str$(1000.0 * deltav1(1), 0, 6) + " meters/second"
print " "
print "y-component of delta-v ", str$(1000.0 * deltav1(2), 0, 6) + " meters/second"
print " "
print "z-component of delta-v ", str$(1000.0 * deltav1(3), 0, 6) + " meters/second"
print " "
print "delta-v magnitude ", str$(1000.0 * vecmag(deltav1()), 0, 6) + " meters/second"
print " "

' opportunity #2 - park orbit state vector

oev_po2(1) = rpmag

oev_po2(2) = 0.0

oev_po2(3) = xinc

oev_po2(4) = 0.0

oev_po2(5) = raan2

oev_po2(6) = arglat2

orb2eci(oev_po2(), rpo2(), vpo2())

uvector(rpo2(), rhat2())

ctheta = vdot(shat(), rhat2())

d = sqr(mu / ((1.0 + ctheta) * rpmag) + vinf * vinf / 4.0)

' velocity vector of hyperbola at injection

for i% = 1 to 3

vhyper2(i%) = (d + 0.5 * vinf) * shat(i%) + (d - 0.5 * vinf) * rhat2(i%)

next i%

' injection delta-v vector

for i% = 1 to 3

deltav2(i%) = vhyper2(i%) - vpo2(i%)

rhyper2(i%) = rpo2(i%)

next i%

eci2orb(rhyper2(), vhyper2(), oev_hyper2())

print " "
print "orbital elements and state vector of park orbit at injection - opportunity #2"
print "-----------------------------------------------------------------------------"
print " "

oeprint(oev_po2())

svprint(rpo2(), vpo2())

print " "
print "orbital elements and state vector of hyperbola at injection - opportunity #2"
print "----------------------------------------------------------------------------"
print " "

oeprint(oev_hyper2())

svprint(rhyper2(), vhyper2())

print " "
print "injection delta-v vector and magnitude - opportunity #2"
print "-------------------------------------------------------"
print " "

print "x-component of delta-v ", str$(1000.0 * deltav2(1), 0, 6) + " meters/second"
print " "
print "y-component of delta-v ", str$(1000.0 * deltav2(2), 0, 6) + " meters/second"
print " "
print "z-component of delta-v ", str$(1000.0 * deltav2(3), 0, 6) + " meters/second"
print " "
print "delta-v magnitude " str$(1000.0 * vecmag(deltav2()), 0, 6) + " meters/second"
print " "

end

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

sub eci2orb(r(), v(), oev())

' convert eci state vector to classical orbital elements

' input

' r() = eci position vector (kilometers)
' v() = eci velocity vector (kilometers/second)

' output

' oev(1) = semimajor axis (kilometers)
' oev(2) = orbital eccentricity (non-dimensional)
' (0 <= eccentricity < 1)
' oev(3) = orbital inclination (radians)
' (0 <= inclination <= pi)
' oev(4) = argument of perigee (radians)
' (0 <= argument of perigee <= 2 pi)
' oev(5) = right ascension of ascending node (radians)
' (0 <= raan <= 2 pi)
' oev(6) = true anomaly (radians)
' (0 <= true anomaly <= 2 pi)

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

local rmag, vmag, sma, p, q, const1

local h, xk, x1, y1, eccm, inc, xlambdat

local raan, argper, tanom

local rhat(3), vhat(3), hhat(3), vtmp(3), ecc(3)

local hv(3), fhat(3), ghat(3)

' position and velocity magnitude

rmag = vecmag(r())

vmag = vecmag(v())

' position and velocity unit vectors

uvector(r(), rhat())

uvector(v(), vhat())

' angular momentum vectors

vcross(r(), v(), hv())

uvector(hv(), hhat())

' eccentricity vector

for i% = 1 to 3

vtmp(i%) = v(i%) / mu

next i%

vcross(vtmp(), hv(), ecc())

for i% = 1 to 3

ecc(i%) = ecc(i%) - rhat(i%)

next i%

' semimajor axis

sma = 1.0 / (2.0 / rmag - vmag^2 / mu)

' keplerian period

if (sma > 0.0) then

period = pi2 / sqr(mu / sma^3)

else

period = 0.0

end if

p = hhat(1) / (1.0 + hhat(3))

q = -hhat(2) / (1.0 + hhat(3))

const1 = 1.0 / (1.0 + p^2 + q^2)

fhat(1) = const1 * (1.0 - p^2 + q^2)
fhat(2) = const1 * 2.0 * p * q
fhat(3) = -const1 * 2.0 * p

ghat(1) = const1 * 2.0 * p * q
ghat(2) = const1 * (1.0 + p^2 - q^2)
ghat(3) = const1 * 2.0 * q

h = vdot(ecc(), ghat())

xk = vdot(ecc(), fhat())

x1 = vdot(r(), fhat())

y1 = vdot(r(), ghat())

' orbital eccentricity

eccm = sqr(h * h + xk * xk)

' orbital inclination

inc = 2.0 * atn(sqr(p * p + q * q))

' true longitude

xlambdat = atan3(y1, x1)

' check for equatorial orbit

if (inc > 1.0e-10) then

raan = atan3(p, q)

else

raan = 0.0

end if

' check for circular orbit

if (eccm > 1.0e-10) then

argper = modulo(atan3(h, xk) - raan)

else

argper = 0.0

end if

' true anomaly

tanom = modulo(xlambdat - raan - argper)

oev(1) = sma
oev(2) = eccm
oev(3) = inc
oev(4) = argper
oev(5) = raan
oev(6) = tanom

end sub

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

sub orb2eci(oev(), r(), v())

' convert orbital elements to state vector subroutine

' input

' oev(1) = semimajor axis (kilometers)
' oev(2) = orbital eccentricity (non-dimensional)
' oev(3) = orbital inclination (radians)
' (0 <= oev(3) <= pi)
' oev(4) = argument of perigee (radians)
' (0 <= oev(4) <= 2 pi)
' oev(5) = right ascension of ascending node or
' east longitude of ascending node (radians)
' (0 <= oev(5) <= 2 pi)
' oev(6) = true anomaly (radians)
' (0 <= oev(6) <= 2 pi)

' output

' r() = geocentric, equatorial position vector (kilometers)
' v() = geocentric, equatorial velocity vector (kilometers/second)

' special notes

' (1) eci coordinates if oev(5) = raan

' (2) ecf coordinates if oev(5) = east longitude

' (3) valid for elliptic and hyperbolic orbits

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

local sma, ecc, inc, argper, nangle, tanom

local slr, rm, arglat, sarglat, carglat

local sinc, cinc, snangle, cnangle, c1, c2, c3

' "unload" orbital elements array

sma = oev(1)
ecc = oev(2)
inc = oev(3)
argper = oev(4)
nangle = oev(5)
tanom = oev(6)

' semiparameter

slr = sma * (1.0 - ecc * ecc)

' radius magnitude

rm = slr / (1.0 + ecc * cos(tanom))

' argument of latitude

arglat = argper + tanom

sarglat = sin(arglat)

carglat = cos(arglat)

sinc = sin(inc)

cinc = cos(inc)

snangle = sin(nangle)

cnangle = cos(nangle)

c1 = sqr(mu / slr)

c2 = ecc * cos(argper) + carglat

c3 = ecc * sin(argper) + sarglat

' x, y and z components of the position vector

r(1) = rm * (cnangle * carglat - snangle * cinc * sarglat)

r(2) = rm * (snangle * carglat + cinc * sarglat * cnangle)

r(3) = rm * sinc * sarglat

' x, y and z components of the velocity vector

v(1) = -c1 * (cnangle * c3 + snangle * cinc * c2)

v(2) = -c1 * (snangle * c3 - cnangle * cinc * c2)

v(3) = c1 * c2 * sinc

end sub

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

sub svprint(r(), v())

' print position and velocity vectors subroutine

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

local rmag, vmag

print "position vector x-component ", str$(r(1), 0, -12) + " kilometers"
print "position vector y-component ", str$(r(2), 0, -12) + " kilometers"
print "position vector z-component ", str$(r(3), 0, -12) + " kilometers"

rmag = vecmag(r())

print " "

print "position vector magnitude ", str$(rmag, 0, -12) + " kilometers"

print " "

print "velocity vector x-component ", str$(v(1), 0, -12) + " kilometers/second"
print "velocity vector y-component ", str$(v(2), 0, -12) + " kilometers/second"
print "velocity vector z-component ", str$(v(3), 0, -12) + " kilometers/second"

vmag = vecmag(v())

print " "

print "velocity vector magnitude ", str$(vmag, 0, -12) + " kilometers/second"

end sub

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

sub oeprint(oev())

' print classical orbital elements subroutine

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

print "semimajor axis ", str$(oev(1), 0, 4) + " kilometers"
print "eccentricity ", str$(oev(2), 0, 10)
print "orbital inclination ", str$(rtd * oev(3), 0, 4) + " degrees"
print "argument of perigee ", str$(rtd * oev(4), 0, 4) + " degrees"
print "RAAN ", str$(rtd * oev(5), 0, 4) + " degrees"
print "true anomaly ", str$(rtd * oev(6), 0, 4) + " degrees"
print " "

end sub

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

function modulo(x) as float

' modulo 2 pi function

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

local a

a = x - pi2 * fix(x / pi2)

if (a < 0.0) then

a = a + pi2

end if

modulo = a

end function

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

function atan3(a, b) as float

' four quadrant inverse tangent function

' input

' a = sine of angle
' b = cosine of angle

' output

' atan3 = angle (0 =< atan3 <= 2 * pi; radians)

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

local c

if (abs(a) < 1.0e-10) then

atan3 = (1.0 - sgn(b)) * pidiv2

exit function

else

c = (2.0 - sgn(a)) * pidiv2

endif

if (abs(b) < 1.0e-10) then

atan3 = c

exit function

else

atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2)

endif

end function

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

function vecmag(a()) as float

' vector magnitude function

' input

' { a } = column vector ( 3 rows by 1 column )

' output

' vecmag = scalar magnitude of vector { a }

vecmag = sqr(a(1) * a(1) + a(2) * a(2) + a(3) * a(3))

end function

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

sub uvector (a(), b())

' unit vector subroutine

' input

' a = column vector (3 rows by 1 column)

' output

' b = unit vector (3 rows by 1 column)

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

local i%, amag

amag = vecmag(a())

for i% = 1 to 3

if (amag <> 0.0) then

b(i%) = a(i%) / amag

else

b(i%) = 0.0

end if

next i%

end sub

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

sub vcross(a(), b(), c())

' vector cross product subroutine

' { c } = { a } x { b }

' input

' { a } = vector a ( 3 rows by 1 column )
' { b } = vector b ( 3 rows by 1 column )

' output

' { c } = { a } x { b } ( 3 rows by 1 column )

c(1) = a(2) * b(3) - a(3) * b(2)
c(2) = a(3) * b(1) - a(1) * b(3)
c(3) = a(1) * b(2) - a(2) * b(1)

end sub

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

function vdot(a(), b()) as float

' vector dot product function

' c = { a } dot { b }

' input

' n% = number of rows
' { a } = column vector with n rows
' { b } = column vector with n rows

' output

' vdot = dot product of { a } and { b }

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

local i%, c = 0.0

for i% = 1 to 3

c = c + a(i%) * b(i%)

next i%

vdot = c

end function

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

function cotan(x) as float

' cotangent function

' x is input value

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

cotan(x) = 1.0 / tan(x)

end function