diff --git a/astrotime.go b/astrotime.go index f2164a7..e60f442 100755 --- a/astrotime.go +++ b/astrotime.go @@ -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 (petegray@ieee.org), 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 (jon@wildducktheories.com) to allow dawn +// and dusk calculations. +// import ( "math" "time" @@ -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 @@ -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 @@ -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) @@ -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 @@ -299,22 +314,41 @@ 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 @@ -322,14 +356,17 @@ func CalcSunrise(t time.Time, latitude float64, longitude float64) time.Time { // 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 } @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 +} diff --git a/astrotime_test.go b/astrotime_test.go new file mode 100644 index 0000000..cd811cf --- /dev/null +++ b/astrotime_test.go @@ -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)) + } + + } +}