#-------------------------------------------------------------------------
# Dati delle orbite dei pianeti del sistema solare
# https://nssdc.gsfc.nasa.gov/planetary/factsheet/planet_table_ratio.html
# semiasse maggiore in unità di quello terrestre
#   Versione con log-plot
#
# GdA marzo 2023
#--------------------------------------------------------------------------

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

# periodo di rivoluzione in unità dell'anno terrestre 
T <- c(0.241, 0.615, 1, 1.88, 11.9, 29.4, 83.7, 163.7)
# 'distanza' dal Sole (-> semiasse maggiore) in funzione dell'U.A.
d <- c(0.387, 0.723,  1 , 1.52, 5.20, 9.57, 19.17, 30.18)
n = length(T)

# limitiamo il nr di pianeti da analizzare (min 2, max 8)
n1 = as.integer(readline("Numero di pianeti (2-8; default 8): "))
if ( is.na(n1) ) n1 = 8
if (n1 > 1 && n1 <= 8) {
    n <- n1
    T <- T[1:n]
    d <- d[1:n]
    cat(sprintf(" => Consideriamo soli i primi %d pianeti \n", n))
} else {
    cat(sprintf(" %d non è valido \n", n1))
    stop()
}

# primo plot, d vs T
plot(T, d, pch=19, xlab='T (y)', ylab='d (R_T)',
     col='blue', main="Pianeti del sistema solare: 'distanza' vs periodo")
grid()
pausa()

# linearizzazione, confidando nella Terza Legge di Keplero
plot(T^2,d^3, pch=19, xlab='T^2 (y^2)', ylab='d^3 (R_T^3)', col='blue',
     main='Pianeti del sistema solare -- linearizzazione mediante Keplero')
grid()
abline(0, 1, col='cyan')
pausa()

# linearizzazione mediante scale logaritmiche ('doppio-log', o 'log-log')
plot(T,d, pch=19, col='blue', log='xy',  xlab='T (y)', ylab='d (R_T)',
     main='Pianeti del sistema solare -- scale logaritmiche')
grid()

# sovrapponiamo quanto atteso dalla Terza Legge di Keplero
Ti <- seq(0.5*T[1], 1.5*T[n], len=100)
points(Ti, Ti^(2/3), ty='l', lty=2, col='cyan')
cat(" -> la linea dà l'andamento secondo Keplero, ovvero d prop. a T^(3/2)\n") 

pausa()

# ricaviamoci il rapporto delle potenze 
# 1) usando primo e ultimo punto:
points(T[1], d[1], pch=1, col='magenta', cex=3)
points(T[n], d[n], pch=1, col='magenta', cex=3)
m1 = ( log(d[n]) - log(d[1]) ) / ( log(T[n]) - log(T[1]) )
cat(sprintf("\n Potenza empirica dell'andamento di d in funzione di T: %.4f\n", m1))
cat(" --> Ricavata facendo uso del primo e dell'ultimo punto\n")
cat("     (inutile fare la retta che passa meglio fra i punti,")
cat(" dati gli scarti esigui)\n")
points(Ti, Ti^m1, ty='l', lty=3, col='magenta')

# volendo, si può anche trovare la slope della migliore
# retta fra i punti, ma gli scostamenti dei punti dalle
#  due rette trovate sono talmente esigui che non ne vale la pena