Finalizing code, report skeleton

This commit is contained in:
Noah L. Schrick 2023-01-19 22:53:41 -06:00
parent ec4e21488c
commit d4d8073cd0
3 changed files with 23 additions and 12 deletions

View File

@ -0,0 +1 @@
,noah,NovaArchSys,19.01.2023 22:53,file:///home/noah/.config/libreoffice/4;

Binary file not shown.

View File

@ -3,9 +3,6 @@
# 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) {
@ -17,14 +14,17 @@ computeBindingEnergy <- function(A, Z){
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)
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
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)))
# 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
@ -39,9 +39,10 @@ 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", main="Number of Nucleons on Binding Energy")
## Part D: find the value of A that gives the maximum binding energy
computeMaxBind <- function(Z, precision=NA){
if (is.na(precision)){precision <- 4}
computeMaxBind <- function(Z){
maxBind <- -Inf
A.maxBind <- -Inf
for (i in seq(Z, 3*Z)){
@ -51,7 +52,16 @@ computeMaxBind <- function(Z, precision=NA){
A.maxBind <- i
}
}
return(c(paste("Energy: ", round(maxBind, precision), "MeV"), paste("Nucleons: ", A.maxBind)))
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)){
# +1 since R indexes at 1
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")