72 lines
1.9 KiB
R
72 lines
1.9 KiB
R
# 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)
|