From dd599c8b669d0477fbed09f3d3071c67eccc087c Mon Sep 17 00:00:00 2001
From: noah <noah@theschricks.com>
Date: Wed, 1 Mar 2023 15:38:29 -0600
Subject: [PATCH] Initial Work for L2 Lagrange orbiter

---
 Schrick-Noah_Homework-4.R | 39 ++++++++++++++++++++++++++++++++++-----
 1 file changed, 34 insertions(+), 5 deletions(-)

diff --git a/Schrick-Noah_Homework-4.R b/Schrick-Noah_Homework-4.R
index e77f840..e925986 100644
--- a/Schrick-Noah_Homework-4.R
+++ b/Schrick-Noah_Homework-4.R
@@ -4,10 +4,10 @@
 # Noah L. Schrick - 1492657
 
 ## 1. Transform the second order ODE into a system of two first order ODEs
-
 # use ode45 or rk4sys to solve for the motion of the plane pendulum 
 if (!require("pracma")) install.packages("pracma")
 library(pracma)
+
 tmin <- 0
 tmax <- 7
 theta0 <- pi/4
@@ -94,7 +94,7 @@ kvar <- 2
 dho.f <- function(t, y, kpass){
   yn <- y[1]
   ydot <- y[2]
-  #         c         k         m
+  #         c               k             m
   ypp <- (-kpass[2]*yn - kpass[3]*ydot)/kpass[1]
   as.matrix(c(ydot,ypp))
 }
@@ -170,7 +170,7 @@ plot(ivp.sol$t, ivp.sol$y[,1] ,type="o",col="blue",
 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
+legend("topright",
        c("ODE45, default","RK4Sys, 100 steps"), # label
        lty=c(1,1), # line 
        lwd=c(2.5,2.5), # weight
@@ -222,7 +222,6 @@ 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))
 }
@@ -270,7 +269,7 @@ 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
+  M_m <- 7.32e22  
   mu_sun <- G*M_S*(86400^2)/1e9  # km^3/days^2
   
   x_earth <- y[1:3]  # xyz position of earth wrt sun
@@ -313,9 +312,39 @@ 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
+# Why: Don't think M_m is used in function
 
 # 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
+findZeroRelax <- function(g, x.guess, tol=1e-6, maxsteps=1e6){
+  steps <- 0
+  x.old <- x.guess
+  isConverged <- FALSE
+  while (!isConverged){
+    x.new <- g(x.old)
+    if(steps >= maxsteps || (abs(x.new-x.old)<tol)){
+      isConverged <- TRUE
+    }
+    else{
+      x.old <- x.new
+    }
+    steps <- steps + 1
+  }
+  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.estimate <- findZeroRelax(orbiter.f, 1000)