# Parameters
r <- 5
center <- c(0,0)

# Simulate points
nTrue <- round(exp(runif(1,6,7)))

x <- runif(nTrue,-r,r) 
y <- runif(nTrue,-sqrt(r^2 - x^2),sqrt(r^2 - x^2))


plot(c(center[1]-r,center[1]+r),c(center[2]-r,center[2]+r),type="n",xlab="x",ylab="y",main=paste("View under Microscope, True Cell Count =",nTrue))
x.microscope <- seq(center[1]-r,center[1]+r,0.005)
lines(x.microscope,sqrt(r^2 - x.microscope^2),col="red")
lines(x.microscope,-sqrt(r^2 - x.microscope^2),col="red")
points(x,y)

# Counting by sampling
nSample <- 1000
r.sample <- r * 0.25
results <- rep(0,nSample)
tol <- r.sample*0.077 #(0.075 - 0.08 gives best results)

center.sample.x <- runif(nSample,center[1]-(r - r.sample),center[1]+ (r - r.sample))
center.sample.y <- runif(nSample,-sqrt((r - r.sample)^2 - center.sample.x^2),sqrt((r -r.sample)^2 - center.sample.x^2))

par(mfrow=c(1,3))
for (i in 1:nSample) {
  # Number of cells in each sample
  index <- which(x < center.sample.x[i] + r.sample - tol & x > center.sample.x[i]  - r.sample + tol & y < center.sample.y[i] + r.sample - tol & y > center.sample.y[i] - r.sample + tol)
  nCells.est <- length(index) * r^2/r.sample^2
  results[i] <- nCells.est
  
  # Replot field of vision
  titleOfPlot <- paste("View under Microscope, Sample #",i)
  plot(c(center[1]-r,center[1]+r),c(center[2]-r,center[2]+r),type="n",xlab="x",ylab="y",main=titleOfPlot)
  x.microscope <- seq(center[1]-r,center[1]+r,0.005)
  lines(x.microscope,sqrt(r^2 - x.microscope^2),col="red")
  lines(x.microscope,-sqrt(r^2 - x.microscope^2),col="red")
  points(x[-index],y[-index])
  points(x[index],y[index],col="blue",pch=19)
  
  # Plot sampling vision
  x.microscope.sample <- seq(center.sample.x[i] - r.sample,center.sample.x[i] + r.sample,0.0001)[-1]
  x.microscope.sample <- x.microscope.sample[-length(x.microscope.sample)]
  lines(x.microscope.sample,sqrt(r.sample^2 - (x.microscope.sample - center.sample.x[i])^2) + center.sample.y[i],col="blue")
  lines(x.microscope.sample,-sqrt(r.sample^2 - (x.microscope.sample - center.sample.x[i])^2) + center.sample.y[i],col="blue")  
  
  # Plot the number of samples
  titleOfPlot <- paste("Estimate of Total Cell Count, True Cell Count =",nTrue)
  plot(seq(1,i,round(nSample*0.005)),results[seq(1,i,round(nSample*0.005))],type="l",xlab="Sample #",ylab="Estimate",main=titleOfPlot,col="blue")
  lines(1:i,cumsum(results[1:i])/(1:i),col="blue",lty=2)
  abline(h = nTrue,col="black")
  legend("bottomright",legend=c("Estimated Value","Mean of Estimated Values","True Cell Count"),lty=c(1,2,1),col=c("blue","blue","black"))
  
  # Plot the sampling distribution of the mean total cell count
  titleOfPlot <- paste("Sampling Distribution of Mean Total Cell Count")
  if (i < 2) {
    plot(c(nTrue-100,nTrue+100),c(0,1),type="n",main=titleOfPlot,xlab="Sample Mean",ylab="Density")
  } else {
    sigmaHat2 = var(results[1:i])
    x.plot = seq(mean(results[1:i])-1*sqrt(sigmaHat2),mean(results[1:i])+1*sqrt(sigmaHat2),length=nSample)
    y.plot = dnorm(x.plot,mean(results[1:i]),sqrt(sigmaHat2/i))
    plot(x.plot,y.plot,type="l",lwd=2,col="blue",main=titleOfPlot,xlab="Sample Mean",ylab="Density")
	abline(v = nTrue,col="black",lwd=2)
	legend("topleft",legend=c(paste("Sample Size =",i),"True Cell Count"),lwd=2,col=c("blue","black"))
  }
}