From 1dd9439d667f459455993ad8efecb285fd06f524 Mon Sep 17 00:00:00 2001
From: noah <noah@theschricks.com>
Date: Wed, 1 Mar 2023 18:16:14 -0600
Subject: [PATCH] L2 Lagrange orbiter

---
 Schrick-Noah_Homework-4.R | 106 +++++++++++++++++++++++++-------------
 1 file changed, 70 insertions(+), 36 deletions(-)

diff --git a/Schrick-Noah_Homework-4.R b/Schrick-Noah_Homework-4.R
index e925986..7085c18 100644
--- a/Schrick-Noah_Homework-4.R
+++ b/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,38 +216,46 @@ y1 <- 1
 tmin <- 0
 tfinal <- 2
 
-pfirst <- 0.5
-psec <- 0.05
-
-dho.pert.f <- function(t, y, kpass){
-  yn <- y[1]
-  ydot <- y[2]
-  ep <- kpass[1]
-  ypp <- (-(1+ep)*yn - ydot)/ep
-  as.matrix(c(ydot,ypp))
+dho.pert.encap.f <- function(ep){
+  dho.pert.f <- function(t, y){
+    yn <- y[1]
+    ydot <- y[2] - yn
+    #ep <- 0.5
+    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), 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
+  ydot_best <- optimize(proj.obj, 
+                        interval=c(1,100), #bisect-esque interval
+                        tol=1e-10, 
+                        y0=0, tfinal=2) # un-optimized obj params
+  ydot_best$minimum  # best ydot
+  
+  best.sol <- rk4sys(dho.pert.f, a=0, b=2, y0=c(0, ydot_best$minimum),
+                     n=20)  # 20 integration stepstmax 
 }
 
-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
-}
+ep <- 0.5
+best.sol <- dho.pert.encap.f(ep)
+plot(best.sol$x, best.sol$y[,1], type="l")
 
-# 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 
+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))
 }
 
+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
+
 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
+  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