# Monte Carlo Simulation of the Tomato Salad Problem
# Author: J. A. Österreicher
# Description: This script simulates dispersoid distributions in TEM micrographs using a Monte Carlo approach.

library(rgl)
library(ggplot2)
library(ggforce)

rgl.clear()

# Sigma/Mu values of lognormal distribution ensuring r_mean = 50
# σ = 0.2 -> μ = 3.892
# σ = 0.5 -> μ = 3.787
# σ = 0.6 -> μ = 3.732
# σ = 0.7 -> μ = 3.667
# σ = 0.8 -> μ = 3.592
# σ = 1.0 -> μ = 3.412


# Simulation parameters
volumeDisp <- 0.01 	# volume fraction of dispersoids (particles)
mu <- 3.787 		# parameter for lognormal distribution
sigma <- 0.5 		# parameter for lognormal distribution
sizeFactor <- 1 	# optional: parameter for multiplying lognormal distribution

# TEM micrograph dimensions (in nm)
widthTEM <- 5000 	# width of TEM micrograph in nm
heightTEM <- 5000 	# height of TEM micrograph in nm
thicknesses <- c(50,100,150)	# thicknesses of TEM foils to be simulated, e.g. c(25,50,75,100,150,200,300,400,500,700,1000)

setwd("J:\\Projekte\\3700094700_AMALFI_2.1\\3.1_Projectwork\\Dispersoide\\dispersoidsimulator\\2025\\")

len_thick <- length(thicknesses)

areaTEM <- widthTEM * heightTEM
totalHeight <- 20000 # height of simulated volume; should be several times the width or height of the TEM image
volumeDispNm <- totalHeight * areaTEM * volumeDisp 	# total volume of dispersoids in nm³


# # # # # # # # # # # # # # # # #
iter <- 30 # number of simulations (=micrographs) per foil thickness
# # # # # # # # # # # # # # # # #

results <- matrix(NA, nrow=iter, ncol=len_thick)
resultsN <- matrix(NA, nrow=iter, ncol=len_thick)

for(a in 1:len_thick){
  for(i in 1:iter){
    
    foilThickness <- thicknesses[a]
    rm(dispersoids)

    sizeCurrent <- rlnorm(1, mu, sigma) * sizeFactor # radius of first dispersoid

    dispersoids <- data.frame(dispRadius=sizeCurrent,xPos=runif(1,min=0+sizeCurrent,max=widthTEM-sizeCurrent),yPos=runif(1,min=0+sizeCurrent,max=heightTEM-sizeCurrent),zPos=runif(1,min=0+sizeCurrent,max=totalHeight-sizeCurrent),   stringsAsFactors = FALSE) # create dataframe with first dispersoid

    volumeDispNmCurrent <- 0

    while(volumeDispNmCurrent <= volumeDispNm) {		# append data frame by one dispersoid at a time
      sizeCurrent <- rlnorm(1, mu, sigma) * sizeFactor # radius of current dispersoid
      New <- c(dispRadius=sizeCurrent,xPos=runif(1,min=0+sizeCurrent,max=widthTEM-sizeCurrent),yPos=runif(1,min=0+sizeCurrent,max=heightTEM-sizeCurrent),zPos=runif(1,min=0+sizeCurrent,max=totalHeight-sizeCurrent))  
      dispersoids = rbind(dispersoids,New)
      volumeDispNmCurrent <- volumeDispNmCurrent + sizeCurrent^3 * 4 * pi / 3
    }
    
    
    attach(dispersoids)
    
    randZ <- round(runif(1,min=1000,max=totalHeight-1000),0) # create position for TEM foil
    
    foil <- subset(dispersoids, zPos+dispRadius > randZ & zPos-dispRadius < randZ + foilThickness , select=c(dispRadius,xPos,yPos,zPos)) # identify all dispersoids that have common volume with the foil
    
	detach(dispersoids)
	
	attach(foil)
	foilCut1 <- subset(foil, (zPos < randZ) | (zPos > (randZ + foilThickness)) , select=c(dispRadius,xPos,yPos,zPos)) # dispersoids that lie outside the foil but are cut
	foilCut2 <- subset(foil, ((zPos >= randZ & zPos < (randZ + dispRadius)) & zPos <= randZ+foilThickness) | (zPos > (randZ+foilThickness-dispRadius) & zPos <= (randZ + foilThickness) & (zPos >= randZ)) , select=c(dispRadius,xPos,yPos,zPos)) # dispersoids that lie inside the foil but are cut
	foilUnCut <- subset(foil, zPos > randZ + dispRadius & zPos < randZ + foilThickness - dispRadius , select=c(dispRadius,xPos,yPos,zPos)) # dispersoids entirely within the foil
    
    detach(foil)
    
    attach(foilCut1)
    
    foilCut1lower <- subset(foilCut1, zPos < randZ,  select=c(dispRadius,xPos,yPos,zPos)) 
    foilCut1upper <- subset(foilCut1, zPos > randZ + foilThickness,  select=c(dispRadius,xPos,yPos,zPos)) 
    
    detach(foilCut1)
    
    # Calculate apparent radii
    foilUnCut$appR <- foilUnCut$dispRadius
    foilCut2$appR <- foilCut2$dispRadius
    
    foilCut1lower$appR <- round(sqrt(foilCut1lower$dispRadius^2-(randZ-foilCut1lower$zPos)^2),0) # calculate apparent radii, round to whole nm
    foilCut1upper$appR <- round(sqrt(foilCut1upper$dispRadius^2-(foilCut1upper$zPos-randZ-foilThickness)^2),0) # calculate apparent radii, round to whole nm
    foilCut1final <- rbind(foilCut1upper,foilCut1lower)
    foilFinal <- rbind(foilCut1final, foilCut2)
    foilFinal <- rbind(foilFinal, foilUnCut)
    results[i,a] <- mean(foilFinal$appR)
    if(i%%10==0) cat("|")	# this is to show progress in the console
    flush.console()
    
    stopifnot(nrow(foil)==nrow(foilFinal))  # if something goes wrong
    resultsN[i,a] <- nrow(foilFinal)
    
    }
    cat(foilThickness, " complete\n")
    flush.console()
  }
  
  
cat("\napparent mean r\n")		# print some results
for(a in 1:len_thick) cat(mean(results[,a],na.rm = TRUE),"\n")
cat("sd\n")
for(a in 1:len_thick) cat(sd(results[,a],na.rm = TRUE),"\n")
cat("number of particles in image r\n")
for(a in 1:len_thick) cat(mean(resultsN[,a],na.rm = TRUE),"\n")
cat("sd of number of particles in image r\n")
for(a in 1:len_thick) cat(sd(resultsN[,a],na.rm = TRUE),"\n")

# Now create and print the histograms of dispersoid sizes
#cat("\nHistogram of dispersoid sizes:\n")
hist(dispersoids$dispRadius, main="Histogram of Dispersoid Sizes", xlab="Dispersoid Radius [nm]", col="lightblue", border="black", breaks=30)

# draw dispersoids in 3D

rgl.spheres(foilCut1$xPos,foilCut1$yPos,foilCut1$zPos,r=foilCut1$dispRadius,color="#e66101") # draw cut dispersoids that appear smaller in TEM
rgl.spheres(foilCut2$xPos,foilCut2$yPos,foilCut2$zPos,r=foilCut2$dispRadius,color="#b2abd2") # draw cut dispersoids that appear with their real diameter in TEM
if(nrow(foilUnCut)) rgl.spheres(foilUnCut$xPos,foilUnCut$yPos,foilUnCut$zPos,r=foilUnCut$dispRadius,color="magenta") # draw uncut dispersoids within the foil

# draw TEM foil in 3D

rgl.quads(x=c(0,5000,5000,0),y=c(0,0,5000,5000),z=c(randZ,randZ,randZ,randZ), alpha=0.5)
rgl.quads(x=c(0,5000,5000,0),y=c(0,0,5000,5000),z=c(randZ+foilThickness,randZ+foilThickness,randZ+foilThickness,randZ+foilThickness),alpha=0.5)

axis3d('z', pos=c( 0 , 0, NA ), col = "black")
title3d(zlab = "[nm]") 
rgl.bg(color="white")


# Plot TEM image with axes in µm
p6 <- ggplot() + 
  theme(plot.margin = margin(0, 0, 0, 0, "pt")) +
  coord_fixed(ratio = 1, xlim = c(0, 5), ylim = c(0, 5)) +  # Adjusted limits to µm scale
  
  geom_circle(data = foilUnCut, 
              aes(x0 = xPos / 1000, y0 = yPos / 1000, r = dispRadius / 1000), 
              color = "magenta", fill = "magenta") +
  
  geom_circle(data = foilCut1final, 
              aes(x0 = xPos / 1000, y0 = yPos / 1000, r = dispRadius / 1000), 
              color = "#e66101", lty = "dashed", fill = "#fdb863") +
			  
  geom_circle(data = foilCut2, 
              aes(x0 = xPos / 1000, y0 = yPos / 1000, r = dispRadius / 1000), 
              color = "#b2abd2", fill = "#b2abd2") +

  geom_circle(data = foilCut1final, 
              aes(x0 = xPos / 1000, y0 = yPos / 1000, r = appR / 1000), 
              color = "#e66101", fill = "#e66101") +

  scale_x_continuous(limits = c(0, 5), labels = scales::number_format(accuracy = 1)) +  
  scale_y_continuous(limits = c(0, 5), labels = scales::number_format(accuracy = 1)) +  
  
  labs(x = "x Position [µm]", y = "y Position [µm]") +  # Updated axis labels

  theme(
    text = element_text(family = "Times New Roman", size = 24),  # Larger overall text
    axis.title = element_text(size = 20),  # Larger axis labels
    axis.text = element_text(size = 20),  # Larger tick labels
    plot.title = element_text(size = 20)  # Larger title
  )

# Save the plot as a PDF
ggsave("TEM_100.pdf", plot = p6, device = cairo_pdf)


# plot b/w TEM image with scalebar

# p6 <- ggplot() + 
# coord_fixed(ratio = 1, xlim = c(0, 5000), ylim = c(0, 5000)) + 
#  geom_circle(data=foilUnCut, aes(x0=xPos, y0=yPos, r=dispRadius), color="black", fill="black") +
#  geom_circle(data=foilCut2, aes(x0=xPos, y0=yPos, r=dispRadius), color="black", fill="black") +
#  geom_circle(data=foilCut1final, aes(x0=xPos, y0=yPos, r=appR), color="black", fill="black") +
#  labs(x = "x-Position [nm]", y = "y-Position [nm]") + 
#  ggtitle(paste("TEM micrograph (foil thickness", foilThickness, "nm)")) +
#  theme(
#    panel.grid.major = element_blank(),     # Remove major grid lines
#    panel.grid.minor = element_blank(),     # Remove minor grid lines
#    axis.line = element_blank(),            # Remove axis lines
#    panel.background = element_rect(fill = "grey90"),  # Grey background
#    plot.background = element_blank()
#  ) +
#  # Add the scale bar as a black rectangle
#  annotate("rect", xmin = 3800, xmax = 4800, ymin = 200, ymax = 250, 
#           fill = "black", color = "black") +
#  annotate("text", x = 4300, y = 400, label = "1 µm", size = 5, hjust = 0.5)



# Mean diameter with error bars (SD)
mean_diameters <- numeric(len_thick)
sd_diameters <- numeric(len_thick)

# Calculate mean and SD of diameters for each foil thickness
for(a in 1:len_thick) {
  mean_diameters[a] <- mean(results[, a], na.rm = TRUE) * 2  # Diameter is 2 times the radius
  sd_diameters[a] <- sd(results[, a], na.rm = TRUE) * 2  # Standard deviation for diameter
}

# Create a data frame for plotting
plot_data <- data.frame(
  Thickness = thicknesses,
  MeanDiameter = mean_diameters,  # Mean diameter
  SD = sd_diameters  # Standard deviation for diameter
)

# Plot observed mean diameter over foil thickness

plot_obj <- ggplot(plot_data, aes(x = Thickness, y = MeanDiameter)) +
  geom_line(color = "steelblue") +  # Line for mean diameter with color
  geom_point(color = "firebrick") +  # Points for mean diameter with color
  geom_errorbar(aes(ymin = MeanDiameter - SD, ymax = MeanDiameter + SD), width = 5, color = "darkgray") +  # Error bars with color
  labs(
    title = "Mean Diameter of Dispersoids vs Foil Thickness",
    x = "Foil Thickness (nm)",
    y = "Observed Mean Diameter (nm)"
  ) +
  scale_y_continuous(
    limits = c(0, NA),  # Set y-axis to start at 0
    breaks = seq(0, max(plot_data$MeanDiameter, na.rm = TRUE), by = 10)  # Y-axis ticks every 10 nm
  ) +
  theme_minimal() +
  theme(
    axis.line = element_line(color = "black"),  # Add black lines to the axes
    axis.ticks = element_line(color = "black"),  # Add ticks to the axes
    axis.text = element_text(color = "black"),  # Axis labels in black
    plot.title = element_text(hjust = 0.5, color = "darkblue", size = 14, face = "bold"),  # Centered title with color and size
    axis.title = element_text(size = 12)  # Axis titles size
  )
# Explicitly print the plots
dev.new()
print(plot_obj)
dev.new()
print(p6)