# 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

## 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
  }
  
  return((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))))
}

## Part B: compute Z=28 and A=58
# Nickel
bind.energy <- computeBindingEnergy(58,28)
cat(round(bind.energy,digits=2),"MeV")

## 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)

plot(carbon_df$A, carbon_df$"B/A", xlab="Number of Nucleons", 
     ylab="Binding Energy [MeV]",
     main="Number of Nucleons on Binding Energy for Carbon")

## Part D: find the value of A that gives the maximum binding energy
computeMaxBind <- function(Z){
  maxBind <- -Inf
  A.maxBind <- -Inf
  for (i in seq(Z, 3*Z)){
    bindE <- computeBindingEnergy(i, Z)/i
    if (bindE > maxBind){
      maxBind <- bindE
      A.maxBind <- i
    }
  }
  return(c(A.maxBind, maxBind))
}

## Part E: for each Z use the A that maximizes E/A
bounds <- 100
binding_df <- data.frame(nrow=bounds, ncol=2)
for (i in seq(1, bounds)){
  binding_df[i,] <- c(i, computeMaxBind(i)[2])
}

plot(binding_df, xlab="Atomic Number", ylab="B/A [MeV]", type="o",
     main="Maximum Binding Energy for Most Stable Isotopes of Each Atom")

# Z with largest binding energy
print(binding_df[which.max(unlist(binding_df[2])),], row.names=FALSE)