# Voglio tracciare il grafico delle ore di luce nei giorni del 2011 e del
# 2012 in varie città. Traccio i grafici rilevando le ore di luce ogni
# 30 giorni. Tengo conto che le date corrispondenti sono:
# anno non bisestile (come il 2011), ogni 30 giorni
# 1/1, 31/1, 2/3, 1/4, 1/5, 31/5, 30/6, 30/7, 29/8, 28/9, 28/10, 27/11, 27/12
# (31+28+31+30+31+30+31+31+30+31+30+31)
# anno bisestile (come il 2012), ogni 30 giorni
# 1/1, 31/1, 1/3, 31/3, 30/4, 30/5, 29/6, 29/7, 28/8, 27/9, 27/10, 26/11, 26/12
# (31+29+31+30+31+30+31+31+30+31+30+31)
# Cerco in WolframAlpha le ore di luce a Genova ogni 30 giorni nel 2011 e 2012
# con Genoa,sun,1/1/2011, Geona,sun,31/1/2011,..., Genoa,sun,27/12/2011,
# Genoa,sun,1/1/2012,..., Genoa,sun,26/12/2012, e trovo per i due anni:
# (potrei trovarle anche per gli altri giorni, ma supponiamo che non sia possibile:
# potrei usare anche altre fonti; ad es.: http://www.aopa.it/sole.asp)
LuceGenova <- c(
8+55/60,9+47/60,11+12/60,12+45/60,14+12/60,15+17/60,15+30/60,
14+44/60,13+25/60,11+54/60,10+24/60,9+13/60,8+51/60,
8+54/60,9+47/60,11+11/60,12+44/60,14+12/60,15+16/60,15+30/60,
14+44/60,13+24/60,11+54/60,10+25/60,9+14/60,8+51/60)
# I giorni, in questi due anni, sono;
giorni1 <-seq(1,365,30); giorni2 <- seq(1,366,30)
giorni <- c(giorni1, giorni2+365)
# Il grafico delle ore di luce a Genova usando una spline cubica
# (f.polin. di 3° grado passante per i punti, ivi derivabile fino al 2° ordine)
dev.new(width=6,height=4)
OreLuce <- LuceGenova; plot(giorni,OreLuce)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
lines(spline(giorni,OreLuce),col="brown")
abline(v=365.5,h=12,lty=2)
# Altra citta', recuperando i dati in modo analogo:
LuceMelbourne <- c(
14+44/60,14+3/60,12+55/60,11+41/60,10+33/60,9+45/60,9+34/60,
10+8/60,11+9/60,12+20/60,13+31/60,14+29/60,14+47/60,
14+44/60,14+3/60,12+55/60,11+42/60,10+34/60,9+45/60,9+34/60,
10+8/60,11+8/60,12+19/60,13+30/60,14+28/60,14+47/60)
OreLuce <- LuceMelbourne; points(giorni,OreLuce)
lines(spline(giorni,OreLuce),col="red")

# Invece che le spline, in una scuola media inferiore avrei potuto usare:
...
lines(giorni,OreLuce,lty=2,col="brown")
...
lines(giorni,OreLuce,lty=2,col="red")
# Stima della durata minima e massima del dì:
min(LuceGenova); max(LuceGenova);min(LuceMelbourne); max(LuceMelbourne)
# 8.85  15.5    9.566667  14.78333
.85*60;.5*60;.566667*60;.78333*60
# 51  30    34.00002  46.9998
# 8 h 51 min; 15 h 30 min; 9 h 34 min; 14 h 47 min

# Vediamo la durata minima e massima del dì a Genova con le spline:
OreLuce <- LuceGenova
min(spline(giorni,OreLuce)$y); max(spline(giorni,OreLuce)$y)
#   8.838997   15.54231
0.838997*60; 0.54231*60
#  50.33982   32.5386
# 8:50 e 15:33 sono le durate minima e massima del dì a Genova
# Volendo calcolare la funzione nei vari giorni posso, invece che
# tracciare solo la spline, esprimerla esplicitamente:
# function(x) splinefun(giorni,OreLuce)(x,n) con n  0, 1 o 2
# è la funzione polinomiale, o la sua derivata 1ª o 2ª
f <- function(x) splinefun(giorni,OreLuce)(x,0)
# Voglio trovare un dì che è durato 12 ore:
k <- function(x) f(x)-12
uniroot(k, c(200, 400))$root
# 269.0581
31+28+31+30+31+30+31+31+26
# 269
# approssimativamente il 26 settembre 2011 il dì è durato 12 ore
# Voglio trovare un dì di durata minima:
g <- function(x) splinefun(giorni,OreLuce)(x,1)
uniroot(g, c(300, 400))$root
# 356.2351
31+28+31+30+31+30+31+31+30+31+30+22
# 356
# Approssimativamente il 22 dicembre 2011 è il dì più breve.
# Il dì è durato:
f(uniroot(g, c(300, 400))$root)
# 8.834599
.834599*60
# 50.07594
# Ritroviamo che è durato 8 h e 50 min
#
# Come segnare i "mesi"
g <- c(31,28,31,30,31,30,31,31,30,31,30,31,31,29,31,30,31,30,31,31,30,31,30,31)
n <- sum(g)
gg <- NULL; gg[1] <- 0; for(i in 1:24) gg[i+1] <- sum(g[1:i])
dev.new(width=6,height=4)
f <- function(x) splinefun(giorni,LuceGenova)(x,0)
plot(f,0,n,xaxt="n")
abline(v=gg, h=axTicks(2), col="blue",lty=3)
f <- function(x) splinefun(giorni,LuceMelbourne)(x,0)
plot(f,0,n,add=TRUE,col="red")
title(main="2011 e 2012, mese per mese")
abline(v=gg[13],lty=2)
# Col comando seguente otterremmo scritte piu' piccole e non in grassetto 
# title(main=list("2011 e 2012, mese per mese",cex=1,font=1))