Skip to content

Instantly share code, notes, and snippets.

@atarp

atarp/kepler22.R

Last active May 17, 2022 21:21
Show Gist options
  • Save atarp/5504e5cbfe82b0602cb5b140a603dd28 to your computer and use it in GitHub Desktop.
Save atarp/5504e5cbfe82b0602cb5b140a603dd28 to your computer and use it in GitHub Desktop.
Download ZIP
New Version of Kepler
Raw
kepler22.R
# years ago I found the original version on R-Bloggers:
# https://www.r-bloggers.com/2014/04/hazardous-and-benign-space-objects-getting-the-data/
# Now the number of NEA-Objekts in the MPC database has doubled an they provide the date
# in different formats
library(tidyverse)
neafile<-"nea.txt"
classes<-c(
"Atira",
"Aten",
"Apollo",
"Amor",
"ObjectQ", # with q < 1.665 AU
"Hungaria",
"Unused", # or internal MPC use only
"Hilda",
"Jupiter Trojan",
"Distant object"
)
classes_abrv<-c(
"ATI",
"ATE",
"APO",
"AMO",
"OBJQ", # with q < 1.665 AU
"HUN",
"UNU", # or internal MPC use only
"HIL",
"TRO",
"DIST"
)
## description: https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html => fortran format string
## 1:7, 9:13, 15:19, 21:25, 27:35, 38:46, 49:57, 60:68, 71:79 81:91 93:103, 106:160, 162:165, 167:194, 195:202
## ObjNum| H | G | Epoch | M || w || Node || i || e | P | a || unused | class |ObjName last
formatF90<-c("A7","X","F5.2","X","F5.2","X","A5","X","F9.5","2X","F9.5","2X","F9.5","2X","F9.5","2X","F9.7","X","F11.8","X","F11.7","2X","A55","X","A4","X","A28","I8")
## wget https://www.minorplanetcenter.net/iau/MPCORB/NEA.txt
delta_t<-as.numeric(strftime(Sys.time(),format = "%d"))-as.numeric(strftime(file.info(neafile)$mtime,format = "%d"))
if(delta_t>2){curl::curl_download(url="https://www.minorplanetcenter.net/iau/MPCORB/NEA.txt",destfile = neafile,quiet = TRUE)}
df<-tibble(read.fortran(neafile,format = formatF90))%>%select(2:11,13:14)
names(df)<-c("H","G","Epoch","M","w","Node","i","e","P","a","class" , "Object")
# there are some caveats with the fortran format e.g. F5.2 reads the values / 100, F5.6 reads values*1e-6 ...
rad<-pi/180.0*1e5
degy_drev<-360.0*1e-8/365.2425 #degrees per day * years per revolution
df<-df%>%
mutate(H=H*1e2, G=G*1e2, i=i*rad, w=w*rad, Node=Node*rad, M=M*rad,e=e*1e7, P=degy_drev/P, a=a*1e7,
hazard=factor(as.character(as.integer(paste("0x",class,sep=''))>8*2^12),levels=c(TRUE,FALSE)),
class=classes_abrv[as.integer(paste("0x",class,sep=''))%%256],Object=str_trim(gsub("\\(|\\)", "", Object)))%>%
select(H:Object,hazard)
df<-df[,c(12,3,10,8,7,5,6,4,9,1,2,11,13)]
# if we need the epoch, we first have to unpack it
#for (c in str_split(df$epoch,"")[[1]]){ if(length(which(LETTERS==c))>0) print(which(LETTERS==c)+9) else print(as.numeric(c))}
head(df)
#
# more or less the same code from here
#
#Statistics
table(df$class)
with(df, table(class, hazard))
#Histogram of Absolute V-Magnitude ( H<30.0 chosen )
ggplot(df[df$H<30.0, ], aes(x = H)) +
geom_histogram(binwidth = 1.0, fill = "lightgray", colour = "black") +
xlab("Absolute V-Magnitude") + ylab("") +
theme_classic()
#ggsave("/artifacts/MagVHist.png")
#Kepler's 3rd law
#show the underlying power law in a loglog plot
ggplot(df, aes(x = a, y = P)) +
geom_point() +
scale_x_log10() +
scale_y_log10() +
xlab("a [AU]") +
ylab("Period [year]") + theme_classic()
#ggsave("/artifacts/powerlaw.png")
#Levenberg-Marquardt verifies a slope of 1.5
stats::lm(log10(P) ~ log10(a), data = df)
#Histogram of Eccentricities
ggplot(df, aes(x = e)) +
geom_histogram(binwidth = 0.05, fill = "lightblue", colour = "black") +
xlab("Eccentricity") + ylab("") +
theme_classic()
#ggsave("/artifacts/eccentr.png")
#Semi Major Axis vs. Eccentricity
#subdivided into the four classes and weather they are
#potentially harzardous or not
ggplot(df[order(df$hazard),], aes(x = e, y = a)) +
geom_point(aes(colour = class), alpha = 0.5, size = 2) +
scale_y_log10() +
scale_colour_discrete("Class", h = c(240, 0), c = 150) +
facet_wrap(~ hazard, nrow = 1) +
xlab("Eccentricity") + ylab("Semi-Major Axis [AU]") +
theme_classic()
#ggsave("/artifacts/AvsEccClasses.png")
library(nleqslv)
library(plyr,warn.conflicts = F)
library(scales,warn.conflicts = F)
# M = E - e * sin E.
# Solve Kepler's Equation to find the eccentric anomaly, E.
# This is a transcendental equation so we need to solve it
# numerically using the Newton-Raphson method.
# Append another column "E"
df = ddply(df, .(Object), mutate,
E = {
kepler <- function(E) {
M - E + e * sin(E)
}
nleqslv(0, kepler, method = "Newton")$x
})
head(df)[, c("Object", "E")]
degrees <- function(x, ...) {parse(text = paste0(x, "*degree"))}
# plot Eccentric Anomaly (E) vs. Mean Anomaly (M)
ggplot(df, aes(x = M / pi * 180, y = E / pi * 180, colour = e)) +
geom_point() +
scale_colour_gradient("Eccentricity (e)", limits=c(0, 1), high = "red", low = "blue") +
xlab("Mean Anomaly (M)") + ylab("Eccentric Anomaly (E)") +
scale_x_continuous(labels = degrees, breaks = seq(0, 360, 90)) +
scale_y_continuous(labels = degrees, breaks = seq(0, 360, 90)) +
annotate("text",x=90,y=270,label="M == E - e%.%sin(E)", parse = TRUE) +
theme_classic()
#ggsave("/artifacts/MeanAnovsEccAno.png")
# compute and fill in the True Anomaly theta
df = transform(df, theta = 2 * atan2(sqrt(1 + e) * sin (E/2), sqrt(1 - e) * cos (E/2)))
#head(orbits)[, c("Object", "theta")]
# plot True Anomaly (E) vs. Mean Anomaly (M)
ggplot(df, aes(x = M / pi * 180, y = theta / pi * 180, colour = e)) +
geom_point() +
scale_colour_gradient("Eccentricity (e)", limits=c(0, 1), high = "red", low = "blue") +
xlab("Mean Anomaly (M)") + ylab(expression(paste("True Anomaly ( ",theta,")"))) +
scale_x_continuous(labels = degrees, breaks = seq(0, 360, 90)) +
scale_y_continuous(labels = degrees, breaks = seq(0, 360, 90)) +
annotate("text",x=90,y=270,label="tan(theta/2) == sqrt(frac(1+e,1-e))%.%tan(E/2)" , parse = TRUE) +
theme_classic()
#ggsave("/artifacts/MeanAnovsTrueAno.png")
#####################################################################
#Part III
#####################################################################
#Transform between heliocentric ecliptic and orbital plane
df = within(df, {
Z = 0
Y = a * sqrt(1 - e**2) * sin(E)
X = a * cos(E) - a * e
})
head(df)
transformation.matrix <- function(w, Node, i) {
a11 = cos(w) * cos(Node) - sin(w) * cos(i) * sin(Node)
a12 = cos(w) * sin(Node) + sin(w) * cos(i) * cos(Node)
a13 = sin(w) * sin(i)
a21 = - sin(w) * cos(Node) - cos(w) * cos(i) * sin(Node)
a22 = - sin(w) * sin(Node) + cos(w) * cos(i) * cos(Node)
a23 = cos(w) * sin(i)
a31 = sin(i) * sin(Node)
a32 = - sin(i) * cos(Node)
a33 = cos(i)
matrix(c(a11, a12, a13, a21, a22, a23, a31, a32, a33), nrow = 3, byrow = TRUE)
}
(A = do.call(transformation.matrix, df[1, c("w", "Node", "i")]))
# We want to invert this transformation, so we need to invert the matrix. However, since this is a
# rotation matrix, the inverse and the transpose are the same... (we will use this fact later!)
A= solve(A)
A %*% matrix(unlist(df[1, c("X", "Y", "Z")]), ncol = 1)
require(plyr)
df = ddply(df, .(Object), function(df2) {
A = transformation.matrix(df2$w, df2$Node, df2$i)
#
# Invert transformation matrix
#
A = t(A)
#
# Transform to reference frame
#
r = A %*% matrix(unlist(df2[, c("X", "Y", "Z")]), ncol = 1)
#
r = matrix(r, nrow = 1)
colnames(r) = c("x", "y", "z")
#
cbind(df2, r)
})
ggplot(df, aes(x = x, y = y)) +
geom_point(alpha = 0.25) +
xlab("x [AU]") + ylab("y [AU]") +
scale_x_continuous(limits = c(-35, 35)) + scale_y_continuous(limits = c(-55, 15)) +
theme_classic() +
theme(aspect.ratio = 1)
#ggsave("/artifacts/distri1.png")
ggplot(df, aes(x = x, y = y)) +
geom_point(aes(colour = hazard), alpha = 0.25) +
scale_colour_discrete("Hazard", h = c(180, 0), c = 150, l = 40) +
xlab("x [AU]") + ylab("y [AU]") +
scale_x_continuous(limits = c(-5, 5)) + scale_y_continuous(limits = c(-5, 5)) +
annotation_custom(grob = grid::circleGrob(r = unit(0.5, "npc"),
gp = grid::gpar(fill = "transparent", col = "black", lwd = 2, lty = "dashed")),
xmin=-1, xmax=1, ymin=-1, ymax=1) +
theme_classic() +
theme(aspect.ratio = 1)
#ggsave("/artifacts/distriTopV.png")
ggplot(df, aes(x = x, y = z)) +
geom_point(aes(colour = hazard), alpha = 0.25) +
scale_colour_discrete("Hazard", h = c(180, 0), c = 150, l = 40) +
xlab("x [AU]") + ylab("z [AU]") +
scale_x_continuous(limits = c(-5, 5)) + scale_y_continuous(limits = c(-5, 5)) +
theme_classic() +
theme(aspect.ratio = 1)
#ggsave("/artifacts/distriPV.png")
library(scatterplot3d)
#graphics.off()
#png(filename = "/artifacts/distri3d.png", width = 8, height = 6, units = "in", res = 300)
#par(mai = c(0.5, 0.5, 0.5, 0.5))
df.3d <- with(df, scatterplot3d(x, y, z, xlim = c(-5, 5), ylim = c(-5, 5), highlight.3d = TRUE))
df.3d$plane3d(c(0, 0, 0))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment