Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/dawn+dusk #4

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Feature/dawn+dusk #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
25806a0
Modification to allow dawn and dusk calculations.
jonseymour Nov 3, 2014
f70b855
Add test data for Sydney, Stockholm and New York.
jonseymour Nov 3, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 83 additions & 21 deletions astrotime.go
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,16 @@
// NAA - NOAA's Astronomical Algorithms
package astrotime

// (JavaScript web page
// (JavaScript web page
// http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html by
// Chris Cornwall, Aaron Horiuchi and Chris Lehman)
// Ported to C++ by Pete Gray ([email protected]), July 2006
// Released as Open Source and can be used in any way, as long as the
// Released as Open Source and can be used in any way, as long as the
// above description remains in place.

//
// Modified by Jon Seymour ([email protected]) to allow dawn
// and dusk calculations.
//
import (
"math"
"time"
Expand All @@ -26,6 +29,17 @@ const (
OneDay = time.Hour * 24
)

const (
ASTRONOMICAL_DAWN = -18.0
ASTRONOMICAL_DUSK = -18.0
NAUTICAL_DAWN = -12.0
NAUTICAL_DUSK = -12.0
CIVIL_DAWN = -6.0
CIVIL_DUSK = -6.0
SUNRISE = 0
SUNSET = 0
)

// CalcJD converts a time.Time object to a julian date
func CalcJD(t time.Time) float64 {
y, m, d, hh, mm, ss, ms := t.Year(), int(t.Month()), t.Day(), t.Hour(), t.Minute(), t.Second(), t.Nanosecond()/1e6
Expand Down Expand Up @@ -235,10 +249,11 @@ func calcSunDeclination(t float64) float64 {
//Return value:
// hour angle of sunrise in radians

func calcHourAngleSunrise(lat float64, solarDec float64) float64 {
func calcHourAngleSunrise(lat float64, solarDec float64, solarElevation float64) float64 {
latRad := DegToRad * lat
sdRad := DegToRad * solarDec
return (math.Acos(math.Cos(DegToRad*90.833)/(math.Cos(latRad)*math.Cos(sdRad)) - math.Tan(latRad)*math.Tan(sdRad)))
correction := calcRefractiveCorrection(solarElevation)
return (math.Acos(math.Cos(DegToRad*(90+correction))/(math.Cos(latRad)*math.Cos(sdRad)) - math.Tan(latRad)*math.Tan(sdRad)))
}

// Name: calcSolNoonUTC
Expand Down Expand Up @@ -273,7 +288,7 @@ func calcSolNoonUTC(t float64, longitude float64) float64 {
// time in minutes from zero Z

// Calculate the UTC sunrise for the given day at the given location
func calcSunriseUTC(jd float64, latitude float64, longitude float64) float64 {
func calcSunriseUTC(jd float64, latitude float64, longitude float64, solarElevation float64) float64 {

t := calcTimeJulianCent(jd)

Expand All @@ -288,7 +303,7 @@ func calcSunriseUTC(jd float64, latitude float64, longitude float64) float64 {

eqTime := calcEquationOfTime(tnoon)
solarDec := calcSunDeclination(tnoon)
hourAngle := calcHourAngleSunrise(latitude, solarDec)
hourAngle := calcHourAngleSunrise(latitude, solarDec, solarElevation)

delta := longitude - RadToDeg*hourAngle
timeDiff := 4 * delta
Expand All @@ -299,37 +314,59 @@ func calcSunriseUTC(jd float64, latitude float64, longitude float64) float64 {
newt := calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0)
eqTime = calcEquationOfTime(newt)
solarDec = calcSunDeclination(newt)
hourAngle = calcHourAngleSunrise(latitude, solarDec)
hourAngle = calcHourAngleSunrise(latitude, solarDec, solarElevation)
delta = longitude - RadToDeg*hourAngle
timeDiff = 4 * delta
timeUTC = 720 + timeDiff - eqTime
return timeUTC
}

// CalcSunrise calculates the sunrise, in local time, on the day t at the
// location specified in longitude and latitude.
func CalcSunrise(t time.Time, latitude float64, longitude float64) time.Time {
// CalcDawn calculates the dawn, in local time, on the day t at the
// location specified in longitude and latitude. If solarElevation is 0, then dawn is identical
// to sunrise. A negative solarElevation specifies that the time calculated is for when the sun
// is -solarElevation degrees below the horizon.
func CalcDawn(t time.Time, latitude float64, longitude float64, solarElevation float64) time.Time {
jd := CalcJD(t)
sunriseUTC := time.Duration(math.Floor(calcSunriseUTC(jd, latitude, longitude)*60) * 1e9)
sunriseUTC := time.Duration(math.Floor(calcSunriseUTC(jd, latitude, longitude, solarElevation)*60) * 1e9)
loc, _ := time.LoadLocation("UTC")
return time.Date(t.Year(), t.Month(), t.Day(), 0, 0, 0, 0, loc).Add(sunriseUTC).In(t.Location())
}

// CalcSunrise calculates the sunrise, in local time, on the day t at the
// location specified in longitude and latitude.
func CalcSunrise(t time.Time, latitude float64, longitude float64) time.Time {
return CalcDawn(t, latitude, longitude, SUNRISE)
}

// Calculates the refractive correction for the specified solar elevation.
// The refractive correction only applies at sunrise and sunset (solarElevation == 0).
// Otherwise the correction is simply the negative of the solarElevation
func calcRefractiveCorrection(solarElevation float64) float64 {
if solarElevation == 0 {
return 0.833
} else {
return -solarElevation
}
}

// Name: calcHourAngleSunset
// Type: Function
// Purpose: calculate the hour angle of the sun at sunset for the
// latitude
// Arguments:
// lat : latitude of observer in degrees
// solarDec : declination angle of sun in degrees
// solarElevation : solarElevation of the sun w.r.t. horizon - 0 for sunset, -6 for civil dusk, etc...
// Return value:
// hour angle of sunset in radians

func calcHourAngleSunset(lat float64, solarDec float64) float64 {
func calcHourAngleSunset(lat float64, solarDec float64, solarElevation float64) float64 {
latRad := DegToRad * lat
sdRad := DegToRad * solarDec

HA := (math.Acos(math.Cos(DegToRad*90.833)/(math.Cos(latRad)*math.Cos(sdRad)) - math.Tan(latRad)*math.Tan(sdRad)))
correction := calcRefractiveCorrection(solarElevation)

HA := (math.Acos(math.Cos(DegToRad*(90+correction))/(math.Cos(latRad)*math.Cos(sdRad)) - math.Tan(latRad)*math.Tan(sdRad)))

return -HA // in radians
}
Expand All @@ -345,7 +382,7 @@ func calcHourAngleSunset(lat float64, solarDec float64) float64 {
// Return value:
// time in minutes from zero Z

func calcSunsetUTC(jd float64, latitude float64, longitude float64) float64 {
func calcSunsetUTC(jd float64, latitude float64, longitude float64, solarElevation float64) float64 {
t := calcTimeJulianCent(jd)

// *** Find the time of solar noon at the location, and use
Expand All @@ -359,7 +396,7 @@ func calcSunsetUTC(jd float64, latitude float64, longitude float64) float64 {

eqTime := calcEquationOfTime(tnoon)
solarDec := calcSunDeclination(tnoon)
hourAngle := calcHourAngleSunset(latitude, solarDec)
hourAngle := calcHourAngleSunset(latitude, solarDec, solarElevation)

delta := longitude - RadToDeg*hourAngle
timeDiff := 4 * delta
Expand All @@ -370,22 +407,29 @@ func calcSunsetUTC(jd float64, latitude float64, longitude float64) float64 {
newt := calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0)
eqTime = calcEquationOfTime(newt)
solarDec = calcSunDeclination(newt)
hourAngle = calcHourAngleSunset(latitude, solarDec)
hourAngle = calcHourAngleSunset(latitude, solarDec, solarElevation)

delta = longitude - RadToDeg*hourAngle
timeDiff = 4 * delta
return 720 + timeDiff - eqTime
}

// CalcSunset calculates the sunset, in local time, on the day t at the
// location specified in longitude and latitude.
func CalcSunset(t time.Time, latitude float64, longitude float64) time.Time {
// CalcDusk calculates the dusk, in local time, on the day t at the
// location specified in longitude and latitude. The solarElevation specifies that the time that the sun will
// be -solarElevation degrees below the horizon.
func CalcDusk(t time.Time, latitude float64, longitude float64, solarElevation float64) time.Time {
jd := CalcJD(t)
sunsetUTC := time.Duration(math.Floor(calcSunsetUTC(jd, latitude, longitude)*60) * 1e9)
sunsetUTC := time.Duration(math.Floor(calcSunsetUTC(jd, latitude, longitude, solarElevation)*60) * 1e9)
loc, _ := time.LoadLocation("UTC")
return time.Date(t.Year(), t.Month(), t.Day(), 0, 0, 0, 0, loc).Add(sunsetUTC).In(t.Location())
}

// CalcSunset calculates the sunset, in local time, on the day t at the
// location specified in longitude and latitude.
func CalcSunset(t time.Time, latitude float64, longitude float64) time.Time {
return CalcDusk(t, latitude, longitude, SUNSET)
}

// NextSunrise returns date/time of the next sunrise after tAfter
func NextSunrise(tAfter time.Time, latitude float64, longitude float64) (tSunrise time.Time) {
tSunrise = CalcSunrise(tAfter, latitude, longitude)
Expand All @@ -395,6 +439,15 @@ func NextSunrise(tAfter time.Time, latitude float64, longitude float64) (tSunris
return
}

// NextDawn returns date/time of the next dawn at the specified solar solarElevation, after the specified time.
func NextDawn(tAfter time.Time, latitude float64, longitude float64, solarElevation float64) (tSunrise time.Time) {
tSunrise = CalcDawn(tAfter, latitude, longitude, solarElevation)
if tAfter.After(tSunrise) {
tSunrise = CalcDawn(tAfter.Add(OneDay), latitude, longitude, solarElevation)
}
return
}

// NextSunset returns date/time of the next sunset after tAfter
func NextSunset(tAfter time.Time, latitude float64, longitude float64) (tSunset time.Time) {
tSunset = CalcSunset(tAfter, latitude, longitude)
Expand All @@ -403,3 +456,12 @@ func NextSunset(tAfter time.Time, latitude float64, longitude float64) (tSunset
}
return
}

// NextDusk returns date/time of the next dusk at the specified solar solarElevation, after the specified tine.
func NextDusk(tAfter time.Time, latitude float64, longitude float64, solarElevation float64) (tSunset time.Time) {
tSunset = CalcDusk(tAfter, latitude, longitude, solarElevation)
if tAfter.After(tSunset) {
tSunset = CalcDusk(tAfter.Add(OneDay), latitude, longitude, solarElevation)
}
return
}
91 changes: 91 additions & 0 deletions astrotime_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
package astrotime

import (
"fmt"
"math"
"testing"
"time"
)

type location struct {
timezone string
latitude float64
longitude float64
}

type times struct {
date string
sunrise string
sunset string
dawn string
dusk string
}

type test struct {
times times
location location
error time.Duration
}

// test data from http://www.timeanddate.com/worldclock/sunrise.html

var (
SYDNEY = location{
timezone: "Australia/Sydney",
latitude: -33.86,
longitude: -151.20,
}
STOCKHOLM = location{
timezone: "Europe/Stockholm",
latitude: 59.33,
longitude: -18.067,
}
NEWYORK = location{
timezone: "America/New_York",
latitude: 40.642,
longitude: 74.017,
}
NOVEMBER = "2014-11-01"
JULY = "2015-07-01"
MINUTE = time.Minute
DATA = []*test{
&test{times{NOVEMBER, "05:55", "19:23", "05:29", "19:49"}, SYDNEY, MINUTE},
&test{times{NOVEMBER, "07:07", "15:55", "06:23", "16:39"}, STOCKHOLM, MINUTE},
&test{times{NOVEMBER, "07:26", "17:52", "06:58", "18:21"}, NEWYORK, MINUTE},
&test{times{JULY, "07:01", "16:57", "06:33", "17:25"}, SYDNEY, MINUTE},
&test{times{JULY, "03:37", "22:06", "02:09", "23:33"}, STOCKHOLM, MINUTE},
&test{times{JULY, "05:29", "20:31", "04:55", "21:04"}, NEWYORK, MINUTE},
}
)

func TestData(t *testing.T) {
for _, datum := range DATA {
if loc, err := time.LoadLocation(datum.location.timezone); err != nil {
t.Fatalf("bad location %s", datum.location.timezone)
} else {
formatted := fmt.Sprintf("%s 00:00:00", datum.times.date)
midnight, err := time.ParseInLocation("2006-01-02 15:04:05", formatted, loc)
if err != nil {
t.Fatalf("failed to parse date/timestamp '%s': %s", formatted, err)
}

checkDate := func(description string, timeOfDay string, calculated time.Time) {
formatted := fmt.Sprintf("%s %s", datum.times.date, timeOfDay)
timestamp, err := time.ParseInLocation("2006-01-02 15:04", formatted, loc)
if err != nil {
t.Fatalf("failed to parse date/timestamp '%s': %s", formatted, err)
}
actualError := math.Abs((float64)(timestamp.Sub(calculated)))
if actualError > float64(datum.error) {
t.Errorf("calculated -> %v, wanted -> %v %d -> (wanted < %d). location=%s date=%s", calculated, timestamp, actualError, datum.error, datum.location.timezone, datum.times.date)
}
}

checkDate("dawn", datum.times.dawn, NextDawn(midnight, datum.location.latitude, datum.location.longitude, CIVIL_DAWN))
checkDate("sunrise", datum.times.sunrise, NextSunrise(midnight, datum.location.latitude, datum.location.longitude))
checkDate("sunset", datum.times.sunset, NextSunset(midnight, datum.location.latitude, datum.location.longitude))
checkDate("dusk", datum.times.dusk, NextDusk(midnight, datum.location.latitude, datum.location.longitude, CIVIL_DUSK))
}

}
}