Introduction

The R programming language can produce Web-based Space Mission Visualizations (WSMV). With the rgl package, a space mission designer can create interactive animated 3D scenes and export JavaScript and WebGL. With the R Markdown language, a WSMV developer can design documents and with the RPubs service provided by the makers of RStudio, the WSMV can be published online.

Scope

The Physics department in College of Saint Benedict+Saint John’s University explains the geometry for 2D and 3D Orbits and provides example Mathematica code. Essentially, those web-pages describe a procedure of specifying a 2D ellipse and then applying 3D rotations to the ellipse using the Keplarian parameters. To gain some insight to the 2D geometry, visit http://www.physics.csbsju.edu/orbit/orbit.2d.html.then click on the link to the 3D orbit web-page.http://www.physics.csbsju.edu/orbit/orbit.3d.html This tutorial applies the methodology of specifying a 2D ellipse and then applying pitch, yaw, and roll rotation matrices to orient the ellipse in three dimensions. The following image depicts these transformations.

Generate a flat surface ellipse

Generating an ellipse involves setting a value for a semi-major axis, a, and eccentricity, e. With an e = 0, the ellipse would be a circle; with an e = 1, the ellipses becomes a line. A value between zero and one gives an ellipse a shape somewhere between a circle and a line. The following chunk of R code generates an ellipse.

#Initialize
a <- 1                  # semi-major axis
e <- 1/sqrt(2)          # eccentricity
b <- a * sqrt(1 - e^2)  # semi-minor axis
c <- e * a              # distance from the center to a focus

# Distance from a focus to the center.
x.c <- -c
y.c <- 0
z.c <- 0

# Generate a numerical sequence
u <- seq(-pi,pi, length.out=80)

# Generate x and y values.
x <- a * cos(u) - e
y <- a * sqrt(1-e^2) * sin(u)
z <- rep(0, times=80)

c.x <- a * cos(u) - c
c.y <- a * sin(u)         # Generate points for an outer circle.
c.z <- rep(0, times=80)

# Plot x & y with a line.
plot(x,y,type="l", ylim=c(-a*sin(pi/2), a*sin(pi/2)),asp=1)
points(c.x,c.y,type="l", col="green")

Transform the flat ellipse into a 3D orbit

Install the R package named rgl. Documentation for rgl is available at https://cran.r-project.org/web/packages/rgl/vignettes/rgl.html. The following R code performs three rotations. Think of an orbit as an airplane. A plane can pitch up or down, it can yaw, i.e., rotate around its vertical axis, and it can roll. In addition to the semi-major axis and eccentricity, other Keplarian parameters include inclination, longitude of the ascending node, ??, and the Right Ascension of the Ascending Node (RAAN). Using inclination for pitch, ?? for yaw, and RAAN for roll, we can rotate the flat ellipse around the X, Y, and Z axes. The following R code initializes three arrays, x.inc, y.inc, z.inc, to hold the values for the inclined ellipse. A for loop steps through each of the values stored in the x, y, z values and applies a 3D rotation matrix. Parameters of the rotate3d function include the object, the angle, and the axis. In this case, the object is one of the Cartesian coordinates, an angle of pi/5, and the 0, 1, 0 represents the Y axis. After calculating the new point, the x, y, and z values are stored in the x.inc, y.inc, and z.inc arrays. The rest of the code on this page conducts similar transformations for yaw, i.e. ??, and roll, i.e. RAAN.

# Pitch ~ inclination, tilt about the y axis.
x.inc <- rep(0, times=80)
y.inc <- rep(0, times=80)
z.inc <- rep(0, times=80)

for (i in 1:80) {
 point <- rotate3d(c(x[i],y[i],z[i]),pi/5,0,1,0)
 x.inc[i] <- point[1]
 y.inc[i] <- point[2]
 z.inc[i] <- point[3]
}
# Transform the focus for the central body.
center <- rotate3d(c(x.c,y.c,z.c),pi/5,0,1,0)

# Yaw ~ longitude of the ascending node (omega),
# so rotate around the z axis.
x.om <- rep(0, times=80)
y.om <- rep(0, times=80)
z.om <- rep(0, times=80)

for (i in 1:80) {
 point <- rotate3d(c(x.inc[i],y.inc[i],z.inc[i]),pi/4,0,0,1)
 x.om[i] <- point[1]
 y.om[i] <- point[2]
 z.om[i] <- point[3]
}
# Transform the focus for the central body.
center <- rotate3d(c(center[1],center[2],center[3]),pi/4,0,0,1)

# Roll ~ Right ascension of the ascending node (RAAN),
# so rotate around the x axis.
x.raan <- rep(0, times=80)
y.raan <- rep(0, times=80)
z.raan <- rep(0, times=80)

for (i in 1:80) {
 point <- rotate3d(c(x.om[i],y.om[i],z.om[i]),pi/4,1,0,0)
 x.raan [i] <- point[1]
 y.raan [i] <- point[2]
 z.raan [i] <- point[3]
}
# Transform the focus for the central body.
center <- rotate3d(c(center[1],center[2],center[3]),pi/4,1,0,0)

Plot the transformed ellipses

Each of transformations stored the values in new variables. The following code sets parameters for a 3D device, establishes layout of one row with three columns, and plots interactive image of the transformed ellipses. The rglwidget function embeds the interactive plot within a web-page.

You can click-drag to rotate and right-click-drag to zoom in and out the following 3D plots.

# Plot the transformed elliptical orbit.
par3d(zoom=0.8, windowRect = c(200, 5, 800, 400))
mfrow3d(nr = 1, nc = 3, sharedMouse = TRUE) 

plot3d(x=x.inc,y=y.inc,z=z.inc, type="l", col="red")
next3d()
plot3d(x=x.om,y=y.om,z=z.om, type="l", col="blue")
next3d()
plot3d(x=x.raan ,y=y.raan ,z=z.raan, type="l", col="green")
rglwidget()

Pitch ~ Inclination,Yaw ~ omega,and Roll ~ RAAN

Now, let’s place a blue sphere at the origin and a little red sphere at the center of the ellipse. You can click-drag to rotate and right-click-drag to zoom in and out the following 3D plots.

open3d(useNULL = rgl.useNULL())          # open another 3d device
## wgl 
##   3
par3d(zoom=0.8)     # set a parameter
plot3d(x.raan,y.raan,z.raan,type="l", col="red")
spheres3d(center[1],center[2],center[3],radius=0.01, col="red")
spheres3d(0,0,0,radius=0.1, col="blue")
rglwidget()

Orbit around a sphere

Propagate a Keplerian orbit

One of Kepler’s laws states that an orbiting object sweeps out equal areas in equal amounts of time. What this law means is that an orbiting object accelerates as it approaches the periapsis and decelerates as it approaches apoapsis. This changinge rate is reflected in Kepler’s equation, E(t) - esin(E) = M(t), which states that the Mean anomaly at time is equal to the Eccentric anomaly at time t minus the eccentricity multiplied by the sine of the Eccentric anomaly. This equation is transcendental because it transcends the simple equations; it must be solved iteratively. Marc Murison developed a practical method for solving Kepler’s equation using a Taylor series expanded to the third order. The following chunk of R code implements Murison’s algorithm.

After the implementation of the algorithm, a function named propagate calls the Keplerian equation solver to generate the True anomaly, which is the current position of an object in an orbit. Code at the end of the chunk uses the propagate function to populate an array of vertices and then implements a playwidget for an interactive 3D model of an orbiting object.

KeplerStart3 <- function(e,M) {
t34 <- e^2
t35 <- e * t34
t33 <- cos(M)

result <- M+(-1/2*t35+e+(t34+3/2 * t33*t35)*t33)*sin(M)
return(result)
}

eps1 <- function(e,M,x) {
 result <- (x - e * sin(x) - M)/(1 - e * cos(x))
 return(result)
}

eps2 <- function(e,M,x) {
  t1 <- -1 + e * cos(x)
  t2 <- e * sin(x)
  t3 <- -x + t2 + M
  result <- t3/(1/2 *t3 * t2/t1 + t1)
  return(result)
}

eps3 <- function(e,M,x) {
   t1 <- cos(x)
   t2 <- -1 + e * t1
   t3 <- sin(x)
   t4 <- e*t3
   t5 <- -x + t4 + M
   t6 <- t5/(1/2 * t5 * t4/t2 + t2)
   result <- t5/((1/2 * t3 - 1/6 * t1 * t6) * e * t6 + t2)
   return(result)
}

KeplerSolve <- function(e, M) {
   tol <- 1.0e-14
   Mnorm <- M %% 2 * pi
   E0 <- KeplerStart3(e, Mnorm)
   dE <- tol + 1
   count <- 0

   while (dE > tol) {
      E <- E0 - eps3(e, Mnorm, E0)
      dE <- abs(E - E0)
      E0 <- E
      count <- count + 1
      
      if (count == 100) {
         msg <- "Astounding! KeplerSolve failed to converge!"
         print(msg)
         break
      } #endif
   } # end while

   return(E)
}

 propagate <- function(clock) {
  T <- 120      # seconds
  n <- 2*pi/T
  tau <- 0      # time of pericenter passage

  M <- n * (clock - tau)
  E <- KeplerSolve(e,M)
  cose <- cos(E) # shows up multiple times so calculate once.

  r = a*(1 - e * cose)
  s.x <- r * ((cose - e)/(1 - e * cose))
  s.y <- r * ((sqrt(1 - e^2)*sin(E))/(1 - e * cose))
  s.z <- 0
  
  point <- rotate3d(c(s.x,s.y,s.z),pi/5,0,1,0)
  point <- rotate3d(c(point[1],point[2],point[3]),pi/4,0,0,1)
  point <- rotate3d(c(point[1],point[2],point[3]),pi/4,1,0,0)

  return(point)
 }
 
par3d(zoom=0.8)     # set a parameter
orbit <- plot3d(x.raan,y.raan,z.raan,type="l", col="black")
spheres3d(0,0,0,radius=0.2, col="blue")

ts <- 120  # Number of time slices.
orbcoords <- cbind(1:ts,0,0) # Allocate a three column array.
for (clock in 1:ts) {
 loc <- propagate(clock)
 orbcoords[clock,1] <- loc[1]
 orbcoords[clock,2] <- loc[2]  # Populate the array of vertices.
 orbcoords[clock,3] <- loc[3]
}
# Put a small sphere in the orbit.
 sphereid <- spheres3d(xyz.coords(x=loc[1],y=loc[2],z=loc[3]), 
             radius = 0.04, col = "green")
 
clock <- 1:ts
rglwidget() %>%
  playwidget(list(
  ageControl(births = 1, ages = clock, vertices = orbcoords, objids = sphereid)),
  start = 1, stop = max(clock), rate = 5,
  components = c("Reverse", "Play", "Slower", "Faster",
                 "Reset", "Slider", "Label"),
  loop = TRUE)

References

  1. Calvert, J.B., Definition of the Ellipse, 2002. http://mysite.du.edu/~jcalvert/math/ellipse.htm
  2. rgl package https://cran.r-project.org/web/packages/rgl/rgl.pdf
  3. Mfrow3d https://www.rdocumentation.org/packages/rgl/versions/0.97.0/topics/mfrow3d
  4. R Markdown Reference http://www.rstudio.com/wp-content/uploads/2015/03/rmarkdown-reference.pdf
  5. knitr chunk options https://yihui.name/knitr/options/
  6. Murison, Marc A., “A Practical Method for Solving the Kepler Equation”, U.S. Naval Observatory, Washington, DC, November 6, 2006.http://alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf