#---------------------------------------------------- -------------------------------
# Volume di un generico solido di rotazione -> uovo di polistirolo
# usa il pacchetto magick:
# https://cran.r-project.org/web/packages/magick/vignettes/intro.html#Read_and_write
#
# Gda, aprile 2023
#------------------------------------------------------------------------------------


library(magick)

pausa <- function() { cat ("\n >> press Enter to continue\n"); scan() }

# funzione che disegna un cerchio (da notare '...' fra gli argomenti!)
cerchio <- function(C=c(0,0), r=1, ...) {
  th <- seq(0, 2*pi, by=0.02)
  points(C[1] + r*cos(th), C[2] + r*sin(th), ty='l', ...)
}

# importa l'immagine jpeg, la trasforma in 'raster' e la mostra in un plot
uovo   <- image_read('uovo_polistirolo.jpeg')
uovo.r <- as.raster(uovo)
plot(uovo.r)
grid(7)
cat("\n Proprietà file 'uovo': \n")
print(uovo)
cat("\n Proprietà file 'uovo.r': \n")
str(uovo.r)

# dimensioni fisiche dell'uovo 
a  <- 12.9  # cm   'asse maggiore'  
b  <- 8.9   # cm   'asse minore' 
RM <- b/2   # raggio della sezione maggiore  (in cm) 

N  <- 50    # numero di dischi con cui affettare l'uovo

cat("\n clicca su 2 punti per selezionare x.min e x.max\n")
p <- locator(2)
xl <- c(p$x[1], p$x[2])
cat(" clicca su 2 punti per selezionare y.min e y.max\n")
p <- locator(2)
yl <- c(p$y[1], p$y[2])

# traccia il rettangolo di contorno 
points(xl, rep(yl[1],2), ty='l', lwd=1)
points(xl, rep(yl[2],2), ty='l', lwd=1)
points(rep(xl[1],2), yl, ty='l', lwd=1)
points(rep(xl[2],2), yl, ty='l', lwd=1)

#  punto di intersezione degli assi
c <- c(mean(p$x) ,mean(yl))   

# traccia gli assi 
points(xl, rep(c[2],2), ty='l', lwd=1, lty=2)
points(rep(c[1], 2), xl, ty='l', lwd=1, lty=2)

# traccia il cerchio che ha centro sull'intersezione degli assi 
# e raggio pari al semiasse minore 
R <- (yl[2]-yl[1])/2   
cerchio(c, R, col='yellow', lty=2)
V.sfera <-  4/3*pi*RM^3
cat(sprintf("\n => Volume della sfera centrata nell'incrocio degli assi: %.0f cm^3\n",
            V.sfera))

# riquadro giallo per terminare il tracciamento dei punti
points(xl[2], yl[2], pch=15, cex=10, col='yellow')
# primo e ultimo punto già determinati 
points(xl[1], c[2], cex=1, col='yellow')
points(xl[2], c[2], cex=1, col='yellow')

cat("\n clicca lungo il contorno superiore dell'uovo \n -> nel quadrato giallo per terminare\n")
cat("    (primo e ultimo punto già determinati e plottati)\n")

n <- 1
# primo punto (già definito))
ux <-  xl[1]     
uy <-  c[2]
# punti intermedi, rilevati cliccando
while(1) {
   p = locator(1)
   if( abs(p$x[1] - xl[2]) < 20 &&  abs(p$y[1] - yl[2]) < 50 ) break
   if(p$x[1] <= ux[n] ||  p$x[1] >= xl[2]) {
       cat("   Att: le x devono essere crescenti! \n")       
   } else {
     points(p, cex=1, col='yellow')
     n <- n + 1
     ux[n] <- p$x[1]
     uy[n] <- p$y[1]
   }
}
n <- n + 1
# ultimo punto (già definito))
ux[n] <- xl[2]
uy[n] <- c[2]
# traccia la spezzata lungo i punti
points(ux, uy, cex=1, ty='l', , lwd=2, col='yellow')

pausa()


# funzione per interpolare fra i punti del contorno superiore dell'uovo
y.fun <- approxfun(ux, uy)

# Plottiamo i rettangolini che rappresentano i dischetti
dx <- (xl[2]-xl[1]) / N          # spessore dei dischetti
xi <- seq(xl[1], xl[2], by=dx)   # punti di separazione fra un dischetto e l'altro

for(i in 1:N) {
   # raggio dell'i-mo dischetto
   r.m <- ( y.fun(xi[i]) + y.fun(xi[i+1])) / 2 - c[2]
   # rappresentazione grafica 
   points( rep(xi[i], 2), c(c[2]+r.m, c[2]-r.m), ty='l', col='cyan' )
   points( rep(xi[i+1], 2), c(c[2]+r.m, c[2]-r.m), ty='l', col='cyan' )
   points( c(xi[i],xi[i+1]), rep(c[2]+r.m, 2), ty='l', col='cyan' )
   points( c(xi[i],xi[i+1]), rep(c[2]-r.m, 2), ty='l', col='cyan' ) 
}

cat("\n Passiamo al calcolo del volume \n")

#-------------------------------------------------------------------
# esprimiamo i punti del contorno superiore (ux e uy) in cm
# per poter valutare il volume

# 1) convertiamo l'ascissa da pixel a cm 
# xl[1] -> 0;  xl[2] = a
x <- (ux - xl[1]) / (xl[2]-xl[1]) *  a            #  coordinata assiale, lungo x

# 2) convertiamo l'ordinata espressa da pixel a cm, a partire da y del centro 
# c[2] -> 0;   yl[1] -> - RM;   yl[2] -> + RM;
# la chiamiamo 'r' perché reppresenta il raggio del dischetto in 'x'
r <- (uy - mean(yl)) / (yl[2]-mean(yl)) * RM       #raggio in funzione di x

# in analogia a y.fun, prepariamo la funzione per interpolare i valori di r

r.fun <- approxfun(x, r)
dx.cm <- a/N                       # chiamiamo 'dx' 'dx.cm' sia per ricordarci
                                   # che sono in cm, sia per non sovrascrivere dx 
xi.cm <- seq(0, a, by=dx.cm)       # idem per 'xi', che diventa ora 'xi.cm'


V = 0
for (i in 1:N) {
    xc <- mean(xi.cm[i], xi.cm[i+1])
    dV <- pi*r.fun(xc)^2 * dx.cm
    V <- V + dV 
}
cat(sprintf("Volume %.0f cm^3\n", V))
ff <- V / V.sfera
cat(sprintf("\nFattore di forma (Volume su volume sfera) %.2f\n", ff))