# Project 1 for the University of Tulsa's CS-7863 Sci-Stat Course
# Nuclear Binding Energy
# Professor: Dr. McKinney, Spring 2023
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Part A: compute binding energy
computeBindingEnergy <- function(A, Z){
  if (A %% 2 == 1) {
    spin <- 0
  }else if (A %% 2 == 0 && Z %% 2 == 0){
    spin <- 12
  }
  else {
    spin <- -12
  }
  
  B <- (15.67*A) - (17.23*(A^(2/3))) - ((0.75)*((Z^2)/(A^(1/3)))) - ((93.2)*(((A-(2*Z))^2)/A)) + (spin/(A^(1/2)))
  
  return(B)
}

## Part B: compute Z=28 and A=58
computeBindingEnergy(58,28)
(15.67*58) - (17.23*(58^(2/3))) - ((0.75)*((28^2)/(58^(1/3)))) - ((93.2)*(((58-(2*28))^2)/58)) + (12/(58^(1/2)))

## Part C: table with two columns: nucleons and binding energy per nucleon
Z=6 # Carbon
carbon_df <- data.frame(nrow=2*Z, ncol=2)
for (i in seq(Z, 3*Z)){
  # +1 since R indexes at 1
  carbon_df[i-Z+1,] <- rbind(i, computeBindingEnergy(i, Z)/i)
}
rownames(carbon_df) <- NULL
colnames(carbon_df) <- c("A", "B/A")
print(carbon_df, row.names=FALSE)
# Get max binding energy row
print(carbon_df[which.max(carbon_df$"B/A"),], row.names=FALSE)

## Part D: find the value of A that gives the maximum binding energy

## Part E: for each Z use the A that maximizes E/A