#-------------------------------------------------------- 
# script per tracciare andamenti di potenza linearizzati
# mediante scale log-log
#  
#   y = a * x^b
#
# GdA marzo 2025
#--------------------------------------------------------

# plot log-log vuoto con range prefissati ---------------
plot(NULL, xlim=c(0.1, 10), ylim=c(1, 100), log='xy',
     xlab='x', ylab='y', xaxs='i', yaxs='i',)
# Nota: non sembra ci sia un modo per avere le decadi in x e y
#       delle stesse dimensioni: -> provare a 'deformare' il plot 

# divisioni fra le tacche principali
abline(v=seq(0.10,0.50,by=0.01), lwd=0.1)
abline(v=seq(0.52,1,by=0.02), lwd=0.1)
abline(v=10*seq(0.10,0.50,by=0.01), lwd=0.1)
abline(v=10*seq(0.52,1,by=0.02), lwd=0.1)
abline(v=seq(0.1,1,by=0.1), lwd=0.15)
abline(v=10*seq(0.1,1,by=0.1), lwd=0.15)

abline(h=10*seq(0.10,0.50,by=0.01), lwd=0.1)
abline(h=10*seq(0.52,1,by=0.02), lwd=0.1)
abline(h=10*10*seq(0.10,0.50,by=0.01), lwd=0.1)
abline(h=10*10*seq(0.52,1,by=0.02), lwd=0.1)
abline(h=10*seq(0.1,1,by=0.1), lwd=0.15)
abline(h=10*10*seq(0.1,1,by=0.1), lwd=0.15)
#-------------------------------------------------------


cat(sprintf("\n Clicca su due punti (coincidenti per uscire)\n"))

xp <- seq(0.1, 10, len=11)  #  ascisse per tracciare la curva

while(1) {
    # leggiamo i punti su cui si è cliccato
    p <- locator(2)

    # se punti circa coincidenti usciamo
    if( abs(log(p$x[2]) - log(p$x[1])) < 0.05 &&        
        abs(log(p$y[2]) - log(p$y[1])) - 0.05 ) break

    # plottiamo i punti letti 
    points(p)

    # valutiamo i parametri dell'andamento di potenza 
    b <- log( p$y[2]/p$y[1] ) / log( p$x[2]/p$x[1] )
    a <- exp( log(p$y[1]) - b*log(p$x[1]) )

    cat(sprintf("P1: (%.2f , %.2f) \n", p$x[1], p$y[1]))
    cat(sprintf("P2: (%.2f , %.2f) \n", p$x[2], p$y[2]))
    cat(sprintf("Legge di potenza:  a = %.3f, b=%.3f \n\n", a, b))  

    # tracciamo le 'rette' 
    points(xp, a*xp^b, ty='l', col=sample(rainbow(14))[1], lwd=2)
}