Skip to content

Instantly share code, notes, and snippets.

@marcoarment
Created January 17, 2021 16:43
Show Gist options
  • Save marcoarment/94b173faba12f6ca3a7de9e3f1a73fa1 to your computer and use it in GitHub Desktop.
Save marcoarment/94b173faba12f6ca3a7de9e3f1a73fa1 to your computer and use it in GitHub Desktop.
Get sunrise, sunset, and the sun's position for a given time and location. Based on NOAA Solar Calculator.
// Sun.swift
// Created by Marco Arment on 1/17/21
//
// Solar-math functions are directly translated from the NOAA Solar Calculator:
// https://www.esrl.noaa.gov/gmd/grad/solcalc/
//
// This is free and unencumbered software released into the public domain.
//
// Anyone is free to copy, modify, publish, use, compile, sell, or
// distribute this software, either in source code form or as a compiled
// binary, for any purpose, commercial or non-commercial, and by any
// means.
//
// In jurisdictions that recognize copyright laws, the author or authors
// of this software dedicate any and all copyright interest in the
// software to the public domain. We make this dedication for the benefit
// of the public at large and to the detriment of our heirs and
// successors. We intend this dedication to be an overt act of
// relinquishment in perpetuity of all present and future rights to this
// software under copyright law.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
// IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//
// For more information, please refer to <http://unlicense.org/>
import Foundation
struct SunPosition {
var date: Date
var azimuth: Double = 0.0
var elevation: Double = 0.0
}
struct SunSchedule {
var position: SunPosition
var startOfDay: SunPosition
var endOfDay: SunPosition
var solarNoon: SunPosition
var solarMidnight: SunPosition
// Optional: some days in some locations never have a sunrise and/or sunset
var sunrise: SunPosition?
var sunset: SunPosition?
}
class Sun {
static func schedule(latitude: Double, longitude: Double, date: Date?) -> SunSchedule? {
return schedule(latitude: latitude, longitude: longitude, date: date, calendar: nil, timeZone: nil)
}
static func schedule(latitude: Double, longitude: Double, date: Date?, calendar: Calendar?, timeZone: TimeZone?) -> SunSchedule? {
let date = date ?? Date()
guard let JT = JulianTime(calendar: calendar ?? Calendar.current, date: date, timeZone: timeZone ?? TimeZone.current) else { return nil }
let currentPosition = azimuthAndElevation(date: date, T: JT.T, positionInDay: JT.positionInDay, latitude: latitude, longitude: longitude, timezoneOffset: JT.timezoneOffset)
let solarNoonDate = JT.startOfDay.addingTimeInterval(solNoon(JD: JT.julianDay, longitude: longitude, timezoneOffset: JT.timezoneOffset) * 60.0)
let solarMidnightDate = JT.startOfDay.addingTimeInterval(solMidnight(JD: JT.julianDay, longitude: longitude, timezoneOffset: JT.timezoneOffset) * 60.0)
guard let solarNoon = position(latitude: latitude, longitude: longitude, date: solarNoonDate, calendar: calendar, timeZone: timeZone),
let solarMidnight = position(latitude: latitude, longitude: longitude, date: solarMidnightDate, calendar: calendar, timeZone: timeZone),
let startOfDay = position(latitude: latitude, longitude: longitude, date: JT.startOfDay, calendar: calendar, timeZone: timeZone),
let endOfDay = position(latitude: latitude, longitude: longitude, date: JT.endOfDay, calendar: calendar, timeZone: timeZone)
else { return nil }
var sunrise: SunPosition? = nil
if let riseMin = sunriseOrSet(rise: true, JD: JT.julianDay, latitude: latitude, longitude: longitude, timezoneOffset: JT.timezoneOffset) {
sunrise = position(latitude: latitude, longitude: longitude, date: JT.startOfDay.addingTimeInterval(riseMin * 60.0), calendar: calendar, timeZone: timeZone)
}
var sunset: SunPosition? = nil
if let setMin = sunriseOrSet(rise: false, JD: JT.julianDay, latitude: latitude, longitude: longitude, timezoneOffset: JT.timezoneOffset) {
sunset = position(latitude: latitude, longitude: longitude, date: JT.startOfDay.addingTimeInterval(setMin * 60.0), calendar: calendar, timeZone: timeZone)
}
return SunSchedule(position: currentPosition, startOfDay: startOfDay, endOfDay: endOfDay, solarNoon: solarNoon, solarMidnight: solarMidnight, sunrise: sunrise, sunset: sunset)
}
static func position(latitude: Double, longitude: Double, date: Date?) -> SunPosition? {
return position(latitude: latitude, longitude: longitude, date: date, calendar: nil, timeZone: nil)
}
static func position(latitude: Double, longitude: Double, date: Date?, calendar: Calendar?, timeZone: TimeZone?) -> SunPosition? {
let date = date ?? Date()
guard let JT = JulianTime(calendar: calendar ?? Calendar.current, date: date, timeZone: timeZone ?? TimeZone.current) else { return nil }
return azimuthAndElevation(date: date, T: JT.T, positionInDay: JT.positionInDay, latitude: latitude, longitude: longitude, timezoneOffset: JT.timezoneOffset)
}
private struct JulianTime {
var julianDay: Double
var T: Double
var positionInDay: Double
var startOfDay: Date
var endOfDay: Date
var timezoneOffset: Double
init?(calendar: Calendar, date: Date, timeZone: TimeZone) {
var midnightBoundary = DateComponents()
midnightBoundary.hour = 0
let dateComponents = calendar.dateComponents([ .year, .month, .day ], from: date)
guard let year = dateComponents.year, let month = dateComponents.month, let day = dateComponents.day,
let dayEnd = calendar.nextDate(after: date, matching: midnightBoundary, matchingPolicy: .nextTime)
else { return nil }
self.endOfDay = dayEnd
self.startOfDay = calendar.startOfDay(for: date)
let dayLength = endOfDay.timeIntervalSince(startOfDay)
self.positionInDay = date.timeIntervalSince(startOfDay) / dayLength
self.timezoneOffset = Double(timeZone.secondsFromGMT(for: date))
self.julianDay = Sun.julianDay(year: year, month: month, day: day)
let total = julianDay + positionInDay + (timezoneOffset / dayLength)
self.T = Sun.julianCent(julianDay: total)
}
}
private static func degreesToRadians(_ degrees: Double) -> Double { return degrees * .pi / 180 }
private static func radiansToDegrees(_ radians: Double) -> Double { return radians * 180 / .pi }
// MARK: NOAA Solar Calculator
// All functions below are directly translated from https://www.esrl.noaa.gov/gmd/grad/solcalc/
// Original JavaScript: https://www.esrl.noaa.gov/gmd/grad/solcalc/main.js
private static func azimuthAndElevation(date: Date, T: Double, positionInDay: Double, latitude: Double, longitude: Double, timezoneOffset: Double) -> SunPosition {
let localtime = positionInDay * 24.0 * 60.0
let eqTime = equationOfTime(T: T) // slight diff
let theta = sunDeclination(T: T)
let solarTimeFix = eqTime + 4.0 * longitude - 60.0 * timezoneOffset / 3600.0
var trueSolarTime = localtime + solarTimeFix
while trueSolarTime > 1440.0 { trueSolarTime -= 1440.0 }
var hourAngle = trueSolarTime / 4.0 - 180.0
if hourAngle < -180.0 { hourAngle += 360.0 }
let haRad = degreesToRadians(hourAngle)
let latitudeRad = degreesToRadians(latitude)
let thetaRad = degreesToRadians(theta)
var csz = sin(latitudeRad) * sin(thetaRad) + cos(latitudeRad) * cos(thetaRad) * cos(haRad)
if csz > 1.0 { csz = 1.0 }
else if csz < -1.0 { csz = -1.0 }
let zenith = radiansToDegrees(acos(csz))
let azDenom = cos(latitudeRad) * sin(degreesToRadians(zenith))
var azimuth = 0.0
if (abs(azDenom) > 0.001) {
var azRad = (( sin(degreesToRadians(latitude)) * cos(degreesToRadians(zenith)) ) - sin(degreesToRadians(theta))) / azDenom
if (abs(azRad) > 1.0) {
if (azRad < 0) {
azRad = -1.0
} else {
azRad = 1.0
}
}
azimuth = 180.0 - radiansToDegrees(acos(azRad))
if hourAngle > 0.0 { azimuth = -azimuth }
} else {
if latitude > 0.0 { azimuth = 180.0 }
}
if azimuth < 0.0 { azimuth += 360.0 }
let exoatmElevation = 90.0 - zenith
let refractionCorrection = refraction(elevation: exoatmElevation)
let solarZen = zenith - refractionCorrection
let elevation = 90.0 - solarZen
return SunPosition(date: date, azimuth: azimuth, elevation: elevation)
}
private static func hourAngleSunrise(latitude: Double, declination: Double) -> Double {
let latRad = degreesToRadians(latitude)
let sdRad = degreesToRadians(declination)
return acos(cos(degreesToRadians(90.833))/(cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad))
}
private static func sunriseOrSetUTC(rise: Bool, JD: Double, latitude: Double, longitude: Double) -> Double {
let t = julianCent(julianDay: JD)
let eqTime = equationOfTime(T: t)
let solarDec = sunDeclination(T: t)
var hourAngle = hourAngleSunrise(latitude: latitude, declination: solarDec)
if !rise { hourAngle = -hourAngle }
let delta = longitude + radiansToDegrees(hourAngle)
return 720 - (4.0 * delta) - eqTime // in minutes
}
private static func sunriseOrSet(rise: Bool, JD: Double, latitude: Double, longitude: Double, timezoneOffset: Double) -> Double? {
let timeUTC = sunriseOrSetUTC(rise: rise, JD: JD, latitude: latitude, longitude: longitude)
let newTimeUTC = sunriseOrSetUTC(rise: rise, JD: JD + timeUTC/1440.0, latitude: latitude, longitude: longitude)
if newTimeUTC.isNaN { return nil } // no sunrise/set on this day in this location (like North/South Pole)
let timezone = timezoneOffset / 3600.0
var timeLocal = newTimeUTC + (timezone * 60.0)
if (timeLocal < 0.0) || (timeLocal >= 1440.0) {
let increment = ((timeLocal < 0.0) ? 1.0 : -1.0)
while ((timeLocal < 0.0)||(timeLocal >= 1440.0)) { timeLocal += increment * 1440.0 }
}
return timeLocal
}
private static func meanObliquityOfEcliptic(T: Double) -> Double { return 23.0 + (26.0 + ((21.448 - T * (46.8150 + T * (0.00059 - T * (0.001813))))/60.0))/60.0 }
private static func obliquityCorrection(T: Double) -> Double { return meanObliquityOfEcliptic(T: T) + 0.00256 * cos(degreesToRadians(125.04 - 1934.136 * T)) }
private static func eccentricityEarthOrbit(T: Double) -> Double { return 0.016708634 - T * (0.000042037 + 0.0000001267 * T) }
private static func geometricMeanAnomalySun(T: Double) -> Double { return 357.52911 + T * (35999.05029 - 0.0001537 * T) }
private static func sunTrueLong(T: Double) -> Double { return geometricMeanLongSun(T: T) + sunEqOfCenter(T: T) }
private static func sunApparentLong(T: Double) -> Double { return sunTrueLong(T: T) - 0.00569 - 0.00478 * sin(degreesToRadians(125.04 - 1934.136 * T)) }
private static func geometricMeanLongSun(T: Double) -> Double {
var L0 = 280.46646 + T * (36000.76983 + T * (0.0003032))
while L0 > 360.0 { L0 -= 360.0 }
while L0 < 0.0 { L0 += 360.0 }
return L0
}
private static func sunEqOfCenter(T: Double) -> Double {
let mrad = degreesToRadians(geometricMeanAnomalySun(T: T))
return sin(mrad) * (1.914602 - T * (0.004817 + 0.000014 * T)) + sin(mrad + mrad) * (0.019993 - 0.000101 * T) + sin(mrad + mrad + mrad) * 0.000289
}
private static func equationOfTime(T: Double) -> Double {
let epsilon = obliquityCorrection(T: T)
let l0 = geometricMeanLongSun(T: T)
let e = eccentricityEarthOrbit(T: T)
let m = geometricMeanAnomalySun(T: T)
var y = tan(degreesToRadians(epsilon)/2.0)
y *= y
let sin2l0 = sin(2.0 * degreesToRadians(l0))
let sinm = sin(degreesToRadians(m))
let cos2l0 = cos(2.0 * degreesToRadians(l0))
let sin4l0 = sin(4.0 * degreesToRadians(l0))
let sin2m = sin(2.0 * degreesToRadians(m))
let Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m
return radiansToDegrees(Etime) * 4.0 // in minutes of time
}
private static func sunDeclination(T: Double) -> Double {
let e = obliquityCorrection(T: T)
let lambda = sunApparentLong(T: T)
let sint = sin(degreesToRadians(e)) * sin(degreesToRadians(lambda))
return radiansToDegrees(asin(sint))
}
private static func refraction(elevation: Double) -> Double {
var correction = 0.0
if (elevation > 85.0) { return 0.0 }
let te = tan(degreesToRadians(elevation))
let te3 = te * te * te
if elevation > 5.0 {
correction = 58.1 / te - 0.07 / te3 + 0.000086 / (te3 * te * te)
} else if elevation > -0.575 {
correction = 1735.0 + elevation * (-518.2 + elevation * (103.4 + elevation * (-12.79 + elevation * 0.711) ) )
} else {
correction = -20.774 / te
}
return correction / 3600.0
}
private static func solNoon(JD: Double, longitude: Double, timezoneOffset: Double) -> Double {
let timezone = timezoneOffset / 3600.0
let tnoon = julianCent(julianDay: JD - longitude/360.0)
var eqTime = equationOfTime(T: tnoon)
let solNoonOffset = 720.0 - (longitude * 4) - eqTime
let newt = julianCent(julianDay: JD + solNoonOffset/1440.0)
eqTime = equationOfTime(T: newt)
var solNoonLocal = 720 - (longitude * 4) - eqTime + (timezone*60.0)// in minutes
while solNoonLocal < 0.0 { solNoonLocal += 1440.0 }
while solNoonLocal >= 1440.0 { solNoonLocal -= 1440.0 }
return solNoonLocal
}
private static func solMidnight(JD: Double, longitude: Double, timezoneOffset: Double) -> Double {
let timezone = timezoneOffset / 3600.0
let tnoon = julianCent(julianDay: JD - longitude/360.0)
var eqTime = equationOfTime(T: tnoon)
let solNoonOffset = 0.0 - (longitude * 4) - eqTime
let newt = julianCent(julianDay: JD + solNoonOffset/1440.0)
eqTime = equationOfTime(T: newt)
var solNoonLocal = 0.0 - (longitude * 4) - eqTime + (timezone*60.0)// in minutes
while solNoonLocal < 0.0 { solNoonLocal += 1440.0 }
while solNoonLocal >= 1440.0 { solNoonLocal -= 1440.0 }
return solNoonLocal
}
private static func julianDay(calendar: Calendar, date: Date) -> Double? {
let dateComponents = calendar.dateComponents([ .year, .month, .day ], from: date)
guard let year = dateComponents.year, let month = dateComponents.month, let day = dateComponents.day else { return nil }
return julianDay(year: year, month: month, day: day)
}
private static func julianDay(year: Int, month: Int, day: Int) -> Double {
var y = Double(year)
var m = Double(month)
let d = Double(day)
if m <= 2.0 {
y -= 1.0
m += 12.0
}
let A = floor(y / 100.0)
let B = 2 - A + floor(A / 4.0)
y += 4716.0
m += 1
return floor(365.25 * y) + floor(30.6001 * m) + (d + B - 1524.5)
}
private static func julianCent(julianDay: Double) -> Double {
return (julianDay - 2451545.0)/36525.0
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment