Sampling Plans – Analysis of Samples – Modeling – Graphing
Below are graphs of altitudes and azimuths of the sun, for differing times of the year, for Stratford, Iowa, which give you an example of the capability of R for use in planning solar installations:
You can download the R function – which plots individual altitude / azimuth graphs, for any place in the world, for any time of year. You can download the file here: angle.plot.
If you call the function ‘angle.plot’, to run the program with default values, you would type ‘angle.plot()’ at the R prompt. To change the defaults, you can manually change the defaults when the function is run; for example, angle.plot(loc=”Hull, MA”, mo=8, da=28, ye=2013, tz=4, la=42.3019, lo=70.9083).
function(loc="Stratford, Iowa", mo=1, da=1, ye=2010, tz=6,
la=42.26782, lo=93.92097) {
# The algorithms for this function come from Astrological Algorithms, 2nd ed.,
# by Jean Meeus, published in 1998, by Willmann-Bell, Inc. The page
# references are to Meeus.
# check the inputs
cat("the location entered is ", loc, "\n")
cat("the date entered is ", mo, "/", da, "/", ye, "\n")
cat("the time zone is ", tz, " hours from UT", "\n")
cat("the latitude is ", la, "degrees and the longitude is ", lo,
" degrees\n", "\n")
cat(" if not correct then enter n ")
tes <- readline()
if(tes == "n") stop("change the value of the arguments")
#
# define tosid - the solar to sidereal time correction factor (p. 408)
#
tosid <- 366.2422 / 365.2422
#
# create local - the local dynamical times for the plot - in degrees
# (see chapter 10 in Meeus)
#
local <- seq(-150, 150, by = .02)
#
# find jd - the julian days for the local dynamical times - in days (p. 61)
#
if(mo == 1 | mo == 2) {
mo <- mo + 12
ye <- ye - 1
}
ga <- floor(ye/100 + 1e-008)
gb <- 2 - ga + floor(ga/4 + 1e-008)
jd <- (floor(365.25 * (ye + 4716)) + floor(30.6001 * (mo + 1)) +
(da +.5) + gb - 1524.5) + local/(15 * 24) + lo/(15 * 24)
#
# find T - the increment from January 1, 2000 at noon TD
# for the local dynamical times - in centuries (p. 163)
#
T <- (jd - 2451545)/36525
Tsq <- T^2
#
# find L0 - the mean longitudes of the sun - in radians (p. 163)
#
L0 <- (280.46646 + 36000.76983 * T + 0.0003032 * Tsq)/180 * pi
#
# find M - the mean anomolies of the sun - in radians (p. 163)
#
M <- (357.52911 + 35999.05029 * T - 0.0001537 * Tsq)/180 * pi
#
# find C - the sun's equations of the center - in radians (p. 164)
#
C <- (1.914602 - 0.004817 * T - 1.4e-005 * Tsq) * sin(M)
C <- C + (0.019993 - 0.000101 * T) * sin(2 * M)
C <- (C + 0.000289 * sin(3 * M))/180 * pi
#
# find omega - the correction factors for finding the apparent longitudes
# of the sun - in radians (p. 164)
#
omega <- (125.04 - 1934.136 * T)/180 * pi
# find lambda - the apparent longitudes of the sun - in radians (p. 164)
#
lambda <- L0 + C - (0.005689999999999997 +
0.004779999999999999 * sin(omega))/180 * pi
#
# find epsilon - the apparent obliquities of the ecliptic - in radians
# (pp. 144-147)
#
Lprime <- (218.3165 + 481267.8813*T) / 180 * pi
deltapsi <- (-17.2*sin(omega) - 1.32*sin(2*L0) -
0.23*sin(2*Lprime) + 0.21*sin(2*omega))/3600
deltaep <- (9.2*cos(omega) + 0.57*cos(2*L0) + 0.10*cos(2*Lprime) -
0.09*cos(2*omega))/3600
epsilon <- (23.439291111 - 0.013004167 * T - 1.638889e-07 * Tsq +
5.036111111e-07 * Tsq * T + deltaep)/180 * pi
# find alpha - the apparent right ascensions of the sun - in radians (p. 165)
#
alpha <- atan2(cos(epsilon) * sin(lambda), cos(lambda))
#
# find delta - the apparent declinations of the sun - in radians (p. 165)
#
delta <- asin(sin(epsilon) * sin(lambda))
#
# find e - the eccentricities of the earth's orbit - unitless (p. 163)
#
e <- 0.016708634 - 0.000042037 * T - 0.0000001267 * Tsq
#
# find deltaT - the factors to convert to mean times from dynamical times
# and vice versa (referenced to solar time) -
# in degrees (p. 78)
#
deltaT <- (65 + 88*T + 25.3*T^2)/3600*15
#
# find jd_UT - the julian day for 0h UT on the day of the plot - in days
# and T_UT - the increment from 1/1/2000 noon UT to jd_UT - in centuries
# (p. 87)
#
jd_UT <- jd[local==0] - lo/(15 * 24) - .5
T_UT <- (jd_UT -2451545)/36525
#
# find THETA0 - the mean sidereal time of the point due south at Greenwich
# at 0h UT on the day of the plot - in radians (p. 87)
#
THETA0 <- 100.46061837 + 36000.770053608*T_UT
THETA0 <- (THETA0 + 0.000387933*T_UT^2 - (T_UT^3)/38710000)/180*pi
#
# find theta - the apparent sidereal times of the points due south at
# the location for the local dynamical times -
# in radians (p. 88)
#
theta0 <- THETA0 + (lo + local + 180 + deltaT)*tosid/180*pi
theta <- (theta0 - lo/180*pi) + (deltapsi*cos(epsilon))/180*pi
#
# find H - the hour angles of the sun for the local dynamical times -
# in radians (referenced to sidereal time) (p. 92)
#
H <- (theta - alpha)%%(2*pi)
H[H>pi] <- H[H>pi]-2*pi
#
# convert the geographical latitude to radians
#
phi
#
# find the altitudes and azimuths of the sun - in degrees (p. 93)
#
azimuth <- atan2(sin(H), (cos(H) * sin(phi) - tan(delta) *
cos(phi)))* 180 / pi
altitude <- asin(sin(phi) * sin(delta) +
cos(phi) * cos(delta) * cos(H)) * 180 / pi
#
# convert local - the local dynamical times - to Time - the civil times -
# in solar hours (p. 77)
#
Time <- (local - deltaT + 180 + lo)/15 - tz
#
# make the plot
#
par(mar = c(6, 6, 7, 7) + 0.1, font.main=1)
plot(Time[altitude>=0], altitude[altitude>=0], axes = F, type = "l",
lty = 1, xlab="", ylab = "",
xaxs = "r", xaxt = "n", yaxt = "n", yaxs = "r", ylim = c(-10, 90),
xlim = c(4, 22), col=2)
axis(1, at = c(4, 6, 8, 10, 12, 14, 16, 18, 20, 22), labels = c("4AM",
"6AM", "8AM", "10AM", "12PM", "2PM", "4PM", "6PM", "8PM",
"10PM"), tck = 1)
axis(2, at = c(-10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90), labels = c(
"-10", "0", "10", "20", "30", "40", "50", "60", "70", "80",
"90"), tck = 1, adj = 1)
legend(3.4,90, c("Altitude", "Azimuth"), lty=1,
col = 2:3, bg="white")
mtext(side = 2, line = 3, "Altitude of the Sun (degrees)")
par(new = "T")
plot(Time[altitude>=0], azimuth[altitude>=0], axes = F, type = "l",
lty = 1,
xaxt = "n", ylim = c(-140, 140), xlim = c(4, 22), ylab = "",
xlab="", col=3)
axis(side = 4, at = c(-120, -90, -60, -30, 0, 30, 60, 90, 120), labels
= c("120 E", "90 E", "60 E", "30 E", "0", "30 W", "60 W", "90 W",
"120 W"), adj = 0, cex.axis=.7)
mtext(side = 4, line = 3, "Azimuth of the Sun (degrees)", srt=270)
mtext(side = 3, line = 2, paste(loc, "\n \n", "Latitude =",
round(la, 5), " Longitude =", round(lo, dig = 5), " "))
mtext(side = 1, line = 3, "Time")
box()
#
# print out some variables
#
for( i in 1:(length(local)-1) ) {
if (altitude[i+1] >= 0 & altitude[i] < 0) {kr
if (altitude[i+1] >= altitude[i]) {kn
if (altitude[i+1] < 0 & altitude[i] >= 0) {ks
}
cat(
"L0 ",(L0[c(kr,kn,ks)]*180/pi)%%360,"\n",
"M ",(M[c(kr,kn,ks)]*180/pi)%%360,"\n",
"C ",(C[c(kr,kn,ks)]*180/pi)%%360,"\n",
"omega ",(omega[c(kr,kn,ks)]*180/pi)%%360,"\n",
"lambda ",(lambda[c(kr,kn,ks)]*180/pi)%%360,"\n",
"epsilon ",(epsilon[c(kr,kn,ks)]*180/pi)%%360,"\n",
"alpha ",(alpha[c(kr,kn,ks)]*180/pi/15)%%24,"\n",
"delta ",delta[c(kr,kn,ks)]*180/pi,"\n",
"e ",e[c(kr,kn,ks)],"\n",
"deltaT ",(deltaT[c(kr,kn,ks)]*180/pi/15)%%24,"\n",
"THETA0 ",(THETA0*180/pi/15)%%24,"\n",
"theta0 ",(theta0[c(kr,kn,ks)]*180/pi/15)%%24,"\n",
"theta ",(theta[c(kr,kn,ks)]*180/pi/15)%%24,"\n",
"H ",(H[c(kr,kn,ks)]*180/pi/15),"\n",
"azimuth ",azimuth[c(kr,kn,ks)],"\n",
"altitude ",altitude[c(kr,kn,ks)],"\n",
"Time ", Time[c(kr,kn,ks)],"\n",
"rise ",kr,"noon ",kn,"set ",ks,"\n")
}
Vanward Statistics 
