#----------------------------------------------------------
# Simulazione analemma da Roma
# da esempio su http://www.r-bloggers.com/analemma-graphs/  
# 
# GdA 18/3/2016
#----------------------------------------------------------

# lon=12.5138  # dipartimento di Fisica
# lon=12.467   # Roma secondo http://www.maxcesar.it/calc_es.shtml

library(oce)    # dipende da gsw

png("analemma_roma_2016.png", width=700, height=700)
    
loc <- list(lon=12.5, lat=41.9) 
t <- seq.POSIXt(as.POSIXct("2015-01-01 11:10:40", tz="UTC"),
                as.POSIXct("2015-12-25 11:10:40", tz="UTC"),
                by="7 d")
sa <- sunAngle(t, lon=loc$lon, lat=loc$lat)
par(mar=c(3, 3, 3, 1), mgp=c(2, 0.7, 0)) # tighten margins
par(bg = "lightblue")
plot(sa$azimuth, sa$altitude, type='p', pch=19, col='yellow',
     xlim=c(170,190), ylim=c(0,75), 
     xlab="Azimuth (°)", ylab="Altitude (°)",
     main="Simulazione da 'Roma' (12.5E, 41.9N) alle 11:10:40 UTC, ogni 7gg")
for (i in 1:(length(t)-1)) {
  arrows(sa$azimuth[i], sa$altitude[i], sa$azimuth[i+1], sa$altitude[i+1],
         length=0.06)
}
grid()
abline(v=180, col='blue')
text(sa$azimuth[1]+0.7, sa$altitude[1] -4, "1/01/2016", pos=2, col='blue')
text(sa$azimuth[length(t)]-0.45, sa$altitude[length(t)] -4, "24/12/2016", pos=4, col='blue')
text(sa$azimuth[16]-1, sa$altitude[16]-1.5, "15/04/2016", pos=2, col='blue')

dev.off()