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

Initial Work for L2 Lagrange orbiter

Browse Source
This commit is contained in:
Noah L. Schrick 2023-03-01 15:38:29 -06:00
parent 511ce16679
commit dd599c8b66
1 changed files with 34 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
Schrick-Noah_Homework-4.R
View File

@ -4,10 +4,10 @@
# Noah L. Schrick - 1492657 # Noah L. Schrick - 1492657
## 1. Transform the second order ODE into a system of two first order ODEs ## 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 # use ode45 or rk4sys to solve for the motion of the plane pendulum
if (!require("pracma")) install.packages("pracma") if (!require("pracma")) install.packages("pracma")
library(pracma) library(pracma)
tmin <- 0 tmin <- 0
tmax <- 7 tmax <- 7
theta0 <- pi/4 theta0 <- pi/4
@ -94,7 +94,7 @@ kvar <- 2
dho.f <- function(t, y, kpass){ dho.f <- function(t, y, kpass){
yn <- y[1] yn <- y[1]
ydot <- y[2] ydot <- y[2]
# c k m # c k m
ypp <- (-kpass[2]*yn - kpass[3]*ydot)/kpass[1] ypp <- (-kpass[2]*yn - kpass[3]*ydot)/kpass[1]
as.matrix(c(ydot,ypp)) as.matrix(c(ydot,ypp))
} }
@ -170,7 +170,7 @@ plot(ivp.sol$t, ivp.sol$y[,1] ,type="o",col="blue",
par(new=T) par(new=T)
lines(ivp.sol.rk4$x, ivp.sol.rk4$y[,1] ,type="o",col="red", lines(ivp.sol.rk4$x, ivp.sol.rk4$y[,1] ,type="o",col="red",
xlab="time",ylab="height") xlab="time",ylab="height")
legend("topright", # coordinates, "topleft" etc legend("topright",
c("ODE45, default","RK4Sys, 100 steps"), # label c("ODE45, default","RK4Sys, 100 steps"), # label
lty=c(1,1), # line lty=c(1,1), # line
lwd=c(2.5,2.5), # weight lwd=c(2.5,2.5), # weight
@ -222,7 +222,6 @@ dho.pert.f <- function(t, y, kpass){
yn <- y[1] yn <- y[1]
ydot <- y[2] ydot <- y[2]
ep <- kpass[1] ep <- kpass[1]
# c m
ypp <- (-(1+ep)*yn - ydot)/ep ypp <- (-(1+ep)*yn - ydot)/ep
as.matrix(c(ydot,ypp)) as.matrix(c(ydot,ypp))
} }
@ -270,7 +269,7 @@ sunearth.f <- function(t,y){
G <- 6.673e-11 # m^3 kg^-1 s^-2 G <- 6.673e-11 # m^3 kg^-1 s^-2
M_S <- 1.9891e30 M_S <- 1.9891e30
M_E <- 5.98e24 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 mu_sun <- G*M_S*(86400^2)/1e9 # km^3/days^2
x_earth <- y[1:3] # xyz position of earth wrt sun 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[,2] - sunmoon.sol$y[,2])^2,
(sunearth.sol$y[,3] - sunmoon.sol$y[,3])^2)))) (sunearth.sol$y[,3] - sunmoon.sol$y[,3])^2))))
# Not correct: giving 60m km - should be ~380k # Not correct: giving 60m km - should be ~380k
# Why: Don't think M_m is used in function
# b. Find eclipses # b. Find eclipses
norm_vec <- function(x) sqrt(rowSums(cbind((x^2)))) norm_vec <- function(x) sqrt(rowSums(cbind((x^2))))
eclipses <- norm_vec(sunearth.sol$y) %*% norm_vec(sunmoon.sol$y) 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 # 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)