# Project 4 for the University of Tulsa's CS-7863 Sci-Stat Course
# Higher Order Differential Equations and Shooting Method
# Professor: Dr. McKinney, Spring 2023
# 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 
tmin <- 0
tmax <- 7
init_angle <- pi/4
init_ang_spd <- 0
g <- 9.81
L <- 1

# a. Plot
plot(pend.sol[,1],pend.sol[,2],type="l",col="blue", 
     ylim=c(-2.3,2.3), xlab="time",ylab="amplitude")
par(new=T)
lines(pend.sol[,1],pend.sol[,3],type="l",col="red")
abline(h=0)
legend("topleft", # coordinates, "topleft" etc
       c("angle","omega"), # label
       lty=c(1,1), # line 
       lwd=c(2.5,2.5), # weight
       #cex=.8, 
       bty="n", # no box
       col=c("blue","red") # color
)

# b. Overlay a plot of the small-angle approx

## 2. Use ode45 to solve the damped harmonic oscillator 
tmin <- 0
tmax <- 10
yo <- 1
ybaro <- 0

# a
m <- 1
c <- 1
k <- 2

# plot pos vs time

# b
m <- 1
c <- 2
k <- 1

# plot pos vs time

## 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
  a <- -g  # could be any f, y'' = f(t,y)
  matrix(c(v,a))
}

proj.obj <- function(v0, y0=0, tfinal){
  # minimize w.r.t. v0
  proj.sol <- ode45(projectile.f, 
                    y=c(y0, v0), t0=0, 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
v_best <- optimize(proj.obj, 
                   interval=c(1,100), #bisect-esque interval
                   tol=1e-10, 
                   y0=0, tfinal=10) # un-optimized obj params
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 

# c. Shooting method for damped oscillator with perturbation parameter
yo <- 0
y1 <- 1
ymin <- 0
ymax <- 2

pfirst <- 0.5
psec <- 0.05

# analytical sol
pthird <- 0

## 4. Position of the earth and moon
# a. Plotly
G <- 6.673e-11   # m^3 kg-1 s^-2
M_S <- 1.9891e30 # sun kg
M_E <- 5.98e24   # earth kg
M_m <- 7.32e22   # moon kg
mu_sun <- G*M_S*(86400^2)/1e9 #  km^3/days^2

# t0: jan 1, 1999 00:00:00am
x0_earth <- c(-27115219762.4, 132888652547.0, 57651255508.0)/1e3 # km
v0_earth <- c(-29794.2199907, -4924.33333333,-2135.59540741)*(86400)/1e3 
# m/s -> km/day

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

# c. keep orbiter at L2 Lagrange point for a year, ignoring the moon's effect