Construct a Halo orbit around L1 of the Sun-Mars system

MontBlanc2012

Active member
Joined
Aug 13, 2017
Messages
138
Reaction score
26
Points
28
This is a second post detailing how to construct the state vectors of a prescribed Halo orbit around a Lagrange point. It is a direct continuation of the post Construct a Hal orbit around L2 of the Earth-Moon system so it's probably a good idea to read that post before reading this one. That post references a bunch of earlier posts detailing specific aspects, so if you are really keen you might want to go back and have a look at those too.

Basically, this post illustrates a fairly mechanical process of using some simple tools to construct a Halo orbit around a Lagrange point (in this case L1 of the Sun-Mars) system. Those tools are written using a Lua script, Excel and Python - so working across the three sets of tools feels a little 'clunky'. However, the process of overall can be rationalised so that it all works in one language (e.g., Lua or C/C++) so that all of this, in all its Byzantine glory, can be squeezed into a single Orbiter compatible Lua script or MFD that works smoothly and efficiently. That process of consolidation is work for another day, and so for the time being, I'll persist with working with the tools at hand - even if the language used to express these tools isn't the same.

As with the Earth-Moon example, we first need to set define an MJD for which we want to calculate the state vectors of a vessel on a prescribed Halo orbit. Rather arbitrarily, let's choose 58849.5000 (1 Jan, 2020). We also note that in the Sun-Mars system, the Sun is the primary and Mars is the secondary gravitating body. With this information, we know modify the Lua script used to interrogate Orbiter's system to read:

Code:
-- set the system MJD
mjd = 58849.500000
oapi.set_simmjd(mjd)

-- set the gravitational constant
G       = 6.67259e-11
primary = oapi.get_objhandle("Sun")
secondary = oapi.get_objhandle("Mars")

-- get the current state vectors and gravitational parameter of the primary gravitating body
q_pri   = oapi.get_globalpos(primary)
p_pri   = oapi.get_globalvel(primary)
m_pri   = oapi.get_mass(primary)
mu_pri  = G * m_pri
	
-- get the current state vectors and gravitational parameter of the secondary gravitating body
q_sec   = oapi.get_globalpos(secondary)
p_sec   = oapi.get_globalvel(secondary)
m_sec   = oapi.get_mass(secondary)
mu_sec  = G * m_sec

-- calculate the orbital eccentricity vector
mu      = mu_pri + mu_sec
X       = vec.sub(q_sec, q_pri)
V       = vec.sub(p_sec, p_pri)
evec    = vec.mul(vec.dotp(V, V) / mu - 1 / vec.length(X), X)
evec    = vec.sub(evec, vec.mul(vec.dotp(V, X) / mu, V))

-- calculate the true anaomaly
nu      = math.acos(vec.dotp(evec, X) / vec.length(evec) / vec.length(X))
if vec.dotp(X, V) < 0 then nu = 2 * math.pi - nu end

-- print out the results
file = io.open("testmb2012.txt", "w")
io.output(file)
file:write(mjd,"\n")
file:write("evec = \{", evec.x, ", ", evec.y, ", ", evec.z, "\}\n")
file:write("ecc  = ", vec.length(evec), "\n")
file:write("nu   = ", nu, "\n")
file:write("q_pri:  ", q_pri.x, "  ", q_pri.z, "  ", q_pri.y, "\n")
file:write("p_pri:  ", p_pri.x, "  ", p_pri.z, "  ", p_pri.y, "\n")
file:write("q_sec:  ", q_sec.x, "  ", q_sec.z, "  ", q_sec.y, "\n")
file:write("p_sec:  ", p_sec.x, "  ", p_sec.z, "  ", p_sec.y, "\n")
file:write("mu_pri: ", mu_pri, "\n")
file:write("mu_sec: ", mu_sec, "\n")
file:close()

And then we run this Lua script in Orbiter and should produce the following output file ('testmb2012.txt'):

Code:
58849.5
evec = {0.085497719952005, -0.0028886029475584, -0.037738624433631}
ecc  = 0.093500844478127
nu   = 4.152440746448
q_pri:  -570133226.34604  1112852003.9499  3460229.4676238
p_pri:  -14.453627375484  -3.4587836356912  0.4001493733008
q_sec:  -197431126884.1  -132172764664.51  2040245433.3189
p_sec:  14476.909918367  -17993.399624731  -732.16524794741
mu_pri: 1.3271243995509e+020
mu_sec: 42828299163780

This file records: the MJD; the eccentricity vector of Mars relative to the Sun; Mars' orbital eccentricity; Mars' true anomaly; the position and velocity vectors of the primary (the Sun) in Orbiter's global coordinate system; the position and velocity vectors of the secondary (Mars) in Orbiter's global coordinate system; and the gravitational parameters for the Sun and Mars. This information is useful in later parts of the calculation.

Next, we open the spreadsheet calculator of A digression on Halo orbits and select the S-Mars L1 tab. Then, set the value of the value of the orbital eccentricity in cell B18 to 0.093500844478127 and set the scale of [MATH]d\alpha_2[/MATH] in cell B19 to a suitable value - 0.00050, say - and rad of the value of the scale parameter [MATH]d\alpha_1[/MATH] from cell B20. This should be 0.00077147.

[OK, so why choose a value of 0.00050 here where in the earlier EM L2 example we chose 0.03? The values [MATH]d\alpha_1[/MATH] and [MATH]d\alpha_2[/MATH] are expressed in terms of a fraction of the distance between the Sun and Mars in this case (whereas in the EM L2 example they were expressed as a fraction of the distance between the Earth and the Moon. Since the Sun-Mars distance averages around 1.5 AU from the Sun or around 225 million km. So, a scale factor of 0.0005 works out around 110,000 km. So, the scale of the Sun-Mars Halo orbit that we are constructing is of the order of twice that - or roughly 220,000 km. So, not small.]

Anyway, we're now ready to run the Python script modified for the Sun-Mars L1 Halo Orbit:

Code:
import openpyxl
import math
import numpy as np

wb = openpyxl.load_workbook('path to file/Sun-Mars L1 - Lissajous - ER3BP.xlsx')

# extract the data from the spreadsheet
wsgeneral = list(wb['general'    ].rows)
wsomega   = list(wb['omega-data' ].rows)
wslambda  = list(wb['lambda-data'].rows)
wsx       = list(wb['x-data'     ].rows)
wsy       = list(wb['y-data'     ].rows)
wsz       = list(wb['z-data'     ].rows)

mu    = wsgeneral[0][1].value
alpha = wsgeneral[1][1].value


def getOmega(e, da1, da2):
    
    omeg = 0
    for idx in range(1, len(wsomega)):
        x     = wsomega[idx]
        omeg += x[4].value * e**x[1].value * da1**x[2].value * da2**x[3].value
    
    return omeg


def getLambda(e, da1, da2):
    
    lamb = 0
    for idx in range(1, len(wslambda)):
        x     = wslambda[idx]
        lamb += x[4].value * e**x[1].value * da1**x[2].value * da2**x[3].value
        
    return lamb


def pos(e, da1, da2, nu, phi1, phi2):
    
    omeg = getOmega (e, da1, da2)
    lamb = getLambda(e, da1, da2)    
   
    theta1 = omeg * nu + phi1
    theta2 = lamb * nu + phi2
    
    eta = 0
    for idx in range(1, len(wsx)):
        x      = wsx[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        eta   += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    kappa = 0
    for idx in range(1, len(wsy)):
        x      = wsy[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        kappa += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    sigma = 0
    for idx in range(1, len(wsz)):
        x      = wsz[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        sigma += x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    deta = 0
    for idx in range(1, len(wsx)):
        x      = wsx[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        deta  += - w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    dkappa = 0
    for idx in range(1, len(wsy)):
        x      = wsy[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        dkappa+= + w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.cos(angle)
        
    dsigma = 0
    for idx in range(1, len(wsz)):
        x      = wsz[idx]
        angle  = nu * x[4].value + theta1 * x[5].value + theta2 * x[6].value
        w      = x[4].value + omeg * x[5].value + lamb * x[6].value
        dsigma+= -w * x[7].value * e**x[1].value * da1**x[2].value * da2**x[3].value * math.sin(angle)
        
    return [[eta + alpha, kappa, sigma], [deta, dkappa, dsigma]]

This is exactly the same as before but we are now access the data in the Sun-Mars L1 file. This file can be downloaded from Sun-Mars General Lissajous Orbit Datasets. To run this code, we need to execute:

Code:
## use the parameter values calculated earlier
e    = 0.093500844478127
nu   = 4.152440746448
da1  = 0.00077147
da2  = 0.00050
phi1 = 0.0
phi2 = math.pi

p = pos(e, da1, da2, nu, phi1, phi2)
print(p)

where we have substituted the values for 'e', 'nu', 'da1' and 'da2' determined earlier. This code snippet should produce:

Code:
[[0.995079955178846, -0.0017203374881442496, 0.0007504952355208614], [-0.0014543153527665142, 0.002918884110869531, 0.0013498758108830297]]

which is a list of two 3-vectors, the first is the position vector of the spacecraft on the prescribed Halo orbit for the given MJD expressed in the rotating-pulsating coordinates of the ET3BP theory; and the second vector is the same for the velocity vector of the vessel. To use these results, we need to convert from the coordinate system of the ER3BP to Orbiter's reference frame. We do this by running the following Python code:

Code:
# input the epoch MJD
mjd = 58849.500000

# input the state vectors of the primary gravitating body - right-handed coordinate system
Qpri  = np.array([-570133226.34604,  1112852003.9499,  3460229.4676238])
Ppri  = np.array([-14.453627375484,  -3.4587836356912,  0.4001493733008])

# input the state vectors of the secondary gravitating body - right-handed coordinate system
Qsec  = np.array([-197431126884.1,  -132172764664.51,  2040245433.3189])
Psec  = np.array([14476.909918367,  -17993.399624731, -732.16524794741])

# set up the gravitational constants
mupri = 1.3271243995509e+020
musec = 42828299163780.0
mu    = mupri + musec

def getECIcoords(a, e, nu, pos, mu):
    
    [[eta, kappa, sigma], [deta, dkappa, dsigma]] = pos
    cs = math.cos(nu)
    sn = math.sin(nu)

    f  = a * (1 - e*e) / (1 + e * cs)
    x  = f * (eta * cs - kappa * sn)
    y  = f * (eta * sn + kappa * cs)
    z  = f * sigma
    
    g  = math.sqrt(mu / a / (1-e*e))
    vx = g * (-sn*eta   - (e + cs)*kappa + (1 + e*cs)*(cs*deta   - sn*dkappa))
    vy = g * (-sn*kappa + (e + cs)*eta   + (1 + e*cs)*(cs*dkappa + sn*deta  ))
    vz = g * ( e*sn*sigma + (1 + e*cs)*dsigma)
    
    return [[x, y, z], [vx, vy, vz]]


# clculate the centre of mass and the centre of velocity
com = (mupri * Qpri + musec * Qsec) / mu
cov = (mupri * Ppri + musec * Psec) / mu

r   = Qsec - Qpri
v   = Psec - Ppri
R   = np.linalg.norm(r)

# calculate a, e, nu
evec= r * np.dot(v,v)/mu - v * np.dot(r,v)/mu - r/R
e   = np.linalg.norm(evec)
a   = mu / (2*mu / R - np.dot(v,v))
nu  = math.acos(np.dot(evec,r)/e/R)
if np.dot(r,v)<0:
    nu = 2 * math.pi - nu
    
xhat= evec / e
zhat= np.cross(r, v)
zhat= zhat / np.linalg.norm(zhat)
yhat= np.cross(zhat,xhat)

# now calculate the ECI state vectors relative to the secondary - right-handed coordinate system
[[Qx, Qy, Qz], [Px, Py, Pz]] = getECIcoords(a, e, nu, p, mu)
Q = Qx * xhat + Qy * yhat + Qz * zhat + com - Qsec
P = Px * xhat + Py * yhat + Pz * zhat + cov - Psec

# print the MJD
print(mjd)

# print the position vector relative to the primary - (Orbiter-freindly left-handed coordinate system, ecliptic, cartesian)
print([Q[0], Q[2], Q[1]])

# print the velocity vector relative to the primary - (Orbiter-friendly left-handed coordinate system, ecliptic, cartesian)
print([P[0], P[2], P[1]])

Here we have set up the MJD, the state vectors of the Sun and Mars as per the earlier Lua output, and the gravitational parameters of the Sun and Mars again as per the earlier Lua output. For convenience, we've also set the state vectors of the vessel relative to Mars rather than the Sun

On running, this produces:

Code:
58849.5
[743773315.5467529, 181031161.34134054, 990567937.2523193]
[-36.1315008765705, 31.09794058234411, 26.062745990879193]

The first line gives the MJD; the second is the position vector (metres) of the vessel relative to Mars; and the third is the velocity vector (m/s) relative to Mars. These results can now be substituted directly into Orbiter. And we can translate these results in a simple scenario file for use with Orbiter:

Code:
BEGIN_ENVIRONMENT
  System Sol
  Date MJD 58849.500000
END_ENVIRONMENT

BEGIN_FOCUS
  Ship GL-NT
END_FOCUS

BEGIN_SHIPS
GL-NT:DeltaGlider
  STATUS Orbiting Mars
  RPOS 743773315.5467529 181031161.34134054 990567937.2523193
  RVEL -36.1315008765705 31.09794058234411 26.062745990879193
  RCSMODE 2
  AFCMODE 7
  PRPLEVEL 0:1.000000 1:1.000000
  NAVFREQ 0 0 0 0
  XPDR 0
  HOVERHOLD 0 1 0.0000e+000 0.0000e+000
  AAP 0:0 0:0 0:0
END
END_SHIPS
 
Top