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
}