3d Plot of earth, moon, and sun over 1 year

Schrick-Noah_Homework-4.R

@@ -96,7 +96,6 @@ dho.f <- function(t, y, kpass){
   ydot <- y[2]
   #         c         k         m
   ypp <- (-kpass[2]*yn - kpass[3]*ydot)/kpass[1]
-  #ypp <- -kpass[2]*y - kpass[3]*ydot
   as.matrix(c(ydot,ypp))
 }
 
@@ -149,18 +148,11 @@ legend("topright", # coordinates, "topleft" etc
 )
 
 ## 3. Use ode45/rk4sys to solve for the traj of a projectile thrown vertically 
-
 # a. IVP
 tmin <- 0
 tmax <- 2.04
 xo <- 0
 vo <- 10
-
-# b. BVP
-yo <- 0
-ymax <- 0
-
-library(pracma)
 projectile.f <- function(t,y){
   g <- 9.81
   v <- y[2]  # y1dot
@@ -168,6 +160,29 @@ projectile.f <- function(t,y){
   matrix(c(v,a))
 }
 
+y.init <- c(xo, vo)
+ivp.sol <- ode45(f=function(t,y){projectile.f(t,y)},
+                 y=y.init, t0=tmin, tfinal=tmax)
+ivp.sol.rk4 <- rk4sys(projectile.f, 0, 2.04, c(0, 10), 100)
+
+plot(ivp.sol$t, ivp.sol$y[,1] ,type="o",col="blue", 
+     ylim = c(-0.5, 5.5), xlab="time",ylab="height")
+par(new=T)
+lines(ivp.sol.rk4$x, ivp.sol.rk4$y[,1] ,type="o",col="red", 
+     xlab="time",ylab="height")
+legend("topright", # coordinates, "topleft" etc
+       c("ODE45, default","RK4Sys, 100 steps"), # label
+       lty=c(1,1), # line 
+       lwd=c(2.5,2.5), # weight
+       #cex=.8, 
+       bty="n", # no box
+       col=c("blue","red", "green") # color
+)
+
+# b. BVP
+yo <- 0
+ymax <- 0
+
 proj.obj <- function(v0, y0=0, tfinal){
   # minimize w.r.t. v0
   proj.sol <- ode45(projectile.f, 
@@ -186,18 +201,53 @@ v_best$minimum  # best v0
 
 best.sol <- rk4sys(projectile.f, a=0, b=10, y0=c(0, v_best$minimum),
                    n=20)  # 20 integration stepstmax 
+plot(best.sol$x, best.sol$y[,1], type="l")
+
+computeHeights <- function(xo, vo, tmin){
+  height <- xo + vo*tmin -0.5*9.81*tmin^2
+}
+
+height <- computeHeights(best.sol$x, best.sol$y[,1], seq(0,10,len=21))
 
 # c. Shooting method for damped oscillator with perturbation parameter
 yo <- 0
 y1 <- 1
-ymin <- 0
+tmin <- 0
-ymax <- 2
+tfinal <- 2
 
 pfirst <- 0.5
 psec <- 0.05
 
+dho.pert.f <- function(t, y, kpass){
+  yn <- y[1]
+  ydot <- y[2]
+  ep <- kpass[1]
+  #         c                           m
+  ypp <- (-(1+ep)*yn - ydot)/ep
+  as.matrix(c(ydot,ypp))
+}
+
+proj.obj <- function(v0, y0=yo, tfinal){
+  # minimize w.r.t. v0
+  proj.sol <- ode45(dho.pert.f, 
+                  y=c(yo, y1), kpass=pfirst, t0=tmin, tfinal=tfinal)
+  final_index <- length(proj.sol$t)
+  yf <- proj.sol$y[final_index,1]  # want equal to right boundary 
+  log(abs(yf)) # minimize this
+}
+
+# user specifies tfinal and yfinal for BVP
+ep_best <- optimize(proj.obj, 
+                   interval=c(1,100), #bisect-esque interval
+                   tol=1e-10, 
+                   y0=0, tfinal=2) # un-optimized obj params
+ep_best$minimum  # best v0
+
+best.sol <- rk4sys(dho.pert.f, a=0, b=2, y0=c(0, ep_best$minimum),
+                   n=20)  # 20 integration stepstmax 
+
 # analytical sol
-pthird <- 0
+# yprime = -y
 
 ## 4. Position of the earth and moon
 # a. Plotly
@@ -216,6 +266,56 @@ x0_moon <- c(-27083318944, 133232649728, 57770257344)/1e3 # km
 v0_moon <- c(-30864.2207031, -4835.03349304, -2042.89546204)*(86400)/1e3 
 # m/s -> km/day
 
-# b. Find eclipses
+sunearth.f <- function(t,y){
+  G <- 6.673e-11 # m^3 kg^-1 s^-2
+  M_S <- 1.9891e30
+  M_E <- 5.98e24
+  M_m <- 7.32e22   # not using right now
+  mu_sun <- G*M_S*(86400^2)/1e9  # km^3/days^2
+  
+  x_earth <- y[1:3]  # xyz position of earth wrt sun
+  v_earth <- y[4:6]  # xyz velocity of earth, also part of ydot
+  r_earth <- sqrt(sum(x_earth^2))  # distance earth to sun
+  acc_earth <- -mu_sun*x_earth/r_earth^3  # part of ydot
+  
+  ydot <- c(v_earth, acc_earth)
+  return(matrix(ydot))
+}
 
+y.init.earth <- c(x0_earth, v0_earth)
+sunearth.sol <- rk4sys(f=sunearth.f, a=0, b=365.25, y0=y.init.earth, n=365)
+sunearth.df <- data.frame(x=sunearth.sol$y[,1],
+                          y=sunearth.sol$y[,2],
+                          z=sunearth.sol$y[,3])
+
+y.init.moon <- c(x0_moon, v0_moon)
+sunmoon.sol <- rk4sys(f=sunearth.f, a=0,b=365.25, y0=y.init.moon, n=365)
+sunmoon.df <- data.frame(x=sunmoon.sol$y[,1],
+                          y=sunmoon.sol$y[,2],
+                          z=sunmoon.sol$y[,3])
+
+if (!require("plotly")) install.packages("plotly")
+library(plotly)
+fig <- plot_ly(sunmoon.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("red"), name="sun")
+fig <- add_trace(fig, x = sunearth.df[1]$x, y = sunearth.df[2]$y, z = sunearth.df[3]$z, mode="markers", color = I("green"), name="earth")
+fig
+
+mean(sqrt(rowSums(cbind(sunearth.sol$y[,1]^2, 
+                        sunearth.sol$y[,2]^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))))
+# Not correct: giving 60m km - should be ~380k
+
+# b. Find eclipses
+norm_vec <- function(x) sqrt(rowSums(cbind((x^2))))
+eclipses <- norm_vec(sunearth.sol$y) %*% norm_vec(sunmoon.sol$y)
+                            
 # c. keep orbiter at L2 Lagrange point for a year, ignoring the moon's effect