program sunsync1

sun-synchronous orbits - j2 solution
------------------------------------

<1> input mean semimajor axis and eccentricity

<2> input mean perigee and apogee altitudes

? 2

please input the mean perigee altitude (kilometers)
? 350

please input the mean apogee altitude (kilometers)
? 1000


program sunsync1

sun-synchronous orbits - j2 solution
------------------------------------

mean perigee altitude 350.0000 kilometers

mean apogee altitude 1000.0000 kilometers

mean semimajor axis 7053.1400 kilometers

mean orbital eccentricity 0.0460787678

mean orbital inclination 98.0571 degrees

number of iterations 3

' program sunsync1.bas August 2, 2017

' mean orbital inclination of sun-synchronous orbits

' Kozai j2 perturbations

' 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

' Earth oblateness coefficient

const xj2 = 0.00108263

' begin main program

print " "
print "program sunsync1"
print " "
print "sun-synchronous orbits - j2 solution"
print "------------------------------------"

do

print " "
print "<1> input mean semimajor axis and eccentricity"
print " "
print "<2> input mean perigee and apogee altitudes"
print " "

input iselect%

if (iselect% = 1 or iselect% = 2) then exit do

loop

if (iselect% = 1) then

print " "
print "please input the mean semimajor axis (kilometers)"
print " "

input sma

print " "
print "please input the mean orbital eccentricity"
print " "

input ecc

' compute perigee and apogee altitudes (kilometers)

hp = sma * (1.0 - ecc) - req

ha = sma * (1.0 + ecc) - req

elseif (iselect% = 2) then

print " "
print "please input the mean perigee altitude (kilometers)"

input hp

print " "
print "please input the mean apogee altitude (kilometers)"

input ha

' compute perigee and apogee geocentric radii

rp = req + hp

ra = req + ha

' compute semimajor axis and eccentricity

sma = 0.5 * (rp + ra)

ecc = (ra - rp) / (ra + rp)

end if

' required nodal regression rate (radians/second)

drdt = dtr * (360.0 / 365.2422) / 86400.0

' Keplerian mean motion (radians/second)

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

' orbital semiparameter (kilometers)

slr = sma * (1.0 - ecc * ecc)

' inclination initial guess

xipar0 = - (2.0 / 3.0) * slr * slr * drdt / (xj2 * req * req * xmm)

' check for possible solution

if (abs(xipar0) > 1.0) then

print " "
print "no sun-synchronous solution !!"
print " "

end

else

' perform iteration

xinc0 = acos(xipar0)

niter% = 0

do

niter% = niter% + 1

' perturbed mean motion (radians/second)

pmm = xmm * (1.0 + 1.5 * xj2 * req * req * sqr(1.0 - ecc * ecc) * (1.0 - 1.5 * sin(xinc0) * sin(xinc0)) / (slr * slr))

' updated inclination

xinc1 = acos(-(2.0 / 3.0) * slr * slr * drdt / (xj2 * req * req * pmm))

' check for convergence or more than 100 iterations

if (abs(xinc1 - xinc0) < 1.0e-9 or niter% > 100) then exit do

' reset inclination

xinc0 = xinc1

loop

end if

' print results

print " "
print " "
print "program sunsync1"
print " "
print "sun-synchronous orbits - j2 solution"
print "------------------------------------"
print " "
print "mean perigee altitude ", str$(hp, 0, 4) + " kilometers"
print " "
print "mean apogee altitude ", str$(ha, 0, 4) + " kilometers"
print " "
print "mean semimajor axis ", str$(sma, 0, 4) + " kilometers"
print " "
print "mean orbital eccentricity ", str$(ecc, 0, 10)
print " "
print "mean orbital inclination ", str$(xinc1 * rtd, 0, 4) + " degrees"
print " "
print "number of iterations ", niter%
print " "

end

program sunsync2

sun-synchronous orbits - j2 + j4 solution
-----------------------------------------

<1> input mean semimajor axis and eccentricity

<2> input mean perigee and apogee altitudes

? 2

please input the mean perigee altitude (kilometers)
? 350

please input the mean apogee altitude (kilometers)
? 1000


program sunsync2

sun-synchronous orbits - j2 + j4 solution
-----------------------------------------

mean perigee altitude 350.0000 kilometers

mean apogee altitude 1000.0000 kilometers

mean semimajor axis 7053.1400 kilometers

mean orbital eccentricity 0.0460787678

mean orbital inclination 98.0306 degrees

' program sunsync2.bas August 6, 2017

' mean orbital inclination of sun-synchronous orbits

' Kozai j2 + j4 iterative solution

' MMX/RasPi/DOS versions of MMBASIC

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

option default float

option base 1

' global constants

dim xmm, drdt0, req2, req4, slr2, slr4, ecc2

' astrodynamic and utility constants

const pi2 = 2.0 * pi

const pidiv2 = 0.5 * 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

' Earth gravity coefficients (non-dimensional)

const xj2 = 0.00108263

const xj4 = 0.0000016109876

const xj22 = 0.0000015747419

' begin main program

print " "
print "program sunsync2"
print " "
print "sun-synchronous orbits - j2 + j4 solution"
print "-----------------------------------------"

do

print " "
print "<1> input mean semimajor axis and eccentricity"
print " "
print "<2> input mean perigee and apogee altitudes"
print " "

input iselect%

if (iselect% = 1 or iselect% = 2) then exit do

loop

if (iselect% = 1) then

print " "
print "please input the mean semimajor axis (kilometers)"
print " "

input sma

print " "
print "please input the mean orbital eccentricity"
print " "

input ecc

' compute perigee and apogee altitudes (kilometers)

hp = sma * (1.0 - ecc) - req

ha = sma * (1.0 + ecc) - req

elseif (iselect% = 2) then

print " "
print "please input the mean perigee altitude (kilometers)"

input hp

print " "
print "please input the mean apogee altitude (kilometers)"

input ha

' compute perigee and apogee geocentric radii

rp = req + hp

ra = req + ha

' compute semimajor axis and eccentricity

sma = 0.5 * (rp + ra)

ecc = (ra - rp) / (ra + rp)

end if

' required nodal regression rate (radians/second)

drdt0 = dtr * (360.0 / 365.2422) / 86400.0

' Keplerian mean motion (radians/second)

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

' orbital semiparameter (kilometers)

slr = sma * (1.0 - ecc * ecc)

' inclination initial guess

xipar0 = - (2.0 / 3.0) * slr * slr * drdt / (xj2 * req * req * xmm)

' check for possible solution

if (abs(xipar0) > 1.0) then

print " "
print "no sun-synchronous solution !!"
print " "

else

' perform iteration

slr2 = slr * slr

slr4 = slr2 * slr2

req2 = req * req

req4 = req2 * req2

ecc2 = ecc * ecc

' compute inclination bounds

xinc1 = xinc0 - 1.0 * dtr

if (xinc1 < pidiv2) then xinc1 = pidiv2

xinc2 = xinc0 + 1.0 * dtr

if (xinc2 > pi) then xinc2 = pi

' solve nonlinear equation

tol = 1.0e-10

root(xinc1, xinc2, tol, xroot, froot)

end if

' print results

print " "
print " "
print "program sunsync2"
print " "
print "sun-synchronous orbits - j2 + j4 solution"
print "-----------------------------------------"
print " "
print "mean perigee altitude ", str$(hp, 0, 4) + " kilometers"
print " "
print "mean apogee altitude ", str$(ha, 0, 4) + " kilometers"
print " "
print "mean semimajor axis ", str$(sma, 0, 4) + " kilometers"
print " "
print "mean orbital eccentricity ", str$(ecc, 0, 10)
print " "
print "mean orbital inclination ", str$(xroot * rtd, 0, 4) + " degrees"
print " "

end

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

sub ss2_func(x, fx)

' delta-raan rate objective function

' input

' x = current orbital inclination (radians)

' output

' fx = objective function

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

local sinc, cinc

local pmm, pmm1, pmm2, pmm3, pmm4, pmm5, pmm6

local drdt, drdt1, drdt2, drdt3, drdt4, drdt5

sinc = sin(x)

cinc = cos(x)

' Kozai j2 + j4 perturbed mean motion (radians/second)

pmm1 = 1.0 + 1.5 * xj2 * req2 * sqr(1.0 - ecc2) * (1.0 - 1.5 * sinc^2) / slr2

pmm2 = (3.0 / 128.0) * xj22 * req4 * sqr(1.0 - ecc2) / slr4

pmm3 = 16.0 * sqr(1.0 - ecc2) + 25.0 * (1.0 - ecc2) - 15.0

pmm4 = (30.0 - 96.0 * sqr(1.0 - ecc2) - 90.0 * (1.0 - ecc2)) * cinc^2

pmm5 = (105.0 + 144.0 * sqr(1.0 - ecc2) + 25.0 * (1.0 - ecc2)) * cinc^4

pmm6 = - (45.0 / 128.0) * xj4 * req4 * sqr(1.0 - ecc2) * ecc2 * (3.0 - 30.0 * cinc^2 + 35.0 * cinc^4) / slr4

pmm = xmm * (pmm1 + pmm2 * (pmm3 + pmm4 + pmm5) + pmm6)

' raan time rate of change

drdt1 = -1.5 * req2 * xj2 * pmm * cinc / slr2

drdt2 = 1.5 * req2 * xj2 / slr2

drdt3 = 1.5 + ecc2 / 6.0 - 2.0 * sqr(1.0 - ecc2)

drdt4 = - ((5.0 / 3.0) - (5.0 / 24.0) * ecc2 - 3.0 * sqr(1.0 - ecc2)) * sinc * sinc

drdt5 = - (35.0 / 8.0) * req4 * xj4 * xmm * (1.0 + 1.5 * ecc2) * ((12.0 - 21.0 * sinc * sinc) / 14.0) * cinc / slr4

drdt = drdt1 * (1.0 + drdt2 * (drdt3 + drdt4)) - drdt5

' calculate the objective function as the difference
' between the desired and current raan rate

fx = drdt - drdt0

end sub

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

sub root(x1, x2, tol, xroot, froot)

' real root of a single non-linear function subroutine

' input

' x1 = lower bound of search interval
' x2 = upper bound of search interval
' tol = convergence criter%ia

' output

' xroot = real root of f(x) = 0
' froot = function value

' note: requires subroutine ss2_func

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

local eps, a, b, c, d, e, fa, fb, fcc, tol1

local xm, p, q, r, s, xmin, tmp

eps = 2.23e-16

e = 0.0

a = x1

b = x2

ss2_func(a, fa)

ss2_func(b, fb)

fcc = fb

for iter% = 1 to 50

if (fb * fcc > 0.0) then

c = a

fcc = fa

d = b - a

e = d

end if

if (abs(fcc) < abs(fb)) then

a = b

b = c

c = a

fa = fb

fb = fcc

fcc = fa

end if

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

xm = 0.5 * (c - b)

if (abs(xm) <= tol1 or fb = 0.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 / fcc

r = fb / fcc

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.0) then q = -q

p = abs(p)

min = abs(e * q)

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

if (min < tmp) then min = tmp

if (2.0 * p < min) 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

ss2_func(b, fb)

next iter%

froot = fb

xroot = b

end sub