L2 Lagrange orbiter

Schrick-Noah_Homework-4.R

@@ -207,7 +207,8 @@ 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))
+height <- computeHeights(yo, v_best$minimum, seq(0,10,len=21))
+plot(seq(0,10,len=21), height, xlab="time", ylab="height")
 
 # c. Shooting method for damped oscillator with perturbation parameter
 yo <- 0
@@ -215,13 +216,11 @@ y1 <- 1
 tmin <- 0
 tfinal <- 2
 
-pfirst <- 0.5
-psec <- 0.05
-
-dho.pert.f <- function(t, y, kpass){
+dho.pert.encap.f <- function(ep){
+  dho.pert.f <- function(t, y){
     yn <- y[1]
-  ydot <- y[2]
-  ep <- kpass[1]
+    ydot <- y[2] - yn
+    #ep <- 0.5
     ypp <- (-(1+ep)*yn - ydot)/ep
     as.matrix(c(ydot,ypp))
   }
@@ -229,24 +228,34 @@ dho.pert.f <- function(t, y, kpass){
   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)
+                      y=c(yo, y1), 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, 
+  ydot_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
+  ydot_best$minimum  # best ydot
   
-best.sol <- rk4sys(dho.pert.f, a=0, b=2, y0=c(0, ep_best$minimum),
+  best.sol <- rk4sys(dho.pert.f, a=0, b=2, y0=c(0, ydot_best$minimum),
                      n=20)  # 20 integration stepstmax 
+}
+
+ep <- 0.5
+best.sol <- dho.pert.encap.f(ep)
+plot(best.sol$x, best.sol$y[,1], type="l")
+
+ep <- 0.05
+best.sol <- dho.pert.encap.f(ep)
+plot(best.sol$x, best.sol$y[,1], type="l")
 
 # analytical sol
 # yprime = -y
+# Same as e^x
 
 ## 4. Position of the earth and moon
 # a. Plotly
@@ -311,6 +320,10 @@ mean(sqrt(rowSums(cbind(sunearth.sol$y[,1]^2,
 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))))
+
+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
 
@@ -336,15 +349,36 @@ findZeroRelax <- function(g, x.guess, tol=1e-6, maxsteps=1e6){
   return(c(x.new, g(x.new), steps))
 }
 
-orbiter.f <- function(l) {
 T <- 365.25
 G <- 6.673e-11 # m^3 kg^-1 s^-2
 M_S <- 1.9891e30
 M_E <- 5.98e24
 x_earth <- c(-27115219762.4, 132888652547.0, 57651255508.0)/1e3 # km
-  
 r_earth <- sqrt(sum(x_earth^2))  # distance earth to sun
-  (T^2)/(4*(pi^2)) * G * (M_S/((r_earth+l)^2) + M_E/(l^2)) - r_earth
+
+orbiter.f <- function(l) {
+  tmp1 <- ((T^2)/(4*(pi^2))) * G
+  tmp2 <- (M_S/((r_earth+l)^2)) + (M_E/(l^2))
+  tmp3 <- tmp1*tmp2
 }
 
-orbiter.estimate <- findZeroRelax(orbiter.f, 1000)
+orbiter.estimate <- findZeroRelax(orbiter.f, 1000000)
+l <- orbiter.estimate[1]-r_earth
+
+x0_orbiter_tmp <- c(l, orbiter.f(l), 0)/1e3 # km
+x0_orbiter <- x0_orbiter + x0_earth
+
+y.init.orbiter <- c(x0_orbiter, v0_earth)
+sunorbit.sol <- rk4sys(f=sunearth.f, a=0, b=365.25, y0=y.init.orbiter, n=365)
+
+sunorbit.df <- data.frame(x=sunorbit.sol$y[,1],
+                          y=sunorbit.sol$y[,2],
+                          z=sunorbit.sol$y[,3])
+
+fig <- plot_ly(sunorbit.df, x = ~x, y = ~y, z = ~z, name="orbiter", 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