' demo_kepler.bas April 11, 2017

' demonstrates how to interact with several MMBASIC
' subroutines which solve Kepler's equation for
' circular, elliptic, parabolic and hyperbolic orbits

' MicroMite eXtreme version

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

option default float

option base 1

const pi2 = 2.0 * pi, pidiv2 = 0.5 * pi

const dtr = pi / 180.0, rtd = 180.0 / pi

print " "
print "numerical solutions of Kepler's equation"
print "----------------------------------------"

do

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

loop until (ecc >= 0.0)

if (ecc <> 1.0) then

do

print " "
print "please input the mean anomaly in degrees"
print "(0 <= mean anomaly <= 360)"
print " "
input manom

loop until (manom >= 0.0 and manom <= 360.0)

manom = dtr * manom

end if

print " "

if (ecc <> 1.0) then

' circular, elliptic and hyperbolic orbits

print "kepler1 solution"
print "----------------"
print " "

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

kprint(ecc, manom, eanom, tanom, niter%)

print "kepler2 solution"
print "----------------"
print " "

kepler2 (manom, ecc, eanom, tanom, niter%)

kprint(ecc, manom, eanom, tanom, niter%)

if (ecc < 1.0) then

print "kepler3 solution"
print "----------------"
print " "

kepler3 (manom, ecc, eanom, tanom)

kprint(ecc, manom, eanom, tanom, 2)

end if

else

' parabolic orbits

do

print " "
print "please input the time relative to perihelion passage (days)"
print " "
input tpp

loop until (tpp >= 0.0)

do

print " "
print "please input the perihelion radius (AU)"
print " "
input rp

loop until (rp > 0.0)

kepler4(tpp, rp, ecc, rmag, tanom, niter%)

print " "
print "kepler4 solution"
print "----------------"
print " "
print "time wrt perihelion ", str$(tpp, 0, 8), " days"
print " "
print "eccentricity ", str$(ecc, 0, 10)
print " "
print "perihelion radius ", str$(rp, 0, 8), " AU"
print " "
print "heliocentric distance ", str$(rmag, 0, 8), " AU"
print " "
print "true anomaly ", str$(rtd * tanom, 0, 8), " degrees"
print " "
print "iterations ", niter%
print " "

end if

end

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

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 kepler2 (manom, ecc, eanom, tanom, niter%)

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

' Danby's method with Mikkola's initial guess

' 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, den, xma, alpha, beta, z, s, ds

local f, fp, fpp, fppp

local delta, deltastar, deltak, sta, cta

' define convergence criterion

ktol = 1.0e-10

if (ecc = 0.0) then

' circular orbit

tanom = manom

eanom = manom

return

end if

den = 1.0 / (4.0 * ecc + 0.5)

if (ecc < 1.0) then

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

alpha = (1.0 - ecc) * den

else

xma = manom

alpha = (ecc - 1.0) * den

end if

beta = 0.5 * xma * den

z = (sqr(alpha * alpha * alpha + beta * beta) + beta) ^ (1.0/3.0)

s = 2 * beta /(z * z + alpha + alpha * alpha / (z * z))

if (ecc > 1.0) then

' hyperbolic orbit

ds = 0.071 * s ^ 5 / (ecc * (1.0 + 0.45 * s ^ 2) * (1.0 + 4.0 * s ^ 2))

else

' elliptic orbit

ds = -0.078 * s * s * s * s * s / (1.0 + ecc)

end if

s = s + ds

' initial guess

if (ecc > 1.0) then

' hyperbolic orbit

eanom = 3 * log(s + sqr(1 + s ^ 2))

else

' elliptic orbit

eanom = xma + ecc * (3.0 * s - 4.0 * s * s * s)

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) 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 kepler2"

stop

end if

' compute true anomaly

if (ecc < 1.0) then

' elliptic orbit

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

cta = cos(eanom) - ecc

else

' hyperbolic orbit

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

cta = ecc - cosh(eanom)

end if

tanom = atan3(sta, cta)

end sub

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

sub kepler3 (manom, ecc, eanom, tanom)

' solve Kepler's equation for elliptic orbits

' Gooding's two iteration method

' input

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

' output

' eanom = eccentric anomaly (radians)
' tanom = true anomaly (radians)

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

local athird, ahalf, e1, emr, e, ee, sta, cta

local w, f, fd, fdd, fddd, dee

athird = 1.0 / 3.0

ahalf = 1.0 / 2.0

' convert to original algorithm argument

e1 = 1.0 - ecc

if (ecc = 0.0) then

' circular orbit

tanom = manom

eanom = manom

return

end if

' mod mean anomaly

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

if (emr < -pi) then

emr = emr + pi2

end if

if (emr > pi) then

emr = emr - pi2

end if

ee = emr

if (ee < 0.0) then

ee = -ee

elseif (ee = 0) then

eanom = ee + (manom - emr)

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

cta = cos(eanom) - ecc

tanom = atan3(sta, cta)

return

end if

e = 1.0 - e1

' solve cubic equation

dcbsol(e, 2.0 * e1, 3.0 * ee, w)

ee = (ee * ee + (pi - ee) * w) / pi

if (emr < 0.0) then

ee = -ee

end if

for i% = 1 to 2

fdd = e * sin(ee)

fddd = e * cos(ee)

if ((e1 + ee * ee / 6.0) >= 0.25) then

f = (ee - fdd) - emr

fd = 1.0 - fddd

else

emkep(e1, ee, etmp)

f = etmp - emr

fd = e1 + 2.0 * e * sin(ahalf * ee) ^ 2

end if

dee = f * fd / (ahalf * f * fdd - fd * fd)

w = fd + ahalf * dee * (fdd + athird * dee * fddd)

fd = fd + dee * (fdd + ahalf * dee * fddd)

ee = ee - (f - dee * (fd - w)) / fd

next i%

eanom = ee + (manom - emr)

' compute true anomaly

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

cta = cos(eanom) - ecc

tanom = atan3(sta, cta)

end sub

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

sub dcbsol (a, b, c, y)

' solve cubic equation function

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

local bsq, d

if (a = 0.0 and b = 0.0 or c = 0) then

y = 0.0

else

bsq = b * b

d = sqr(a) * abs(c)

d = dcubrt(d + sqr(b * bsq + d * d)) ^ 2

y = 2.0 * c / (d + b + bsq / d)

end if

end sub

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

function dcubrt(x) as float

' cube root function

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

local athird, yy, tmp

athird = 1.0 / 3.0

yy = abs(x)

if (x = 0.0) then

tmp = 0.0

else

tmp = yy ^ athird

tmp = tmp - athird * (tmp - yy / tmp ^ 2)

tmp = sgn(x) * tmp

end if

dcubrt = tmp

end function

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

sub emkep(e1, ee, y)

' solves ee - (1 - e1) * sin(ee)
' when (e1, ee) is close to (1, 0)

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

local ee2, term, d, y0

y = e1 * ee

ee2 = -ee * ee

term = ee * (1.0 - e1)

d = 0.0

do

d = d + 2

term = term * ee2 / (d * (d + 1))

y0 = y

y = y - term

if (y0 = y) then

exit do

end if

loop

end sub

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

sub kepler4(t, q, ecc, r, tanom, niter%)

' Kepler's equation for heliocentric,
' parabolic and near-parabolic orbits

' input

' t = time relative to perihelion passage (days)
' q = perihelion radius (AU)
' ecc = orbital eccentricity (non-dimensional)

' output

' r = heliocentric distance (AU)
' tanom = true anomaly (radians)
' niter% = number of algorithm iterations

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

local e2, fac, tau, e20, a, b, u, u2

local c1, c2, c3, x, y

const kgauss = 0.01720209895

e2 = 0.0

fac = 0.5 * ecc

niter% = 0

tau = kgauss * t

do

niter% = niter% + 1

e20 = e2

a = 1.5 * sqr(fac / (q * q * q)) * tau

b = (sqr(a * a + 1.0) + a)^(1.0/3.0)

u = b - 1.0 / b

u2 = u * u

e2 = u2 * (1.0 - ecc) / fac

stumpff(e2, c1, c2, c3)

fac = 3.0 * ecc * c3

' check for convergence

if (abs(e2 - e20) < 1.0e-9) then exit do

loop

if (niter% > 15) then

print " "

print "more than 15 iterations in kepler4"

stop

end if

' heliocentric distance (AU)

r = q * (1.0 + u2 * c2 * ecc / fac)

' x and y coordinates

x = q * (1.0 - u2 * c2 / fac)

y = q * sqr((1.0 + ecc) / fac) * u * c1

' true anomaly (radians)

tanom = atan3(y, x)

end sub

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

sub stumpff(e2, c1, c2, c3)

' Stumpff functions for argument E^2

' input

' e2 = eccentric anomaly squared

' output

' c1, c2, c3 = Stumpff functions

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

local deltac, n%

c1 = 0.0

c2 = 0.0

c3 = 0.0

deltac = 1.0

n% = 1

do

c1 = c1 + deltac

deltac = deltac / (2.0 * n%)

c2 = c2 + deltac

deltac = deltac / (2.0 * n% + 1.0)

c3 = c3 + deltac

deltac = -e2 * deltac

n% = n% + 1

' check for convergence

if (abs(deltac) < 1.0e-12) then exit do

loop

end sub

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

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

' print results subroutine

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

local kerror

if (ecc < 1.0) then

' circular and elliptic orbits

print "eccentricity = ", str$(ecc, 0, 8)
print " "
print "mean anomaly = ", str$(rtd * manom, 0, 8), " degrees"
print " "
print "eccentric anomaly = ", str$(rtd * eanom, 0, 8), " degrees"
print " "
print "true anomaly = ", str$(rtd * tanom, 0, 8), " degrees"
print " "
print "iterations = ", niter%
print " "

kerror = eanom - ecc * sin(eanom) - manom

print "solution error = ", str$(kerror, 0, 10)
print " "

else

' hyperbolic orbits

print "eccentricity = ", str$(ecc, 0, 8)
print " "
print "hyperbolic anomaly = ", str$(rtd * manom, 0, 8), " degrees"
print " "
print "eccentric anomaly = ", str$(rtd * eanom, 0, 8), " degrees"
print " "
print "true anomaly = ", str$(rtd * tanom, 0, 8), " degrees"
print " "
print "iterations = ", niter%
print " "

kerror = ecc * sinh(eanom) - eanom - manom

print "solution error = ", str$(kerror, 0, 10)
print " "

end if

end sub

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

function cosh(x) as float

' hyperbolic cosine function

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

cosh = 0.5 * (exp(x) + exp(-x))

end function

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

function sinh(x) as float

' hyperbolic sine function

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

sinh = 0.5 * (exp(x) - exp(-x))

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