noah/CS-7863-Sci-Stat-Proj-4
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Finalizing orbital mechanics

Browse Source
This commit is contained in:
Noah L. Schrick 2023-03-01 20:52:04 -06:00
parent 1dd9439d66
commit a503e9fffe
1 changed files with 32 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
Schrick-Noah_Homework-4.R
View File

@ -259,11 +259,6 @@ plot(best.sol$x, best.sol$y[,1], type="l")
## 4. Position of the earth and moon
# a. Plotly
G <- 6.673e-11 # m^3 kg-1 s^-2
M_S <- 1.9891e30 # sun kg
M_E <- 5.98e24 # earth kg
M_m <- 7.32e22 # moon kg
mu_sun <- G*M_S*(86400^2)/1e9 # km^3/days^2
# t0: jan 1, 1999 00:00:00am
x0_earth <- c(-27115219762.4, 132888652547.0, 57651255508.0)/1e3 # km
@ -317,15 +312,39 @@ mean(sqrt(rowSums(cbind(sunearth.sol$y[,1]^2,
sunearth.sol$y[,3]^2))))
# 150 million km
mean(sqrt(rowSums(cbind((sunearth.sol$y[,1] - sunmoon.sol$y[,1])^2,
(sunearth.sol$y[,2] - sunmoon.sol$y[,2])^2,
(sunearth.sol$y[,3] - sunmoon.sol$y[,3])^2))))
sunmoon.f <- function(t,y){
G <- 6.673e-11 # m^3 kg^-1 s^-2
M_E <- 5.98e24
M_m <- 7.32e22
mu_sun <- G*M_E*(86400^2)/1e9 # km^3/days^2
mean(sqrt(rowSums(cbind(sunmoon.sol$y[,1]^2,
sunmoon.sol$y[,2]^2,
sunmoon.sol$y[,3]^2))))
# Not correct: giving 60m km - should be ~380k
# Why: Don't think M_m is used in function
x_moon <- y[1:3] # xyz position of moon wrt earth
v_moon <- y[4:6] # xyz velocity of moon, also part of ydot
r_moon <- sqrt(sum(x_moon^2)) # distance moon to earth
acc_moon <- -mu_sun*x_moon/r_moon^3 # part of ydot
ydot <- c(v_moon, acc_moon)
return(matrix(ydot))
}
x0_moonearth <- x0_moon - x0_earth
v0_moonearth <- v0_moon - v0_earth
y.init.moonearth <- c(x0_moonearth, v0_moonearth)
earthmoon.sol <- rk4sys(f=sunmoon.f, a=0,b=365.25, y0=y.init.moonearth, n=365)
earthmoon.df <- data.frame(x=earthmoon.sol$y[,1],
y=earthmoon.sol$y[,2],
z=earthmoon.sol$y[,3])
mean(sqrt(rowSums(cbind(earthmoon.sol$y[,1]^2,
earthmoon.sol$y[,2]^2,
earthmoon.sol$y[,3]^2))))
fig <- plot_ly(earthmoon.df, x = ~x, y = ~y, z = ~z, name="moon", type = "scatter3d", mode="markers")
fig <- fig %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
fig <- add_trace(fig, x = 0, y = 0, z = 0, mode="markers", color = I("green"), name="earth")
fig
# b. Find eclipses
norm_vec <- function(x) sqrt(rowSums(cbind((x^2))))