Compute & Plot (approx.) Solar Declination Function f(d), R

Language: Layout:

#R version 3.3.2
# Computes and plots the (approximate) Solar Declination Function f(d),
#   where d is the day of the year.
library(ggplot2)
# Define some constants first:
a=23.4368 # axial tilt of the Earth
pi=3.1416
b=(2*pi)/365.2422 # 2*pi/days per year
c=-79.5*b # for spring equinox d=~79.5
# Now define vectors for abscissa and ordinate:
d=seq(0,365,1)
y=a*sin((b*d)+c) # general formula for a modified sine wave
# Now plot the function:
title1="Solar Declination by Day of Year"
qplot(d,y,geom="path",main=title1,xlab="Day of Year",ylab="Solar Declination")+
scale_x_continuous(breaks=seq(0,366,by=61))+scale_y_continuous(breaks=seq(-25,25,5))+
# (Yes, on the graph d goes to 366, but just to make the days and months approximately 
# coincide well enough. Using 73 instead of 61 as the by= increment in seq() would even 
# the x-axis out but would also make the graph harder to read.)
# Next put in colored vertical lines at equinoxes and solstices
geom_line(aes(x=80.5,y=y),color="green")+geom_line(aes(x=263,y=y),color="orange")+
geom_line(aes(x=172,y=y),color="yellow")+geom_line(aes(x=354,y=y),color="blue")+
# Add some descriptive text to the graph
annotate("text",x=c(36,126,217,308),y=-5,label=c("Winter","Spring","Summer","Fall"))
# Done.

[ + ] Show input
Absolute running time: 10.15 sec, cpu time: 6.17 sec, memory peak: 50 Mb, absolute service time: 10,24 sec