Introduction
Noon is 12 o’clock at daytime. However, the chronological time does not coincide with solar noon, when the sun reaches its highest point in the sky. This is due to three factors:
- Longitude
- The equation of time
- Daylight savings
In this section, daylight savings are not taken into consideration. But the other factors are explored to calculate the time of local solar noon. Furthermore, the timing of sun rise and sun set are estimated and graphically presented.
Longitude correction
Latitude and longitude are used to describe a location on earth. They can be expressed in degrees, minutes and seconds (DMS notation) or in digital degrees (DD notation). In software programs, it is easiest to use the DD notation. It is however straight forward to convert from one format to the other using the following functions:
R download: convert_dms_dd
R download: convert_dd_dms
In Britain, the time corresponds to the time at Greenwich (GMT) at the null meridian. When located to the east, solar time is earlier than GMT and to the west later. To calculate solar time, it is necessary to correct for this difference in longitude. The correction required is 240 seconds for every degree (1 day per 360 degrees = 24 × 60 × 60 / 360 = 240 seconds per degree).
Equation of Time
The equation of time describes the difference between mean solar time (chronological time as shown on the clock) and the apparent solar time (as shown on a sundial). This is mainly caused by two factors:
- The obliquity of the ecliptic plane, which is inclined by approximately 23 degrees relative to the earth’s equator
- The eccentricity of the earth’s orbit around the sun
The effect of both factors can be summarised in the following formula (or equation of time) with reasonable accuracy (other approximations have been described):
where,
Day number 81 is 22nd March.
For leap years, divide be 366
R download: equation_of_time
With this formula, the timing of solar noon can be estimated. Subsequently, if the declination of the sun is know, the timing of sun set and sun rise can also be estimated as described below.
Sun declination
The declination of the sun is estimated using the following formula:
Day 284 is 11th October
R download: sun_declination
Sun rise formula
Knowing the solar declination on a particular day, the sun rise and sun set times can be estimated using the following formula (all angles in radians):
where the declination can be estimated using the formula described above. To convert to radians:
To convert the latitude to radians:
If the altitude is zero, the formula gives the timing that the centre of the solar disc is at the horizon. To estimate when the top of the solar disc reaches the horizon, subtract 16 arc minutes for the semi-diameter of the solar disc and 34 arc minutes for atmospheric refraction (different under different atmospheric conditions). In total, subtract 50 arc minutes:
The parameter s corrects for the solar altitude at different timings. For sun rise and sun set, s should be set at zero. Twilights can be estimated by varying s:
- s = 0 sun rise / sun set
- s = -6 civil twilight
- s = -12 nautical twilight
- s = -18 astronomical twilight
The formula returns a value w, that can be converted to an hour angle (in radians) by taking the arccos (inverse function of the cosine) of w:
which in turn can be converted to degrees:
To calculate the sun rise and sun set times before and after solar noon (in seconds):
(1 day per 360 degrees = 24 × 60 × 60 / 360 = 240 seconds per degree).
Sun rise can be found by:
and sun set:
R download: sun_rise_formula
Graphical presentation
To allow visualisation in R, the following packages should be installed as described:
- dplyR package: 1, to manipulate data frames
- ggplot2 package: 2, to create plots
- lubridate package: 3, to work with dates
The following functions, as described above, should be copied and pasted into R and executed as described:
Location
As example the location is set at the Robert Jones and Agnes Hunt Orthopaedic Hospital NHS Foundation Trust. The coordinates are (DMS):
Longitude: 3 degrees, 1 minute and 57.8 seconds West
Latitude: 52 degrees, 53 minutes and 3.7 seconds North
To convert to DD and show the result (North and East are +, South and West are -):
library(dplyr)
library(ggplot2)
library(lubridate)
convert_dms_dd <- function(deg, min, sec){
result <- deg + min / 60 + sec / 3600
}
longitude <- convert_dms_dd(3, 1, 57.8)
# East is negative, so make negative:
longitude <- – longitude
latitude <- convert_dms_dd(52, 53, 3.7)
# show the values:
longitude
[1] -3.032722
latitude
[1] 52.88436
Create a vector of day numbers:
days <- rep(1:365, by = 1)
Load the formulas defined above:
sun_rise_formula <- function(day_number, declination, latitude, s){
# altitude correction needed; otherwise centre of the sun disk.
# – 16 arc minutes for semi-diameter of the sun disk
# horizontal refraction: – 34 arc minutes
# total – 50 arc minutes = – 50 / 60 deg
altitude_rad <- – 50 / 60 * pi / 180 + s * pi / 180
latitude_rad <- latitude * pi / 180
declination_rad <- declination * pi / 180
w <- (sin(altitude_rad) – (sin(latitude_rad) * sin(declination_rad))) /
(cos(latitude_rad) * cos(declination_rad))
hour_angle_rad <- acos(w)
hour_angle <- hour_angle_rad * 180 / pi
rise_set_seconds <- hour_angle * 240
result <- rise_set_seconds
}
sun_declination <- function(day_number){
result <- 23.45 * sin(2 * pi / 365 * (day_number + 284))
}
equation_of_time <- function(day_number){
# day_number is the number of the day in the year
# 81 is the 22nd March
# B in radians
B <- 2 * pi * (day_number – 81) / 365
# correction in minutes
result <- 9.87 * sin(2 * B) – 7.53 * cos(B) – 1.5 * sin(B)
# correction in seconds
result <- result * 60
}
Calculate the correction (equation of time), the solar declination and the time of sun rise / sun set before / after solar noon:
correction <- equation_of_time(day_number = days)
declination <- sun_declination(day_number = days)
rise_set_seconds <- sun_rise_formula(day_number = days, declination = declination, latitude = latitude, s = 0)
Convert to a data frame (called solar_time):
solar_time <- data.frame(“day” = days, “correction” = correction, “declination” = declination, “rise_set” = rise_set_seconds)
Calculate the solar altitude from the solar declination and add to the data frame:
solar_time$altitude <- 90 – latitude + solar_time$declination
Add the longitude correction and the total time correction (longitude and equation of time) to the data frame:
solar_time$longitude_correction <- longitude * 240
solar_time$total_correction <- solar_time$longitude_correction + solar_time$correction
Also add the start time of the day and the time of chronological noon to the data frame:
solar_time$start <- ymd_hms(“2017-12-31 00:00:00″, tz = “GMT”) + ddays(solar_time$day)
solar_time$noon <- ymd_hms(“2017-12-31 12:00:00″, tz = “GMT”) + ddays(solar_time$day)
Calculate apparent solar noon (as seen on a sun dial) and add to the data frame:
solar_time$noon_apparent <- solar_time$noon – dseconds(solar_time$total_correction)
Calculate sun set, sun rise and the total day length and add to the data frame:
solar_time$sun_rise <- solar_time$noon_apparent – dseconds(solar_time$rise_set)
solar_time$sun_set <- solar_time$noon_apparent + dseconds(solar_time$rise_set)
solar_time$day_length <- solar_time$sun_set – solar_time$sun_rise
Similarly, calculate the timings of civil, nautical and astronomical twilight and add to the data frame:
# add civil twilight s = – 6
# add nautical twilight s = – 12
# add astronomical twilight, s = – 18
solar_time$civil_seconds <- sun_rise_formula(day_number = days, declination = declination, latitude = latitude, s = – 6)
solar_time$nautical_seconds <- sun_rise_formula(day_number = days, declination = declination, latitude = latitude, s = – 12)
solar_time$astro_seconds <- sun_rise_formula(day_number = days, declination = declination, latitude = latitude, s = – 18)
Warning message:
In acos(w) : NaNs produced
solar_time$civil_rise <- solar_time$noon_apparent – dseconds(solar_time$civil_seconds)
solar_time$civil_set <- solar_time$noon_apparent + dseconds(solar_time$civil_seconds)
solar_time$nautical_rise <- solar_time$noon_apparent – dseconds(solar_time$nautical_seconds)
solar_time$nautical_set <- solar_time$noon_apparent + dseconds(solar_time$nautical_seconds)
solar_time$astro_rise <- solar_time$noon_apparent – dseconds(solar_time$astro_seconds)
solar_time$astro_set <- solar_time$noon_apparent + dseconds(solar_time$astro_seconds)
The warning message shows that NaNs are introduced when there is no astronomical twilight in the height of summer.
The data frame is now complete and can be reviewed with:
str(solar_time)
head(solar_time)
tail(solar_time)
solar_time
Output omitted for brevity.
Retrieve the following from the data frame and define them for later use:
- earliest sun rise
- latest sun rise
- earliest sun set
- latest sun set
- shortest day
- longest day
earliest_sun_rise <- solar_time %>%
dplyr::mutate(sun_rise_hour = hour(sun_rise) + minute(sun_rise) / 60 + second(sun_rise) / 3600) %>%
dplyr::filter(sun_rise_hour == min(sun_rise_hour))
latest_sun_rise <- solar_time %>%
dplyr::mutate(sun_rise_hour = hour(sun_rise) + minute(sun_rise) / 60 + second(sun_rise) / 3600) %>%
dplyr::filter(sun_rise_hour == max(sun_rise_hour))
earliest_sun_set <- solar_time %>%
dplyr::mutate(sun_set_hour = hour(sun_set) + minute(sun_set) / 60 + second(sun_set) / 3600) %>%
dplyr::filter(sun_set_hour == min(sun_set_hour))
latest_sun_set <- solar_time %>%
dplyr::mutate(sun_set_hour = hour(sun_set) + minute(sun_set) / 60 + second(sun_set) / 3600) %>%
dplyr::filter(sun_set_hour == max(sun_set_hour))
shortest_day <- solar_time %>%
dplyr::filter(day_length == min(day_length))
longest_day <- solar_time %>%
dplyr::filter(day_length == max(day_length))
earliest_sun_rise
day 168, 17/06
latest_sun_rise
day 364, 30/12
earliest_sun_set
day 346, 12/12
latest_sun_set
day 176, 25/6
shortest_day
day 355, 21/12
longest_day
day 172, 21/06
Output omitted, but summarised for readability.
Although the longest day is 21st June, the earliest sun rise is a few days earlier (17th June) and the latest sun set a few days later ((25th June). Similarly, around the shortest day (21st December) the latest sun rise is a few days later (30th of December) and the earliest sun set a few days earlier (12th December). This is better shown graphically in a plot that follows later.
Equation of time
It is now straight forward to plot the equation of time with ggplot2 2:
dev.new()
solar_time %>%
ggplot(aes(x = day, y = correction)) +
geom_line(colour = “blue”) +
scale_x_continuous(“Day Number”) +
scale_y_continuous(“Correction [seconds]”) +
ggtitle(“Equation of Time”) +
theme_bw()
Or add the current day to the plot:
dev.new()
solar_time %>%
ggplot(aes(x = day, y = correction)) +
geom_line(colour = “blue”) +
geom_point(data = subset(solar_time, solar_time$day == yday(now(tzone = “GMT”))), aes(x = day, y = correction), colour = “black”) +
geom_text(data = subset(solar_time, solar_time$day == yday(now(tzone = “GMT”))), aes(x = day, y = correction), label = today(), hjust = 1.1) +
scale_x_continuous(“Day Number”) +
scale_y_continuous(“Correction [seconds]”) +
ggtitle(“Equation of Time”) +
theme_bw()
To show a plot of the sun rise in red, sun set in orange and day length in green:
dev.new()
solar_time %>%
ggplot() +
geom_line(aes(x = day, y = (sun_rise – start)), colour = “red”) +
geom_line(aes(x = day, y = (sun_set – start)), colour = “orange”) +
geom_line(aes(x = day, y = day_length), colour = ‘darkgreen’) +
geom_vline(xintercept = longest_day$day, colour = “darkgreen”) +
geom_vline(xintercept = shortest_day$day, colour = “darkgreen”) +
geom_vline(xintercept = earliest_sun_rise$day, colour = “red”, linetype = 2) +
geom_vline(xintercept = latest_sun_rise$day, colour = “red”, linetype = 2) +
geom_vline(xintercept = earliest_sun_set$day, colour = “orange”, linetype = 2) +
geom_vline(xintercept = latest_sun_set$day, colour = “orange”, linetype = 2) +
ggtitle(“Sun Rise (red), Sun Set (orange) and Day Length (green)”) +
scale_x_continuous(“Day Number”, breaks = c(1, 100, latest_sun_set$day, longest_day$day, earliest_sun_rise$day, 200, 300, latest_sun_rise$day, shortest_day$day, earliest_sun_set$day), limits = c(0, 385)) +
geom_text(data = latest_sun_set, aes(x = day, y = 10, label = date(start)), colour = “orange”, hjust = – 0.1) +
geom_text(data = earliest_sun_rise, aes(x = day, y = 10, label = date(start)), colour = “red”, hjust = 1.1) +
geom_text(data = earliest_sun_set, aes(x = day, y = 10, label = date(start)), colour = “orange”, hjust = 1.1) +
geom_text(data = latest_sun_rise, aes(x = day, y = 10, label = date(start)), colour = “red”, hjust = – 0.1) +
geom_text(data = longest_day, aes(x = day, y = 17.5 , label = date(start)), colour = “darkgreen”, hjust = – 0.1) +
geom_text(data = shortest_day, aes(x = day, y = 17.5, label = date(start)), colour = “darkgreen”, hjust = – 0.1) +
scale_y_continuous(“Hours after Midnight”) +
theme_bw() +
theme(axis.text.x=element_text(angle = 90, hjust = 0, vjust = 0)) # after theme_bw call otherwise it doesnt’ work!
An analemna, shows the position of the sun in the sky every day at midday (chronological noon):
dev.new()
solar_time %>%
ggplot(aes(x = correction, y = declination)) +
geom_point(colour = “orange”) +
scale_x_continuous(“Correction [seconds]”) +
scale_y_continuous(“Declination [degrees]”) +
ggtitle(“Analemna – Solar Declination”) +
theme_bw()
Annotate different days on the plot:
dev.new()
solar_time %>%
ggplot(aes(x = correction, y = declination)) +
geom_point(colour = “grey”) +
scale_x_continuous(“Correction [seconds]”) +
scale_y_continuous(“Declination [degrees]”, limits = c(-25, 25)) +
ggtitle(“Analemna – Solar Declination”) +
# geom_text(data = subset(solar_time, solar_time$day == 1), aes(x = correction, y = declination, label = day), vjust = 2) +
# geom_point(data = subset(solar_time, solar_time$day == 1), aes(x = correction, y = declination), colour = “black”) +
geom_text(data = longest_day, aes(x = correction, y = declination, label = day), colour = “darkgreen”, vjust = – 1) +
geom_point(data = longest_day, aes(x = correction, y = declination), colour = “darkgreen”) +
geom_text(data = shortest_day, aes(x = correction, y = declination, label = day), colour = “darkgreen”, vjust = 2) +
geom_point(data = shortest_day, aes(x = correction, y = declination), colour = “darkgreen”) +
geom_text(data = latest_sun_set, aes(x = correction, y = declination, label = day), colour = “orange”, hjust = 0.7, vjust = 2) +
geom_point(data = latest_sun_set, aes(x = correction, y = declination), colour = “orange”) +
geom_text(data = latest_sun_rise, aes(x = correction, y = declination, label = day), colour = “red”, hjust = 0.5, vjust = 2) +
geom_point(data = latest_sun_rise, aes(x = correction, y = declination), colour = “red”) +
geom_text(data = earliest_sun_set, aes(x = correction, y = declination, label = day), colour = “orange”, vjust = 2) +
geom_point(data = earliest_sun_set, aes(x = correction, y = declination), colour = “orange”) +
geom_text(data = earliest_sun_rise, aes(x = correction, y = declination, label = day), colour = “red”, hjust = 0.1 , vjust = 2) +
geom_point(data = earliest_sun_rise, aes(x = correction, y = declination), colour = “red”) +
geom_point(data = subset(solar_time, solar_time$day == 104), aes(x = correction, y = declination), colour = “black”) +
geom_text(data = subset(solar_time, solar_time$day == 104), aes(x = correction, y = declination, label = day), hjust = 0, vjust = 2) +
geom_point(data = subset(solar_time, solar_time$day == 241), aes(x = correction, y = declination), colour = “black”) +
geom_text(data = subset(solar_time, solar_time$day == 241), aes(x = correction, y = declination, label = day), hjust = 0.2, vjust = -1) +
# add equinox in blue, table shows equinox on 20th March: day 79 and 23rd September day 266
geom_text(data = subset(solar_time, solar_time$day == 79), aes(x = correction, y = declination, label = day), colour = “blue”, hjust = – 0.3, vjust = 1) +
geom_point(data = subset(solar_time, solar_time$day == 79), aes(x = correction, y = declination), colour = “blue”) +
geom_text(data = subset(solar_time, solar_time$day == 266), aes(x = correction, y = declination, label = day), colour = “blue”, hjust = – 0.2, vjust = 0) +
geom_point(data = subset(solar_time, solar_time$day == 266), aes(x = correction, y = declination), colour = “blue”) +
# add the extreme corrections:
geom_point(data = subset(solar_time, solar_time$day == 44), aes(x = correction, y = declination), colour = “salmon”) +
geom_text(data = subset(solar_time, solar_time$day == 44), aes(x = correction, y = declination, label = day), colour = “salmon”, hjust = – 0.5, vjust = 0) +
geom_point(data = subset(solar_time, solar_time$day == 303), aes(x = correction, y = declination), colour = “salmon”) +
geom_text(data = subset(solar_time, solar_time$day == 303), aes(x = correction, y = declination, label = day), colour = “salmon”, hjust = 1.4, vjust = 0) +
theme_bw()
The plot is easier to read with text:
dev.new()
solar_time %>%
ggplot(aes(x = correction, y = declination)) +
geom_point(colour = “grey”) +
scale_x_continuous(“Correction [seconds]”) +
scale_y_continuous(“Declination [degrees]”, limits = c(- 30, 30)) +
ggtitle(“Analemna – Solar Declination”) +
# geom_text(data = subset(solar_time, solar_time$day == 1), aes(x = correction, y = declination, label = day), vjust = 2) +
# geom_point(data = subset(solar_time, solar_time$day == 1), aes(x = correction, y = declination), colour = “black”) +
geom_text(data = longest_day, aes(x = correction, y = declination), colour = “darkgreen”, label = “longest\nday”, vjust = – 0.5) +
geom_point(data = longest_day, aes(x = correction, y = declination), colour = “darkgreen”) +
geom_text(data = shortest_day, aes(x = correction, y = declination), colour = “darkgreen”, label = “shortest\nday”, vjust = 1.5) +
geom_point(data = shortest_day, aes(x = correction, y = declination), colour = “darkgreen”) +
geom_text(data = latest_sun_set, aes(x = correction, y = declination), colour = “orange”, label = “latest\nsun set”, hjust = 0.7, vjust = 1.5) +
geom_point(data = latest_sun_set, aes(x = correction, y = declination), colour = “orange”) +
geom_text(data = latest_sun_rise, aes(x = correction, y = declination), colour = “red”, label = “latest\nsun rise”, hjust = 0.5, vjust = 1.5) +
geom_point(data = latest_sun_rise, aes(x = correction, y = declination), colour = “red”) +
geom_text(data = earliest_sun_set, aes(x = correction, y = declination, label = day), colour = “orange”, label = “earliest\nsun set”, vjust = 1.5) +
geom_point(data = earliest_sun_set, aes(x = correction, y = declination), colour = “orange”) +
geom_text(data = earliest_sun_rise, aes(x = correction, y = declination), colour = “red”, label = “earliest\nsun rise”, hjust = 0.1 , vjust = 1.5) +
geom_point(data = earliest_sun_rise, aes(x = correction, y = declination), colour = “red”) +
geom_text(data = subset(solar_time, solar_time$day == 79), aes(x = correction, y = declination), colour = “blue”, label = “vernal equinox”, hjust = – 0.1, vjust = 1) +
geom_point(data = subset(solar_time, solar_time$day == 79), aes(x = correction, y = declination), colour = “blue”) +
geom_text(data = subset(solar_time, solar_time$day == 266), aes(x = correction, y = declination), colour = “blue”, label = “autumnal equinox”, hjust = – 0.1, vjust = 0) +
geom_point(data = subset(solar_time, solar_time$day == 266), aes(x = correction, y = declination), colour = “blue”) +
theme_bw()
Rather than the solar declination, plot the solar altitude of the location (RJAH Hospital) on the y-axis and indicate the current day with the timing of solar noon:
dev.new()
solar_time %>%
ggplot(aes(x = total_correction / 60, y = altitude)) + # correction in mintues
geom_point(colour = “red”, alpha = 0.3) +
scale_x_continuous(“Correction [minutes]”) +
scale_y_continuous(“Altitude [degrees]”) +
ggtitle(“Alanemna\nSolar Altitude at RJAH Hospital”) +
geom_text(data = subset(solar_time, solar_time$day == yday(now(tzone = “GMT”))), aes(x = total_correction / 60, y = altitude, label = noon_apparent), hjust = 1, vjust = – 0.5) +
geom_point(data = subset(solar_time, solar_time$day == yday(now(tzone = “GMT”))), aes(x = total_correction / 60, y = altitude), colour = “black”) +
theme_bw()
To show the timings of the different twilights (civil, nautical and astronomical):
dev.new()
solar_time %>%
ggplot() +
geom_line(aes(x = day, y = (sun_rise – start)), colour = “red”) +
geom_line(aes(x = day, y = (sun_set – start)), colour = “orange”) +
# geom_line(aes(x = day, y = day_length), colour = ‘darkgreen’) +
geom_line(aes(x = day, y = (astro_rise – start) / 60), colour = “darkblue”) + # gives difference in minutes
geom_line(aes(x = day, y = (astro_set – start)), colour = “darkblue”) +
geom_line(aes(x = day, y = (nautical_rise – start)), colour = “blue”) +
geom_line(aes(x = day, y = (nautical_set – start)), colour = “blue”) +
geom_line(aes(x = day, y = (civil_rise – start)), colour = “steelblue”) +
geom_line(aes(x = day, y = (civil_set – start)), colour = “steelblue”) +
# add longest / shortest day earliest / latest sunrise / sunset
geom_vline(xintercept = longest_day$day, colour = “darkgreen”) +
geom_vline(xintercept = shortest_day$day, colour = “darkgreen”) +
geom_vline(xintercept = earliest_sun_rise$day, colour = “red”, linetype = 2) +
geom_vline(xintercept = latest_sun_rise$day, colour = “red”, linetype = 2) +
geom_vline(xintercept = earliest_sun_set$day, colour = “orange”, linetype = 2) +
geom_vline(xintercept = latest_sun_set$day, colour = “orange”, linetype = 2) +
# equinoxes
geom_vline(xintercept = 79, colour = “blue”, linetype = 3) +
geom_vline(xintercept = 266, colour = “blue”, linetype = 3) +
# astro twilight limits
geom_vline(xintercept = 132, colour = “darkblue”, linetype = 2) +
geom_vline(xintercept = 212, colour = “darkblue”, linetype = 2) +
scale_x_continuous(“Day Number”, limits = c(0, 375), breaks = c()) +
scale_x_continuous(“Day Number”, breaks = c(1, 79, 132, 212, 266, latest_sun_set$day, longest_day$day, earliest_sun_rise$day, latest_sun_rise$day, shortest_day$day, earliest_sun_set$day), limits = c(0, 375)) +
scale_y_continuous(“Hours after Midnight”) +
ggtitle(“Civil, Nautical and Astronomical Twilight”) +
theme_bw() +
theme(axis.text.x=element_text(angle = 90, hjust = 0, vjust = 0)) # after theme_bw call otherwise it doesnt’ work!
Scale for ‘x’ is already present. Adding another scale for ‘x’, which will replace the existing scale.
It is, of course, also possible to calculate the solar time from the computer time and estimate the timing of apparent solar noon:
###### calculate solar time based on computer time.
# use GMT, ignore British Summer Time
computer_time <- now(tzone = “GMT”)
year_day <- yday(computer_time)
correction_today <- equation_of_time(year_day) + longitude * 240
solar_time_today <- computer_time – correction_today
mean_solar_noon_today <- today() + dhours(12)
tz(mean_solar_noon_today) <- “GMT”
apparent_solar_noon_today <- mean_solar_noon_today – correction_today
computer_time
[1] “2018-04-13 15:37:41 GMT”
year_day
[1] 103
correction_today
[1] -774.042
solar_time_today
[1] “2018-04-13 15:50:35 GMT”
mean_solar_noon_today
[1] “2018-04-13 12:00:00 GMT”
apparent_solar_noon_today
[1] “2018-04-13 12:12:54 GMT”