-
-
Save BurgerZ/62cf95504c6004ecd442f66f49e34660 to your computer and use it in GitHub Desktop.
Extended moon phase utils
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/** Extended moon phase utils | |
Ported from https://github.com/solarissmoke/php-moon-phase/blob/master/Solaris/MoonPhase.php | |
Maksym Huk 2017 | |
*/ | |
import Foundation | |
import SwiftMoment | |
class Moon { | |
enum Phase : Int { | |
case newMoon | |
case firstQuarter | |
case fullMoon | |
case lastQuarter | |
case nextNewMoon | |
case nextFirstQuarter | |
case nextFullMoon | |
case nextLastQuarter | |
} | |
var timestamp: Double = 0 | |
var currentPhase: Double = 0 // Phase (0 to 1) | |
var illumination: Double = 0 // Illuminated fraction (0 to 1) | |
var age: Double = 0 // Age of moon (days) | |
var distance: Double = 0 // Distance (kilometres) | |
var angularDiameter: Double = 0 // Angular diameter (degrees) | |
var sunAngularDiameter: Double = 0 // Sun angular diameter (degrees) | |
var sunDistance: Double = 0 // Distance to Sun (kilometres) | |
var synodicMonth: Double = 0 // Synodic month (new Moon to new Moon) | |
private var quarters: [Date] = [] | |
init(_ date: Date) { | |
var pdate = date.timeIntervalSince1970 | |
self.timestamp = pdate | |
/* Astronomical constants */ | |
let epoch = 2444238.5 // 1980 January 0.0 | |
/* Constants defining the Sun's apparent orbit */ | |
let elonge = 278.833540 // Ecliptic longitude of the Sun at epoch 1980.0 | |
let elongp = 282.596403 // Ecliptic longitude of the Sun at perigee | |
let eccent = 0.016718 // Eccentricity of Earth's orbit | |
let sunsmax = 1.495985e8 // Semi-major axis of Earth's orbit, km | |
let sunangsiz = 0.533128 // Sun's angular size, degrees, at semi-major axis distance | |
/* Elements of the Moon's orbit, epoch 1980.0 */ | |
let mmlong = 64.975464 // Moon's mean longitude at the epoch | |
let mmlongp = 349.383063 // Mean longitude of the perigee at the epoch | |
let mecc = 0.054900 // Eccentricity of the Moon's orbit | |
let mangsiz = 0.5181 // Moon's angular size at distance a from Earth | |
let msmax = 384401.0 // Semi-major axis of Moon's orbit in km | |
let synmonth = 29.53058868 // Synodic month (new Moon to new Moon) | |
// pdate is coming in as a UNIX timstamp, so convert it to Julian | |
pdate = pdate / 86400 + 2440587.5 | |
/* Calculation of the Sun's position */ | |
let Day = pdate - epoch // Date within epoch | |
let N = fixangle((360 / 365.2422) * Day) // Mean anomaly of the Sun | |
let M = fixangle(N + elonge - elongp) // Convert from perigee co-ordinates to epoch 1980.0 | |
var Ec = kepler(M, eccent) // Solve equation of Kepler | |
Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2) | |
Ec = 2 * radToDeg(atan(Ec)) // True anomaly | |
let Lambdasun = fixangle(Ec + elongp) // Sun's geocentric ecliptic longitude | |
let F = ((1 + eccent * cos(degToRad(Ec))) / (1 - eccent * eccent)); // Orbital distance factor | |
let SunDist = sunsmax / F // Distance to Sun in km | |
let SunAng = F * sunangsiz // Sun's angular size in degrees | |
/* Calculation of the Moon's position */ | |
let ml = fixangle(13.1763966 * Day + mmlong) // Moon's mean longitude | |
let MM = fixangle(ml - 0.1114041 * Day - mmlongp) // Moon's mean anomaly | |
let Ev = 1.2739 * sin(degToRad(2 * (ml - Lambdasun) - MM)) // Evection | |
let Ae = 0.1858 * sin(degToRad(M)) // Annual equation | |
let A3 = 0.37 * sin(degToRad(M)) // Correction term | |
let MmP = MM + Ev - Ae - A3 // Corrected anomaly | |
let mEc = 6.2886 * sin(degToRad(MmP)) // Correction for the equation of the centre | |
let A4 = 0.214 * sin(degToRad(2 * MmP)) // Another correction term | |
let lP = ml + Ev + mEc - Ae + A4 // Corrected longitude | |
let V = 0.6583 * sin(degToRad(2 * (lP - Lambdasun))) // Variation | |
let lPP = lP + V // True longitude | |
/* Calculation of the phase of the Moon */ | |
let MoonAge = lPP - Lambdasun // Age of the Moon in degrees | |
let MoonPhase = (1 - cos(degToRad(MoonAge))) / 2 // Phase of the Moon | |
// Distance of moon from the centre of the Earth | |
let MoonDistPart1 = (msmax * (1 - mecc * mecc)) | |
let MoonDistPart2 = (1 + mecc * cos(degToRad(MmP + mEc))) | |
let MoonDist = MoonDistPart1 / MoonDistPart2 | |
let MoonDFrac = MoonDist / msmax | |
let MoonAng = mangsiz / MoonDFrac // Moon's angular diameter | |
// store results | |
self.synodicMonth = synmonth | |
self.currentPhase = fixangle(MoonAge) / 360 // Phase (0 to 1) | |
self.illumination = MoonPhase // Illuminated fraction (0 to 1) | |
self.age = synmonth * currentPhase // Age of moon (days) | |
self.distance = MoonDist // Distance (kilometres) | |
self.angularDiameter = MoonAng // Angular diameter (degrees) | |
self.sunDistance = SunDist // Distance to Sun (kilometres) | |
self.sunAngularDiameter = SunAng // Sun's angular diameter (degrees) | |
} | |
private func degToRad(_ x: Double) -> Double { | |
return x * M_PI / 180.0 | |
} | |
private func radToDeg(_ x: Double) -> Double { | |
return x * 180 / M_PI | |
} | |
private func fixangle(_ a: Double) -> Double { | |
return ( a - 360 * floor(a / 360) ) | |
} | |
private func kepler(_ m: Double, _ ecc: Double) -> Double { // Solve the equation of Kepler. | |
let epsilon = 0.000001 // 1E-6 | |
let m = degToRad(m) | |
var e = m | |
var delta: Double | |
repeat { | |
delta = e - ecc * sin(e) - m | |
e -= delta / ( 1 - ecc * cos(e) ) | |
} while abs(delta) > epsilon | |
return e | |
} | |
/** Calculates time of the mean new Moon for a given | |
base date. This argument K to this func is the | |
precomputed synodic month index, given by: | |
K = (year - 1900) * 12.3685 | |
where year is expressed as a year and fractional year. | |
*/ | |
private func meanPhase(_ sdate: Double, _ k: Double) -> Double { | |
// Time in Julian centuries from 1900 January 0.5 | |
let t = ( sdate - 2415020.0 ) / 36525 | |
let t2 = t * t | |
let t3 = t2 * t | |
let nt1 = 2415020.75933 + synodicMonth * k | |
+ 0.0001178 * t2 | |
- 0.000000155 * t3 | |
+ 0.00033 * sin( degToRad( 166.56 + 132.87 * t - 0.009173 * t2 ) ) | |
return nt1 | |
} | |
/** Given a K value used to determine the mean phase of | |
the new moon, and a phase selector (0.0, 0.25, 0.5, | |
0.75), obtain the true, corrected phase time. | |
*/ | |
private func truePhase(_ k: Double, _ phase: Double) -> Double { | |
var apcor = false | |
let k = k + phase // Add phase to new moon time | |
let t = k / 1236.85 // Time in Julian centuries from 1900 January 0.5 | |
let t2 = t * t // Square for frequent use | |
let t3 = t2 * t // Cube for frequent use | |
var pt = 2415020.75933 // Mean time of phase | |
+ synodicMonth * k | |
+ 0.0001178 * t2 | |
- 0.000000155 * t3 | |
+ 0.00033 * sin( degToRad( 166.56 + 132.87 * t - 0.009173 * t2 ) ) | |
let m = 359.2242 + 29.10535608 * k - 0.0000333 * t2 - 0.00000347 * t3 // Sun's mean anomaly | |
let mprime = 306.0253 + 385.81691806 * k + 0.0107306 * t2 + 0.00001236 * t3 // Moon's mean anomaly | |
let f = 21.2964 + 390.67050646 * k - 0.0016528 * t2 - 0.00000239 * t3 // Moon's argument of latitude | |
if ( phase < 0.01 || abs( phase - 0.5 ) < 0.01 ) { | |
// Corrections for New and Full Moon | |
pt += (0.1734 - 0.000393 * t) * sin( degToRad( m ) ) | |
+ 0.0021 * sin( degToRad( 2 * m ) ) | |
- 0.4068 * sin( degToRad( mprime ) ) | |
+ 0.0161 * sin( degToRad( 2 * mprime) ) | |
- 0.0004 * sin( degToRad( 3 * mprime ) ) | |
+ 0.0104 * sin( degToRad( 2 * f ) ) | |
- 0.0051 * sin( degToRad( m + mprime ) ) | |
- 0.0074 * sin( degToRad( m - mprime ) ) | |
+ 0.0004 * sin( degToRad( 2 * f + m ) ) | |
- 0.0004 * sin( degToRad( 2 * f - m ) ) | |
- 0.0006 * sin( degToRad( 2 * f + mprime ) ) | |
+ 0.0010 * sin( degToRad( 2 * f - mprime ) ) | |
+ 0.0005 * sin( degToRad( m + 2 * mprime ) ) | |
apcor = true | |
} else if ( abs( phase - 0.25 ) < 0.01 || abs( phase - 0.75 ) < 0.01 ) { | |
pt += (0.1721 - 0.0004 * t) * sin( degToRad( m ) ) | |
+ 0.0021 * sin( degToRad( 2 * m ) ) | |
- 0.6280 * sin( degToRad( mprime ) ) | |
+ 0.0089 * sin( degToRad( 2 * mprime) ) | |
- 0.0004 * sin( degToRad( 3 * mprime ) ) | |
+ 0.0079 * sin( degToRad( 2 * f ) ) | |
- 0.0119 * sin( degToRad( m + mprime ) ) | |
- 0.0047 * sin( degToRad ( m - mprime ) ) | |
+ 0.0003 * sin( degToRad( 2 * f + m ) ) | |
- 0.0004 * sin( degToRad( 2 * f - m ) ) | |
- 0.0006 * sin( degToRad( 2 * f + mprime ) ) | |
+ 0.0021 * sin( degToRad( 2 * f - mprime ) ) | |
+ 0.0003 * sin( degToRad( m + 2 * mprime ) ) | |
+ 0.0004 * sin( degToRad( m - 2 * mprime ) ) | |
- 0.0003 * sin( degToRad( 2 * m + mprime ) ) | |
if phase < 0.5 { // First quarter correction | |
pt += 0.0028 - 0.0004 * cos( degToRad( m ) ) + 0.0003 * cos( degToRad( mprime ) ) | |
} else { // Last quarter correction | |
pt += -0.0028 + 0.0004 * cos( degToRad( m ) ) - 0.0003 * cos( degToRad( mprime ) ) | |
} | |
apcor = true | |
} | |
if !apcor { // func was called with an invalid phase selector | |
return 0 | |
} | |
return pt | |
} | |
/** Find time of phases of the moon which surround the current date. | |
Five phases are found, starting and | |
ending with the new moons which bound the current lunation. | |
*/ | |
private func phaseHunt() { | |
let sdate = utcToJulian( timestamp ) | |
var adate = sdate - 45 | |
let ats = timestamp - 86400 * 45 | |
let atsMoment = moment(Date(timeIntervalSince1970: ats)) | |
let yy = Double(atsMoment.year) | |
let mm = Double(atsMoment.month) | |
var k1 = floor( ( yy + ( ( mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 ) | |
var nt1 = meanPhase( adate, k1 ) | |
adate = nt1 | |
var k2: Double = 0 | |
var nt2: Double = 0 | |
while (true) { | |
adate += synodicMonth | |
k2 = k1 + 1 | |
nt2 = meanPhase( adate, k2 ) | |
// if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate | |
if( abs( nt2 - sdate ) < 0.75 ) { | |
nt2 = truePhase( k2, 0.0 ) | |
} | |
if ( nt1 <= sdate && nt2 > sdate ) { | |
break | |
} | |
nt1 = nt2 | |
k1 = k2 | |
} | |
// results in Julian dates | |
let data = [ | |
truePhase( k1, 0.0 ), | |
truePhase( k1, 0.25 ), | |
truePhase( k1, 0.5 ), | |
truePhase( k1, 0.75 ), | |
truePhase( k2, 0.0 ), | |
truePhase( k2, 0.25 ), | |
truePhase( k2, 0.5 ), | |
truePhase( k2, 0.75 ) | |
] | |
quarters = [] | |
for v in data { | |
quarters.append(Date(timeIntervalSince1970:(v - 2440587.5) * 86400)) // convert to UNIX time | |
} | |
} | |
/* Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction). */ | |
private func utcToJulian(_ ts: TimeInterval) -> Double { | |
return ts / 86400 + 2440587.5 | |
} | |
func phase(_ phase: Phase) -> Date { | |
if quarters.isEmpty { | |
phaseHunt() | |
} | |
return quarters[phase.rawValue] | |
} | |
func nextMajorPhase(after date: Date) -> (Phase, Date) { | |
let phases: [(Phase, Date)] = [ | |
(.newMoon, phase(.newMoon)), | |
(.nextNewMoon, phase(.nextNewMoon)), | |
(.fullMoon, phase(.fullMoon)), | |
(.nextFullMoon, phase(.nextFullMoon)), | |
] | |
var earliestPhaseIndex = 0 | |
for (index, phase) in phases.enumerated() { | |
let date = phase.1 | |
let earliestDate = phases[earliestPhaseIndex].1 | |
if date.compare(earliestDate) == .orderedAscending { | |
earliestPhaseIndex = index | |
} | |
} | |
return phases[earliestPhaseIndex] | |
} | |
func currentPhaseName() -> String { | |
let names = [ "New Moon", "Waxing Crescent", "First Quarter", "Waxing Gibbous", "Full Moon", "Waning Gibbous", "Third Quarter", "Waning Crescent", "New Moon" ] | |
// There are eight phases, evenly split. A "New Moon" occupies the 1/16th phases either side of phase = 0, and the rest follow from that. | |
return names[ Int(floor( ( currentPhase + 0.0625 ) * 8 )) ] | |
} | |
} | |
var age: Double = 0 // Age of moon (days) | |
var distance: Double = 0 // Distance (kilometres) | |
var angularDiameter: Double = 0 // Angular diameter (degrees) | |
var sunAngularDiameter: Double = 0 // Sun angular diameter (degrees) | |
var sunDistance: Double = 0 // Distance to Sun (kilometres) | |
var synodicMonth: Double = 0 // Synodic month (new Moon to new Moon) | |
private var quarters: [Date] = [] | |
init(_ date: Date) { | |
var pdate = date.timeIntervalSince1970 | |
self.timestamp = pdate | |
/* Astronomical constants */ | |
let epoch = 2444238.5 // 1980 January 0.0 | |
/* Constants defining the Sun's apparent orbit */ | |
let elonge = 278.833540 // Ecliptic longitude of the Sun at epoch 1980.0 | |
let elongp = 282.596403 // Ecliptic longitude of the Sun at perigee | |
let eccent = 0.016718 // Eccentricity of Earth's orbit | |
let sunsmax = 1.495985e8 // Semi-major axis of Earth's orbit, km | |
let sunangsiz = 0.533128 // Sun's angular size, degrees, at semi-major axis distance | |
/* Elements of the Moon's orbit, epoch 1980.0 */ | |
let mmlong = 64.975464 // Moon's mean longitude at the epoch | |
let mmlongp = 349.383063 // Mean longitude of the perigee at the epoch | |
let mecc = 0.054900 // Eccentricity of the Moon's orbit | |
let mangsiz = 0.5181 // Moon's angular size at distance a from Earth | |
let msmax = 384401.0 // Semi-major axis of Moon's orbit in km | |
let synmonth = 29.53058868 // Synodic month (new Moon to new Moon) | |
// pdate is coming in as a UNIX timstamp, so convert it to Julian | |
pdate = pdate / 86400 + 2440587.5 | |
/* Calculation of the Sun's position */ | |
let Day = pdate - epoch // Date within epoch | |
let N = fixangle((360 / 365.2422) * Day) // Mean anomaly of the Sun | |
let M = fixangle(N + elonge - elongp) // Convert from perigee co-ordinates to epoch 1980.0 | |
var Ec = kepler(M, eccent) // Solve equation of Kepler | |
Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2) | |
Ec = 2 * radToDeg(atan(Ec)) // True anomaly | |
let Lambdasun = fixangle(Ec + elongp) // Sun's geocentric ecliptic longitude | |
let F = ((1 + eccent * cos(degToRad(Ec))) / (1 - eccent * eccent)); // Orbital distance factor | |
let SunDist = sunsmax / F // Distance to Sun in km | |
let SunAng = F * sunangsiz // Sun's angular size in degrees | |
/* Calculation of the Moon's position */ | |
let ml = fixangle(13.1763966 * Day + mmlong) // Moon's mean longitude | |
let MM = fixangle(ml - 0.1114041 * Day - mmlongp) // Moon's mean anomaly | |
let Ev = 1.2739 * sin(degToRad(2 * (ml - Lambdasun) - MM)) // Evection | |
let Ae = 0.1858 * sin(degToRad(M)) // Annual equation | |
let A3 = 0.37 * sin(degToRad(M)) // Correction term | |
let MmP = MM + Ev - Ae - A3 // Corrected anomaly | |
let mEc = 6.2886 * sin(degToRad(MmP)) // Correction for the equation of the centre | |
let A4 = 0.214 * sin(degToRad(2 * MmP)) // Another correction term | |
let lP = ml + Ev + mEc - Ae + A4 // Corrected longitude | |
let V = 0.6583 * sin(degToRad(2 * (lP - Lambdasun))) // Variation | |
let lPP = lP + V // True longitude | |
/* Calculation of the phase of the Moon */ | |
let MoonAge = lPP - Lambdasun // Age of the Moon in degrees | |
let MoonPhase = (1 - cos(degToRad(MoonAge))) / 2 // Phase of the Moon | |
// Distance of moon from the centre of the Earth | |
let MoonDistPart1 = (msmax * (1 - mecc * mecc)) | |
let MoonDistPart2 = (1 + mecc * cos(degToRad(MmP + mEc))) | |
let MoonDist = MoonDistPart1 / MoonDistPart2 | |
let MoonDFrac = MoonDist / msmax | |
let MoonAng = mangsiz / MoonDFrac // Moon's angular diameter | |
// store results | |
self.synodicMonth = synmonth | |
self.currentPhase = fixangle(MoonAge) / 360 // Phase (0 to 1) | |
self.illumination = MoonPhase // Illuminated fraction (0 to 1) | |
self.age = synmonth * currentPhase // Age of moon (days) | |
self.distance = MoonDist // Distance (kilometres) | |
self.angularDiameter = MoonAng // Angular diameter (degrees) | |
self.sunDistance = SunDist // Distance to Sun (kilometres) | |
self.sunAngularDiameter = SunAng // Sun's angular diameter (degrees) | |
} | |
private func degToRad(_ x: Double) -> Double { | |
return x * M_PI / 180.0 | |
} | |
private func radToDeg(_ x: Double) -> Double { | |
return x * 180 / M_PI | |
} | |
private func fixangle(_ a: Double) -> Double { | |
return ( a - 360 * floor(a / 360) ) | |
} | |
private func kepler(_ m: Double, _ ecc: Double) -> Double { // Solve the equation of Kepler. | |
let epsilon = 0.000001 // 1E-6 | |
let m = degToRad(m) | |
var e = m | |
var delta: Double | |
repeat { | |
delta = e - ecc * sin(e) - m | |
e -= delta / ( 1 - ecc * cos(e) ) | |
} while abs(delta) > epsilon | |
return e | |
} | |
/** Calculates time of the mean new Moon for a given | |
base date. This argument K to this func is the | |
precomputed synodic month index, given by: | |
K = (year - 1900) * 12.3685 | |
where year is expressed as a year and fractional year. | |
*/ | |
func meanPhase(_ sdate: Double, _ k: Double) -> Double { | |
// Time in Julian centuries from 1900 January 0.5 | |
let t = ( sdate - 2415020.0 ) / 36525 | |
let t2 = t * t | |
let t3 = t2 * t | |
let nt1 = 2415020.75933 + synodicMonth * k | |
+ 0.0001178 * t2 | |
- 0.000000155 * t3 | |
+ 0.00033 * sin( degToRad( 166.56 + 132.87 * t - 0.009173 * t2 ) ) | |
return nt1 | |
} | |
/** Given a K value used to determine the mean phase of | |
the new moon, and a phase selector (0.0, 0.25, 0.5, | |
0.75), obtain the true, corrected phase time. | |
*/ | |
func truePhase(_ k: Double, _ phase: Double) -> Double { | |
var apcor = false | |
let k = k + phase // Add phase to new moon time | |
let t = k / 1236.85 // Time in Julian centuries from 1900 January 0.5 | |
let t2 = t * t // Square for frequent use | |
let t3 = t2 * t // Cube for frequent use | |
var pt = 2415020.75933 // Mean time of phase | |
+ synodicMonth * k | |
+ 0.0001178 * t2 | |
- 0.000000155 * t3 | |
+ 0.00033 * sin( degToRad( 166.56 + 132.87 * t - 0.009173 * t2 ) ) | |
let m = 359.2242 + 29.10535608 * k - 0.0000333 * t2 - 0.00000347 * t3 // Sun's mean anomaly | |
let mprime = 306.0253 + 385.81691806 * k + 0.0107306 * t2 + 0.00001236 * t3 // Moon's mean anomaly | |
let f = 21.2964 + 390.67050646 * k - 0.0016528 * t2 - 0.00000239 * t3 // Moon's argument of latitude | |
if ( phase < 0.01 || abs( phase - 0.5 ) < 0.01 ) { | |
// Corrections for New and Full Moon | |
pt += (0.1734 - 0.000393 * t) * sin( degToRad( m ) ) | |
+ 0.0021 * sin( degToRad( 2 * m ) ) | |
- 0.4068 * sin( degToRad( mprime ) ) | |
+ 0.0161 * sin( degToRad( 2 * mprime) ) | |
- 0.0004 * sin( degToRad( 3 * mprime ) ) | |
+ 0.0104 * sin( degToRad( 2 * f ) ) | |
- 0.0051 * sin( degToRad( m + mprime ) ) | |
- 0.0074 * sin( degToRad( m - mprime ) ) | |
+ 0.0004 * sin( degToRad( 2 * f + m ) ) | |
- 0.0004 * sin( degToRad( 2 * f - m ) ) | |
- 0.0006 * sin( degToRad( 2 * f + mprime ) ) | |
+ 0.0010 * sin( degToRad( 2 * f - mprime ) ) | |
+ 0.0005 * sin( degToRad( m + 2 * mprime ) ) | |
apcor = true | |
} else if ( abs( phase - 0.25 ) < 0.01 || abs( phase - 0.75 ) < 0.01 ) { | |
pt += (0.1721 - 0.0004 * t) * sin( degToRad( m ) ) | |
+ 0.0021 * sin( degToRad( 2 * m ) ) | |
- 0.6280 * sin( degToRad( mprime ) ) | |
+ 0.0089 * sin( degToRad( 2 * mprime) ) | |
- 0.0004 * sin( degToRad( 3 * mprime ) ) | |
+ 0.0079 * sin( degToRad( 2 * f ) ) | |
- 0.0119 * sin( degToRad( m + mprime ) ) | |
- 0.0047 * sin( degToRad ( m - mprime ) ) | |
+ 0.0003 * sin( degToRad( 2 * f + m ) ) | |
- 0.0004 * sin( degToRad( 2 * f - m ) ) | |
- 0.0006 * sin( degToRad( 2 * f + mprime ) ) | |
+ 0.0021 * sin( degToRad( 2 * f - mprime ) ) | |
+ 0.0003 * sin( degToRad( m + 2 * mprime ) ) | |
+ 0.0004 * sin( degToRad( m - 2 * mprime ) ) | |
- 0.0003 * sin( degToRad( 2 * m + mprime ) ) | |
if phase < 0.5 { // First quarter correction | |
pt += 0.0028 - 0.0004 * cos( degToRad( m ) ) + 0.0003 * cos( degToRad( mprime ) ) | |
} else { // Last quarter correction | |
pt += -0.0028 + 0.0004 * cos( degToRad( m ) ) - 0.0003 * cos( degToRad( mprime ) ) | |
} | |
apcor = true | |
} | |
if !apcor { // func was called with an invalid phase selector | |
return 0 | |
} | |
return pt | |
} | |
/** Find time of phases of the moon which surround the current date. | |
Five phases are found, starting and | |
ending with the new moons which bound the current lunation. | |
*/ | |
func phaseHunt() { | |
let sdate = utcToJulian( timestamp ) | |
var adate = sdate - 45 | |
let ats = timestamp - 86400 * 45 | |
let atsMoment = moment(Date(timeIntervalSince1970: ats)) | |
let yy = Double(atsMoment.year) | |
let mm = Double(atsMoment.month) | |
var k1 = floor( ( yy + ( ( mm - 1 ) * ( 1 / 12 ) ) - 1900 ) * 12.3685 ) | |
var nt1 = meanPhase( adate, k1 ) | |
adate = nt1 | |
var k2: Double = 0 | |
var nt2: Double = 0 | |
while (true) { | |
adate += synodicMonth | |
k2 = k1 + 1 | |
nt2 = meanPhase( adate, k2 ) | |
// if nt2 is close to sdate, then mean phase isn't good enough, we have to be more accurate | |
if( abs( nt2 - sdate ) < 0.75 ) { | |
nt2 = truePhase( k2, 0.0 ) | |
} | |
if ( nt1 <= sdate && nt2 > sdate ) { | |
break | |
} | |
nt1 = nt2 | |
k1 = k2 | |
} | |
// results in Julian dates | |
let data = [ | |
truePhase( k1, 0.0 ), | |
truePhase( k1, 0.25 ), | |
truePhase( k1, 0.5 ), | |
truePhase( k1, 0.75 ), | |
truePhase( k2, 0.0 ), | |
truePhase( k2, 0.25 ), | |
truePhase( k2, 0.5 ), | |
truePhase( k2, 0.75 ) | |
] | |
quarters = [] | |
for v in data { | |
quarters.append(Date(timeIntervalSince1970:(v - 2440587.5) * 86400)) // convert to UNIX time | |
} | |
} | |
/* Convert UNIX timestamp to astronomical Julian time (i.e. Julian date plus day fraction). */ | |
func utcToJulian(_ ts: TimeInterval) -> Double { | |
return ts / 86400 + 2440587.5 | |
} | |
func phase(_ phase: Phase) -> Date { | |
if quarters.isEmpty { | |
phaseHunt() | |
} | |
return quarters[phase.rawValue] | |
} | |
func nextMajorPhase(after date: Date) -> (Phase, Date) { | |
let phases: [(Phase, Date)] = [ | |
(.newMoon, phase(.newMoon)), | |
(.nextNewMoon, phase(.nextNewMoon)), | |
(.fullMoon, phase(.fullMoon)), | |
(.nextFullMoon, phase(.nextFullMoon)), | |
] | |
var earliestPhaseIndex = 0 | |
for (index, phase) in phases.enumerated() { | |
let date = phase.1 | |
let earliestDate = phases[earliestPhaseIndex].1 | |
if date.compare(earliestDate) == .orderedAscending { | |
earliestPhaseIndex = index | |
} | |
} | |
return phases[earliestPhaseIndex] | |
} | |
func currentPhaseName() -> String { | |
let names = [ "New Moon", "Waxing Crescent", "First Quarter", "Waxing Gibbous", "Full Moon", "Waning Gibbous", "Third Quarter", "Waning Crescent", "New Moon" ] | |
// There are eight phases, evenly split. A "New Moon" occupies the 1/16th phases either side of phase = 0, and the rest follow from that. | |
return names[ Int(floor( ( currentPhase + 0.0625 ) * 8 )) ] | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment