Skip to content

Instantly share code, notes, and snippets.

@neilt

neilt/gist:6bc284ac7ce00d566002bc45bc0d86dd

Created June 16, 2019 16:08
Show Gist options
  • Save neilt/6bc284ac7ce00d566002bc45bc0d86dd to your computer and use it in GitHub Desktop.
Save neilt/6bc284ac7ce00d566002bc45bc0d86dd to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.txt
//
//
// NTSolar.swift
//
// Created by Neil Tiffin on 5/8/19.
// Copyright © 2019 Performance Champions, Inc.
// Copyright © 2019 Neil Tiffin.
//
// Released to the public domain by Neil Tiffin, May 2019
// Released to the public domain by Performance Champions, Inc., May 2019
//
import Foundation
import CoreLocation
/// Class to calculate sunrise sunset.
///
/// C code originally from: [http://stjarnhimlen.se/comp/sunriset.c](http://stjarnhimlen.se/comp/sunriset.c)
class NTSolar {
// MARK: - Public Swift Interface
/// Calculate the sun rise and set times.
///
/// - Parameters:
/// - forDate: The date for the calculation. You should ensure that the date, which is stored as UTC, is the date
/// really wanted in the given time zone. It will be converted to the given time zone before being used.
/// - atLocation: The latitude and longitude for the calculation.
/// - inTimeZone: The time zone for the resulting date and times.
/// - Returns: If the sun both rises and sets on the day requested and at the
/// latitude requested then return sun rise and set times rounded down to the minute, nil othewise.
class func sunRiseAndSet(forDate: Date,
atLocation: CLLocationCoordinate2D,
inTimeZone: TimeZone) -> (sunrise: Date, sunset: Date)? {
var calendar = Calendar(identifier: .gregorian)
calendar.timeZone = inTimeZone
var comp = calendar.dateComponents([.day, .year, .month], from: forDate)
comp.calendar = calendar
guard let year = comp.year,
let month = comp.month,
let day = comp.day else {
gLog.error("Failed to find date components.")
return nil
}
let (riseUTC, setUTC, code) = NTSolar.sun_rise_set(year: year,
month: month,
day: day,
lon: atLocation.longitude,
lat: atLocation.latitude)
if code != .RiseAndSet {
gLog.error("Failed to find rise and set.")
return nil
}
// Calc sunrise
var riseLocalHrs = riseUTC + (Double(inTimeZone.secondsFromGMT()) / 3600.0)
if riseLocalHrs > 24.0 {
riseLocalHrs -= 24.0
}
let riseHoursInt = Int(riseLocalHrs)
comp.hour = riseHoursInt
comp.minute = Int( (( riseLocalHrs - Double(riseHoursInt) ) * 60.0).rounded() )
comp.second = 0
comp.nanosecond = 0
guard let riseTime = comp.date else {
gLog.error("\nFailed to calculate rise time hrs: \(riseLocalHrs).\n\(comp)\n")
return nil
}
// Calc sunset
var setLocalHrs = setUTC + (Double(inTimeZone.secondsFromGMT()) / 3600.0)
if setLocalHrs > 24.0 {
setLocalHrs -= 24.0
}
let setHoursInt = Int(setLocalHrs)
comp.hour = setHoursInt
comp.minute = Int( ((setLocalHrs - Double(setHoursInt) ) * 60.0).rounded() )
comp.second = 0
comp.nanosecond = 0
guard let setTime = comp.date else {
gLog.error("\nFailed to calculate set time hrs: \(setLocalHrs).\n\(comp)\n")
return nil
}
return (riseTime, setTime)
}
// MARK: - SUNRISET.C
// C code originally from: [http://stjarnhimlen.se/comp/sunriset.c](http://stjarnhimlen.se/comp/sunriset.c)
//
// The conversion process removed pointers in favor of return tuples, converted macros to function calls,
// added return code enum, and converted
// comments to work with Xcode. Of course the C code had to be converted to Swift, but that was minimal.
//
// As much as possible the original code was left intact in order to not introduce bugs.
// In other words, some code was left unfashionable by today's standards.
enum ReturnCode: Int {
/// Sun is below the specified "horizon" 24 hours
/// "Day" length = 0 hours, trise and tset are
/// both set to the time when the sun is at south.
case SunAlwaysBelow = -1
/// Sun rises/sets this day, times stored at rise and set.
case RiseAndSet = 0
/// Sun above the specified "horizon" 24 hours.
/// trise set to time when the sun is at south,
/// minus 12 hours while tset is set to the south
/// time plus 12 hours. "Day" length = 24 hours
case SunAlwaysAbove = 1
}
/* +++Date last modified: 05-Jul-1997 */
/* Updated comments, 05-Aug-2013 */
/*
SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
the length of the day at any date and latitude
Written as DAYLEN.C, 1989-08-16
Modified to SUNRISET.C, 1992-12-01
(c) Paul Schlyter, 1989, 1992
Released to the public domain by Paul Schlyter, December 1992
*/
/// A macro to compute the number of days elapsed since 2000 Jan 0.0
/// (which is equal to 1999 Dec 31, 0h UT)
private class func days_since_2000_Jan_0(y:Int, m:Int, d:Int) -> Int {
return (367*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530)
}
/* Some conversion factors between radians and degrees */
private static let PI = 3.1415926535897932384
private static let RADEG = ( 180.0 / PI )
private static let DEGRAD = ( PI / 180.0 )
/* The trigonometric functions in degrees */
private class func sind(x: Double) -> Double { return sin((x)*DEGRAD) }
private class func cosd(x: Double) -> Double { return cos((x)*DEGRAD) }
private class func tand(x: Double) -> Double { return tan((x)*DEGRAD) }
private class func atand(x: Double) -> Double { return (RADEG*atan(x)) }
private class func asind(x: Double) -> Double { return (RADEG*asin(x)) }
private class func acosd(x: Double) -> Double { return (RADEG*acos(x)) }
private class func atan2d(y: Double,x: Double)-> Double { return (RADEG*atan2(y,x)) }
/* Following are some macros around the "workhorse" function __daylen__ */
/* They mainly fill in the desired values for the reference altitude */
/* below the horizon, and also selects whether this altitude should */
/* refer to the Sun's center or its upper limb. */
/** This macro computes the length of the day, from sunrise to sunset. */
/** Sunrise/set is considered to occur when the Sun's upper limb is */
/** 35 arc minutes below the horizon (this accounts for the refraction */
/** of the Earth's atmosphere). */
class func day_length(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> Double {
return daylen( year: year, month: month, day: day, lon: lon, lat: lat, altit: -35.0/60.0, upper_limb: 1 )
}
/** This macro computes the length of the day, including civil twilight. */
/** Civil twilight starts/ends when the Sun's center is 6 degrees below */
/** the horizon. */
class func day_civil_twilight_length(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> Double {
return daylen( year: year, month: month, day: day, lon: lon, lat: lat, altit: -6.0, upper_limb: 0 )
}
/** This macro computes the length of the day, incl. nautical twilight. */
/** Nautical twilight starts/ends when the Sun's center is 12 degrees */
/** below the horizon. */
class func day_nautical_twilight_length(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> Double {
return daylen( year: year, month: month, day: day, lon: lon, lat: lat, altit: -12.0, upper_limb: 0 )
}
/** This macro computes the length of the day, incl. astronomical twilight. */
/** Astronomical twilight starts/ends when the Sun's center is 18 degrees */
/** below the horizon. */
class func day_astronomical_twilight_length(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> Double {
return daylen( year: year, month: month, day: day, lon: lon, lat: lat, altit: -18.0, upper_limb: 0 )
}
/** This macro computes times for sunrise/sunset. */
/** Sunrise/set is considered to occur when the Sun's upper limb is */
/** 35 arc minutes below the horizon (this accounts for the refraction */
/** of the Earth's atmosphere). */
class func sun_rise_set(year: Int, month: Int, day: Int , lon: Double, lat: Double) -> (trise: Double, tset: Double, code: ReturnCode) {
let (start, end, code) = sunriset( year: year, month: month, day: day, lon: lon, lat: lat, altit: -35.0/60.0, upper_limb: 1)
return (start, end, code)
}
/** This macro computes the start and end times of civil twilight. */
/** Civil twilight starts/ends when the Sun's center is 6 degrees below */
/** the horizon. */
class func civil_twilight(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> (trise: Double, tset: Double, code: ReturnCode) {
let (start, end, code) = sunriset( year: year, month: month, day: day, lon: lon, lat: lat, altit: -6.0, upper_limb: 0)
return (start, end, code)
}
/** This macro computes the start and end times of nautical twilight. */
/** Nautical twilight starts/ends when the Sun's center is 12 degrees */
/** below the horizon. */
class func nautical_twilight(year: Int, month: Int ,day: Int, lon: Double, lat: Double) -> (trise: Double, tset: Double, code: ReturnCode) {
let (start, end, code) = sunriset( year: year, month: month, day: day, lon: lon, lat: lat, altit: -12.0, upper_limb: 0)
return (start, end, code)
}
/** This macro computes the start and end times of astronomical twilight. */
/** Astronomical twilight starts/ends when the Sun's center is 18 degrees */
/** below the horizon. */
class func astronomical_twilight(year: Int, month: Int, day: Int, lon: Double, lat: Double) -> (trise: Double, tset: Double, code: ReturnCode) {
let (start, end, code) = sunriset( year: year, month: month, day: day, lon: lon, lat: lat, altit: -18.0, upper_limb: 0)
return (start, end, code)
}
/// The "workhorse" function for sun rise/set times
///
/// - Parameters:
/// - year: calendar date, 1801-2099 only.
/// - month: calendar date, 1801-2099 only.
/// - day: calendar date, 1801-2099 only.
/// - lon: Eastern longitude positive, Western longitude negative. The longitude value IS critical in this function!
/// - lat: Northern latitude positive, Southern latitude negative
/// - altit: the altitude which the Sun should cross. Set to -35/60 degrees for rise/set, -6 degrees
/// for civil, -12 degrees for nautical and -18 degrees for astronomical twilight.
/// - upper_limb: non-zero -> upper limb, zero -> center
/// Set to non-zero (e.g. 1) when computing rise/set
/// times, and to zero when computing start/end of twilight.
/// - Returns: rise, set, code.
///
/// Both times in hours UT are relative to the specified altitude,
/// and thus this function can be used to compute
/// various twilight times, as well as rise/set times.
///
/// Code 0 = sun rises/sets this day, times stored at rise and set.
///
/// Code +1 = sun above the specified "horizon" 24 hours.
/// *trise set to time when the sun is at south,
/// minus 12 hours while *tset is set to the south
/// time plus 12 hours. "Day" length = 24 hours
///
/// Code -1 = sun is below the specified "horizon" 24 hours
/// "Day" length = 0 hours, *trise and *tset are
/// both set to the time when the sun is at south.
private class func sunriset(year: Int,
month: Int,
day: Int,
lon: Double,
lat: Double,
altit: Double,
upper_limb: Int) -> (trise: Double, tset: Double, code: ReturnCode) {
var altit = altit
var d: Double /* Days since 2000 Jan 0.0 (negative before) */
var sr: Double /* Solar distance, astronomical units */
var sRA: Double /* Sun's Right Ascension */
var sdec: Double /* Sun's declination */
var sradius: Double /* Sun's apparent radius */
var t: Double /* Diurnal arc */
var tsouth: Double /* Time when Sun is at south */
var sidtime: Double /* Local sidereal time */
var rc: ReturnCode = ReturnCode.RiseAndSet /* Return cde from function - usually 0 */
/* Compute d of 12h local mean solar time */
d = Double(days_since_2000_Jan_0(y: year,m: month,d: day)) + 0.5 - lon/360.0;
/* Compute the local sidereal time of this moment */
sidtime = revolution( x: GMST0(d: d) + 180.0 + lon )
/* Compute Sun's RA, Decl and distance at this moment */
(sRA, sdec, sr) = sun_RA_dec(d: d)
/* Compute time when Sun is at south - in hours UT */
tsouth = 12.0 - rev180(x: sidtime - sRA)/15.0
/* Compute the Sun's apparent radius in degrees */
sradius = 0.2666 / sr
/* Do correction to upper limb, if necessary */
if upper_limb != 0 {
altit -= sradius
}
/* Compute the diurnal arc that the Sun traverses to reach */
/* the specified altitude altit: */
do {
let cost: Double = ( sind(x: altit) - sind(x: lat) * sind(x: sdec) ) / ( cosd(x: lat) * cosd(x: sdec) );
if ( cost >= 1.0 ) {
rc = ReturnCode.SunAlwaysBelow
t = 0.0 /* Sun always below altit */
}
else if ( cost <= -1.0 ) {
rc = ReturnCode.SunAlwaysAbove
t = 12.0 /* Sun always above altit */
}
else {
t = acosd(x: cost)/15.0 /* The diurnal arc, hours */
}
}
/* Store rise and set times - in hours UT */
let trise = tsouth - t;
let tset = tsouth + t;
return (trise, tset, rc)
} /* __sunriset__ */
/// The "workhorse" function
///
/// - Parameters:
/// - year: year,month,date = calendar date, 1801-2099 only.
/// - month: year,month,date = calendar date, 1801-2099 only.
/// - day: year,month,date = calendar date, 1801-2099 only.
/// - lon: Eastern longitude positive, Western longitude negative
/// - lat: Northern latitude positive, Southern latitude negative
/// - altit: altit = the altitude which the Sun should cross
/// Set to -35/60 degrees for rise/set, -6 degrees
/// for civil, -12 degrees for nautical and -18
/// degrees for astronomical twilight.
/// - upper_limb: upper_limb: non-zero -> upper limb, zero -> center
/// Set to non-zero (e.g. 1) when computing day length
/// and to zero when computing day+twilight length.
/// - Returns: Day number
private class func daylen(year: Int, month: Int, day: Int, lon: Double, lat: Double,
altit: Double, upper_limb: Int ) -> Double {
/**********************************************************************/
/* Note: year,month,date = calendar date, 1801-2099 only. */
/* Eastern longitude positive, Western longitude negative */
/* Northern latitude positive, Southern latitude negative */
/* The longitude value is not critical. Set it to the correct */
/* longitude if you're picky, otherwise set to to, say, 0.0 */
/* The latitude however IS critical - be sure to get it correct */
/* altit = the altitude which the Sun should cross */
/* Set to -35/60 degrees for rise/set, -6 degrees */
/* for civil, -12 degrees for nautical and -18 */
/* degrees for astronomical twilight. */
/* upper_limb: non-zero -> upper limb, zero -> center */
/* Set to non-zero (e.g. 1) when computing day length */
/* and to zero when computing day+twilight length. */
/**********************************************************************/
var altit = altit
var d: Double /* Days since 2000 Jan 0.0 (negative before) */
var obl_ecl: Double /* Obliquity (inclination) of Earth's axis */
var sr: Double /* Solar distance, astronomical units */
var slon: Double /* True solar longitude */
var sin_sdecl: Double /* Sine of Sun's declination */
var cos_sdecl: Double /* Cosine of Sun's declination */
var sradius: Double /* Sun's apparent radius */
var t: Double /* Diurnal arc */
/* Compute d of 12h local mean solar time */
d = Double(days_since_2000_Jan_0(y: year, m: month, d: day)) + 0.5 - lon/360.0;
/* Compute obliquity of ecliptic (inclination of Earth's axis) */
obl_ecl = 23.4393 - 3.563E-7 * d;
/* Compute Sun's ecliptic longitude and distance */
(slon, sr) = sunpos( d: d )
/* Compute sine and cosine of Sun's declination */
sin_sdecl = sind(x: obl_ecl) * sind(x: slon);
cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
/* Compute the Sun's apparent radius, degrees */
sradius = 0.2666 / sr;
/* Do correction to upper limb, if necessary */
if upper_limb != 0 {
altit -= sradius
}
/* Compute the diurnal arc that the Sun traverses to reach */
/* the specified altitude altit: */
do {
let cost: Double = ( sind(x: altit) - sind(x: lat) * sin_sdecl ) / ( cosd(x: lat) * cos_sdecl );
if cost >= 1.0 {
t = 0.0 /* Sun always below altit */
}
else if cost <= -1.0 {
t = 24.0; /* Sun always above altit */
}
else {
t = (2.0/15.0) * acosd(x: cost); /* The diurnal arc, hours */
}
}
return t;
} /* __daylen__ */
/// This function computes the Sun's position at any instant.
private class func sunpos(d: Double) -> (lon: Double, r: Double ) {
/******************************************************/
/* Computes the Sun's ecliptic longitude and distance */
/* at an instant given in d, number of days since */
/* 2000 Jan 0.0. The Sun's ecliptic latitude is not */
/* computed, since it's always very near 0. */
/******************************************************/
var M: Double /* Mean anomaly of the Sun */
var w: Double /* Mean longitude of perihelion */
/* Note: Sun's mean longitude = M + w */
var e: Double /* Eccentricity of Earth's orbit */
var E: Double /* Eccentric anomaly */
var x: Double
var y: Double /* x, y coordinates in orbit */
var v: Double /* True anomaly */
/* Compute mean elements */
M = revolution( x: 356.0470 + 0.9856002585 * d );
w = 282.9404 + 4.70935E-5 * d;
e = 0.016709 - 1.151E-9 * d;
/* Compute true longitude and radius vector */
E = M + e * RADEG * sind(x: M) * ( 1.0 + e * cosd(x: M) );
x = cosd(x: E) - e;
y = sqrt( 1.0 - e*e ) * sind(x: E);
let r = sqrt( x*x + y*y ); /* Solar distance */
v = atan2d( y: y, x: x ); /* True anomaly */
var lon = v + w /* True solar longitude */
if lon >= 360.0 {
lon -= 360.0 /* Make it 0..360 degrees */
}
return (lon, r)
}
private class func sun_RA_dec( d: Double ) -> (RA: Double, dec: Double, r: Double ) {
/******************************************************/
/* Computes the Sun's equatorial coordinates RA, Decl */
/* and also its distance, at an instant given in d, */
/* the number of days since 2000 Jan 0.0. */
/******************************************************/
var obl_ecl: Double
var x: Double
var y: Double
var z: Double
/* Compute Sun's ecliptical coordinates */
let (lon, r) = sunpos( d: d )
/* Compute ecliptic rectangular coordinates (z=0) */
x = r * cosd(x: lon)
y = r * sind(x: lon)
/* Compute obliquity of ecliptic (inclination of Earth's axis) */
obl_ecl = 23.4393 - 3.563E-7 * d
/* Convert to equatorial rectangular coordinates - x is unchanged */
z = y * sind(x: obl_ecl)
y = y * cosd(x: obl_ecl)
/* Convert to spherical coordinates */
let RA = atan2d( y: y, x: x )
let dec = atan2d( y: z, x: sqrt(x*x + y*y) )
return (RA, dec, r)
} /* sun_RA_dec */
private static let INV360: Double = ( 1.0 / 360.0 )
/*******************************************************************/
/** This function reduces any angle to within the first revolution */
/** by subtracting or adding even multiples of 360.0 until the */
/** result is >= 0.0 and < 360.0 */
/*******************************************************************/
private class func revolution( x: Double ) -> Double {
/*****************************************/
/* Reduce angle to within 0..360 degrees */
/*****************************************/
return( x - 360.0 * floor( x * INV360 ) );
} /* revolution */
private class func rev180(x: Double ) -> Double {
/*********************************************/
/* Reduce angle to within +180..+180 degrees */
/*********************************************/
return( x - 360.0 * floor( x * INV360 + 0.5 ) );
} /* revolution */
/********************************************************************/
/** This function computes GMST0, the Greenwich Mean Sidereal Time */
/** at 0h UT (i.e. the sidereal time at the Greenwhich meridian at */
/** 0h UT). GMST is then the sidereal time at Greenwich at any */
/** time of the day. I've generalized GMST0 as well, and define it */
/** as: GMST0 = GMST - UT -- this allows GMST0 to be computed at */
/** other times than 0h UT as well. While this sounds somewhat */
/** contradictory, it is very practical: instead of computing */
/** GMST like: */
/** */
/** GMST = (GMST0) + UT * (366.2422/365.2422) */
/** */
/** where (GMST0) is the GMST last time UT was 0 hours, one simply */
/** computes: */
/** */
/** GMST = GMST0 + UT */
/** */
/** where GMST0 is the GMST "at 0h UT" but at the current moment! */
/** Defined in this way, GMST0 will increase with about 4 min a */
/** day. It also happens that GMST0 (in degrees, 1 hr = 15 degr) */
/** is equal to the Sun's mean longitude plus/minus 180 degrees! */
/** (if we neglect aberration, which amounts to 20 seconds of arc */
/** or 1.33 seconds of time) */
/** */
/********************************************************************/
private class func GMST0( d: Double ) -> Double {
var sidtim0: Double
/* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr */
/* L = M + w, as defined in sunpos(). Since I'm too lazy to */
/* add these numbers, I'll let the C compiler do it for me. */
/* Any decent C compiler will add the constants at compile */
/* time, imposing no runtime or code overhead. */
sidtim0 = revolution( x: ( 180.0 + 356.0470 + 282.9404 ) +
( 0.9856002585 + 4.70935E-5 ) * d )
return sidtim0
} /* GMST0 */
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment