program repeat3

repeating ground track - Wagner's method
----------------------------------------

please input the mean orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0

please input the mean orbital inclination (degrees)
(0 <= inclination < 180)
? 108

please input the number of orbits in the repeat cycle
? 271

please input the number of nodal days in the repeat cycle
? 19

program repeat3.bas

repeating ground track - Wagner's method
----------------------------------------

mean semimajor axis 7192.2303 kilometers

mean eccentricity 0.0000000000

mean inclination 108.0000 degrees


orbits to repeat 271

solar days to repeat 19.0548


Keplerian period 101.1708 minutes

nodal period 101.2507 minutes


length of nodal day 1444.1545 minutes

fundamental interval 25.2399 degrees

' program repeat3.bas August 3, 2017

' determine mean semimajor axis required
' for a repeating ground track orbit

' Wagner's algorithm

' MMX/RasPi/DOS versions of MMBASIC

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

option default float

option base 1

' astrodynamic and utility constants

const pi2 = 2.0 * pi

' conversion factor - degrees to radians

const dtr = pi / 180.0

' conversion factor - radians to degrees

const rtd = 180.0 / pi

' Earth gravitational constant (km**3/sec**2)

const xmu = 398600.4415

' equatorial radius of the Earth (kilometers)

const req = 6378.14

' first zonal coefficient (non-dimensional)

const xj2 = 1.08262668355e-3

' Earth rotation rate (radians/second)

const omega = 7.29211564186e-05

print " "
print "program repeat3"
print " "
print "repeating ground track - Wagner's method"
print "----------------------------------------"

' request mean orbital elements

do

print " "
print "please input the mean orbital eccentricity (non-dimensional)"
print "(0 <= eccentricity < 1)"

input ecc

if (ecc >= 0.0 and ecc < 1.0) then exit do

loop

do

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

input xdata

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

loop

xinc = dtr * xdata

do

print " "
print "please input the number of orbits in the repeat cycle"

input xnorbits

if (xnorbits > 0) then exit do

loop

do

print " "
print "please input the number of nodal days in the repeat cycle"

input xndays

if (xndays > 0) then exit do

loop

c20 = -xj2

' calculate initial guess for semimajor axis (kilometers)

a0 = xmu^(1.0 / 3.0) * ((xnorbits / xndays) * omega)^(-2.0 / 3.0)

aold = a0 * (1.0 - c20 * (req / a0)^2 * (4.0 * cos(xinc)^2 - (xnorbits / xndays) * cos(xinc) - 1.0))

do

slr = aold * (1.0 - ecc * ecc)

tmp1 = xmu^(1.0 / 3.0) * (xnorbits * omega / xndays)^(-2.0 / 3.0)

tmp2 = (1.0 - 1.5 * c20 * (req / slr)^2 * (1.0 - 1.5 * sin(xinc)^2))^(2.0 / 3.0)

tmp3 = (1.0 + c20 * (req / slr)^2 * (1.5 * (xnorbits / xndays) * cos(xinc) - 0.75 * (5.0 * cos(xinc)^2 - 1.0)))^(2.0 / 3.0)

anew = tmp1 * tmp2 * tmp3

if (abs(anew - aold) < 1.0e-6) then

' converged to within tolerance

exit do

else

' continue iteration

aold = anew

end if

loop

sma = anew

' Keplerian period (seconds)

tkepler = pi2 * sma * sqr(sma / xmu)

' keplerian mean motion (rad/sec)

xmm = sqr(xmu / (sma * sma * sma))

' orbital semiparameter (km)

slr = sma * (1.0 - ecc * ecc)

b = sqr(1.0 - ecc * ecc)

c = req / slr

d = c * c

e = sin(xinc) * sin(xinc)

' perturbed mean motion (rad/sec)

pmm = xmm * (1.0 + 1.5 * xj2 * d * b * (1.0 - 1.5 * e))

' argument of perigee perturbation (rad/sec)

apdot = 1.5 * xj2 * pmm * d * (2.0 - 2.5 * e)

' raan perturbation (rad/sec)

raandot = -1.5 * xj2 * pmm * d * cos(xinc)

' nodal period - time between nodal crossings

tnode = pi2 / (apdot + pmm)

' delta longitude per nodal period (radians)

dlong = tnode * (omega - raandot)

' length of nodal day (seconds)

tnday = pi2 / (omega - raandot)

' number of solar days to repeat

xnsdays = tnode * xnorbits / 86400.0

' print results

print " "
print "program repeat3.bas"
print " "
print "repeating ground track - Wagner's method"
print "----------------------------------------"
print " "

print "mean semimajor axis ", str$(sma, 0, 4) + " kilometers"
print " "
print "mean eccentricity ", str$(ecc, 0, 10)
print " "
print "mean inclination ", str$(rtd * xinc, 0, 4) + " degrees"
print " "
print " "
print "orbits to repeat ", str$(xnorbits)
print " "
print "solar days to repeat ", str$(xnsdays)
print " "
print " "
print "Keplerian period ", str$(tkepler / 60.0, 0, 4) + " minutes"
print " "
print "nodal period ", str$(tnode / 60.0, 0, 4) + " minutes"
print " "
print " "
print "length of nodal day ", str$(tnday / 60.0, 0, 4) + " minutes"
print " "
print "fundamental interval ", str$(dlong * rtd, 0, 4) + " degrees"
print " "

end

program repeat1 - time to repeat ground track - analytic solution

please input the semimajor axis (kilometers)
(semimajor axis > 0)
? 8000

please input the orbital eccentricity (non-dimensional)
(0 <= eccentricity < 1)
? 0.025

please input the orbital inclination (degrees)
(0 <= inclination <= 180)
? 28.5

please input the argument of perigee (degrees)
(0 <= argument of perigee <= 360)
? 120

please input the closure tolerance (degrees)
(a value between 0.1 and 0.5 degrees is recommended)
? 0.1

program repeat1 - time to repeat ground track - analytic solution

mean semimajor axis 8000.000000 kilometers

mean eccentricity (nd) 0.02500000

mean inclination 28.500000 degrees

mean argument of perigee 120.000000 degrees

number of orbits to repeat 2027

number of solar days to repeat 166.705086

Keplerian period 118.684693 minutes

nodal period 118.428872 minutes

length of nodal day 1420.446851 minutes

fundamental interval 30.014776 degrees

closure tolerance 0.100000 degrees

actual closure 0.049141 degrees

' repeat1.bas August 7, 2017

' time to repeat ground track

' Kozai analytic orbit propagation

' MMBASIC for MMX, RasPi and DOS

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

option default float

option base 1

dim ioev%(6), oev1(7), oev2(7), rs(3), vs(3)

' global constants and variables

dim mm, pmm, apdot, raandot, sinc, cinc, slr

' zonal gravity constant (nd)

const xj2 = 0.00108263

' Earth equatorial radius (kilometers)

const req = 6378.14

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

const xmu = 398600.4415

' Earth inertial rotation rate (rad/sec)

const omega = 7.2921151467e-5

const pi2 = 2.0 * pi

const dtr = pi / 180.0

const rtd = 180.0 / pi

print " "
print "program repeat1 - time to repeat ground track - analytic solution"

'''''''''''''''''''''''''''''''''''''''''
' request mean classical orbital elements
'''''''''''''''''''''''''''''''''''''''''

ioev%(1) = 1
ioev%(2) = 1
ioev%(3) = 1
ioev%(4) = 1
ioev%(5) = 0
ioev%(6) = 0

getoe(ioev%(), oev1())

do

print " "

print "please input the closure tolerance (degrees)"

print "(a value between 0.1 and 0.5 degrees is recommended)"

input tol

if (tol > 0.0) then exit do

loop

' set true anomaly to ascending node

oev1(6) = -oev1(4)

' calculate initial mean anomaly (radians)

a = sqr(1.0 - oev1(2) * oev1(2)) * sin(oev1(6))

b = cos(oev1(6)) + oev1(2)

eanom = atan3(a, b)

oev1(6) = modulo(eanom - oev1(2) * sin(eanom))

' compute Kozai mean motion and perturbations

t = 0.0

kozai1(1, t, oev1(), oev2(), rs(), vs())

' Keplerian period (seconds)

tkepler = pi2 / mm

' nodal period => time between nodal crossings

tnode = pi2 / (apdot + pmm)

' delta longitude per nodal period (radians)

dlong = tnode * (omega - raandot)

' length of nodal day (seconds)

tnday = pi2 / (omega - raandot)

' initialize

xlong = 0.0

norbits = 0

' increment nodal crossing and look for closure

do

' increment current longitude

xlong = xlong + dlong

' modulo longitude

if (xlong >= pi2) then xlong = xlong - pi2

' increment number of orbits

norbits = norbits + 1

' check for ground track closure

if ((abs(xlong - pi2) <= dtr * tol) or (abs(xlong) <= dtr * tol)) then exit do

loop

' number of days to repeat ground track

ndays = tnode * norbits / 86400.0

' actual closure delta longitude (radians)

delta = abs(xlong - pi2)

if (abs(xlong) < delta) then delta = abs(xlong)

' print results

print " "
print "program repeat1 - time to repeat ground track - analytic solution"
print " "

print "mean semimajor axis ", str$(oev1(1), 0, 6) + " kilometers"
print " "
print "mean eccentricity (nd) ", str$(oev1(2), 0, 8)
print " "
print "mean inclination ", str$(rtd * oev1(3), 0, 6) + " degrees"
print " "
print "mean argument of perigee ", str$(rtd * oev1(4), 0, 6) + " degrees"
print " "
print "number of orbits to repeat ", norbits
print " "
print "number of solar days to repeat ", str$(ndays, 0, 6)
print " "
print "Keplerian period ", str$(tkepler / 60.0, 0, 6) + " minutes"
print " "
print "nodal period ", str$(tnode / 60.0, 0, 6) + " minutes"
print " "
print "length of nodal day ", str$(tnday / 60.0, 0, 6) + " minutes"
print " "
print "fundamental interval ", str$(dlong * rtd, 0, 6) + " degrees"
print " "
print "closure tolerance ", str$(tol, 0, 6) + " degrees"
print " "
print "actual closure ", str$(delta * rtd, 0, 6) + " degrees"
print " "

end

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

sub kozai1(iniz%, t, oev1(), oev2(), r(), v())

' analytic orbit propagation

' Kozai's method - ECI version

' input

' iniz = initialization flag (1 = initialize, 0 = bypass)
' t = elapsed simulation time (seconds)

' initial orbital elements

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

' output

' final orbital elements and state vector at time = t

' oev2(1) = semimajor axis (kilometers)
' oev2(2) = orbital eccentricity (non-dimensional)
' (0 <= eccentricity < 1)
' oev2(3) = orbital inclination (radians)
' (0 <= oev2(3) <= pi)
' oev2(4) = updated argument of perigee (radians)
' (0 <= oev2(4) <= 2 pi)
' oev2(5) = updated right ascension of ascending node (radians)
' (0 <= oev2(5) <= 2 pi)
' oev2(6) = updated mean anomaly (radians)
' (0 <= oev2(6) <= 2 pi)
' oev2(7) = updated true anomaly (radians)
' (0 <= oev2(6) <= 2 pi)
' r = eci position vector (kilometers)
' v = eci velocity vector (kilometers/second)

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

local b, c, d, e, c1, c2, c3, rmag

local argperi, raani, manomi

' set final orbital elements equal to initial

for i% =1 to 6

oev2(i%) = oev1(i%)

next i%

if (iniz% = 1) then

sinc = sin(oev1(3))
cinc = cos(oev1(3))

' keplerian mean motion

mm = sqr(xmu / (oev1(1) * oev1(1) * oev1(1)))

' orbital semiparameter

slr = oev1(1) * (1.0 - oev1(2) * oev1(2))

b = sqr(1.0 - oev1(2) * oev1(2))
c = req / slr
d = c * c
e = sinc * sinc

' perturbed mean motion

pmm = mm * (1.0 + 1.5 * xj2 * d * b * (1.0 - 1.5 * e))

' argument of perigee perturbation

apdot = 1.5 * xj2 * pmm * d * (2.0 - 2.5 * e)

' raan perturbation

raandot = -1.5 * xj2 * pmm * d * cinc

end if

' set initial values of argument of perigee,
' raan and mean anomaly

argperi = oev1(4)
raani = oev1(5)
manomi = oev1(6)

' propagate orbit

argper = modulo(argperi + apdot * t)

raan = modulo(raani + raandot * t)

manom = modulo(manomi + pmm * t)

oev2(4) = argper
oev2(5) = raan
oev2(6) = manom

' solve Kepler's equation using Danby's method

ecc = oev1(2)

kepler1(manom, ecc, ea, tanom)

oev2(7) = tanom

' position magnitude

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

sargper = sin(argper)
cargper = cos(argper)

' argument of latitude

arglat = modulo(argper + tanom)

sarglat = sin(arglat)
carglat = cos(arglat)

sraan = sin(raan)
craan = cos(raan)

' eci position vector

r(1) = rmag * (craan * carglat - cinc * sarglat * sraan)
r(2) = rmag * (sraan * carglat + cinc * sarglat * craan)
r(3) = rmag * sinc * sarglat

' eci velocity vector

c1 = sqr(xmu / slr)
c2 = ecc * cargper + carglat
c3 = ecc * sargper + sarglat

v(1) = -c1 * (craan * c3 + sraan * cinc * c2)
v(2) = -c1 * (sraan * c3 - craan * cinc * c2)
v(3) = c1 * c2 * sinc

end sub

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

sub kepler1 (manom, ecc, eanom, tanom, niter%)

' solve Kepler's equation for circular,
' elliptic and hyperbolic orbits

' Danby's method

' input

' manom = mean anomaly (radians)
' ecc = orbital eccentricity (non-dimensional)

' output

' eanom = eccentric anomaly (radians)
' tanom = true anomaly (radians)
' niter% = number of algorithm iterations

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

local ktol, xma, s, c, f, fp, fpp, fppp

local delta, deltastar, deltak, sta, cta

' define convergence criterion

ktol = 1.0e-10

xma = manom - pi2 * fix(manom / pi2)

' initial guess

if (ecc = 0.0) then

' circular orbit

tanom = xma

eanom = xma

return

elseif (ecc < 1.0) then

' elliptic orbit

eanom = xma + 0.85 * sgn(sin(xma)) * ecc

else

' hyperbolic orbit

eanom = log(2.0 * xma / ecc + 1.8)

end if

' perform iterations

niter% = 0

do

if (ecc < 1.0) then

' elliptic orbit

s = ecc * sin(eanom)

c = ecc * cos(eanom)

f = eanom - s - xma

fp = 1.0 - c

fpp = s

fppp = c

else

' hyperbolic orbit

s = ecc * sinh(eanom)

c = ecc * cosh(eanom)

f = s - eanom - xma

fp = c - 1.0

fpp = s

fppp = c

end if

niter% = niter% + 1

' check for convergence

if (abs(f) <= ktol or niter% > 20) then exit do

' update eccentric anomaly

delta = -f / fp

deltastar = -f / (fp + 0.5 * delta * fpp)

deltak = -f / (fp + 0.5 * deltastar * fpp + deltastar * deltastar * fppp / 6.0)

eanom = eanom + deltak

loop

if (niter% > 20) then

print " "

print "more than 20 iterations in kepler1"

stop

end if

' compute true anomaly

if (ecc < 1.0) then

' elliptic orbit

sta = sqr(1.0 - ecc * ecc) * sin(eanom)

cta = cos(eanom) - ecc

else

' hyperbolic orbit

sta = sqr(ecc * ecc - 1.0) * sinh(eanom)

cta = ecc - cosh(eanom)

end if

tanom = atan3(sta, cta)

end sub

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

sub getoe(ioev%(), oev())

' interactive request of classical orbital elements

' input

' ioev%(1) = request semimajor axis (1 = yes, 0 = no)
' ioev%(2) = request orbital eccentricity (1 = yes, 0 = no)
' ioev%(3) = request orbital inclination (1 = yes, 0 = no)
' ioev%(4) = request argument of perigee (1 = yes, 0 = no)
' ioev%(5) = request right ascension of the ascending node (1 = yes, 0 = no)
' ioev%(6) = request true anomaly (1 = yes, 0 = no)

' output

' oev(1) = semimajor axis (kilometers)
' oev(2) = orbital eccentricity
' oev(3) = orbital inclination (radians)
' oev(4) = argument of perigee (radians)
' oev(5) = right ascension of the ascending node (radians)
' oev(6) = true anomaly (radians)

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

local sma, ecc, orbinc, argper, raan, tanom

' initialize all orbital elements

for i% = 1 to 6

oev(i%) = 0.0

next i%

if (ioev%(1) = 1) then

' semimajor axis (kilometers)
' (semimajor axis > req)

do

print " "
print "please input the semimajor axis (kilometers)"
print "(semimajor axis > 0)"

input sma

if (sma > 0.0) then exit do

loop

oev(1) = sma

end if

if (ioev%(2) = 1) then

' orbital eccentricity (non-dimensional)
' (0 <= eccentricity < 1)

do
print " "
print "please input the orbital eccentricity (non-dimensional)"
print "(0 <= eccentricity < 1)"

input ecc

if (ecc >= 0.0 and ecc < 1.0) then exit do

loop

oev(2) = ecc

end if

if (ioev%(3) = 1) then

' orbital inclination (degrees)
' (0 <= inclination <= 180)

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

input orbinc

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

loop

oev(3) = dtr * orbinc

end if

if (ioev%(4) = 1) then

' argument of perigee (degrees)
' (0 <= argument of perigee <= 360)

if (ecc > 0.0) then

do
print " "
print "please input the argument of perigee (degrees)"
print "(0 <= argument of perigee <= 360)"

input argper

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

loop

else

argper = 0.0

end if

oev(4) = dtr * argper

end if

if (ioev%(5) = 1) then

' raan (degrees)
' (0 <= raan <= 360)

do
print " "
print "please input the right ascension of the ascending node (degrees"
print "(0 <= raan <= 360)"

input raan

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

loop

oev(5) = dtr * raan

end if

if (ioev%(6) = 1) then

' true anomaly (degrees)
' (0 <= true anomaly <= 360)

do

print " "
print "please input the true anomaly (degrees)"
print "(0 <= true anomaly <= 360)"

input tanom

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

loop

oev(6) = dtr * tanom

end if

end sub

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

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 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