calendar: Add sunrise, sunset, solar-noon

Closes #2510

diff --git a/basis/calendar/calendar-docs.factor b/basis/calendar/calendar-docs.factor
index 2e79a60a9e512574e4ac9c4f63e24e339624d4fb..87d5197e839d93d5c10d4b0b73eb635076115e53 100644 (file)
--- a/basis/calendar/calendar-docs.factor
+++ b/basis/calendar/calendar-docs.factor
@@ -493,6 +493,29 @@ HELP: since-1970
     { "timestamp" timestamp } }
 { $description "Adds the duration to the beginning of Unix time and returns the result as a timestamp." } ;
 
+HELP: sunrise
+{ $values
+    { "timestamp" timestamp }
+    { "latitude" real }
+    { "longitude" real }
+    { "new-timestamp" timestamp } }
+{ $description "Calculates the time of sunrise on the given day at the given location in the given timezone." } ;
+
+HELP: sunset
+{ $values
+    { "timestamp" timestamp }
+    { "latitude" real }
+    { "longitude" real }
+    { "new-timestamp" timestamp } }
+{ $description "Calculates the time of sunset on the given day at the given location in the given timezone." } ;
+
+HELP: solar-noon
+{ $values
+    { "timestamp" timestamp }
+    { "longitude" real }
+    { "new-timestamp" timestamp } }
+{ $description "Calculates solar noon of the given day at the given longitude in the given timezone." } ;
+
 ARTICLE: "calendar" "Calendar"
 "The " { $vocab-link "calendar" } " vocabulary defines two data types and a set of operations on them:"
 { $subsections
@@ -526,6 +549,13 @@ ARTICLE: "calendar" "Calendar"
 }
 "Both " { $link timestamp } "s and " { $link duration } "s implement the " { $link "math.order" } "."
 $nl
+"Solar position calculations:"
+{ $subsections
+    sunrise
+    sunset
+    solar-noon
+}
+$nl
 "Metadata about the calendar:"
 { $subsections "calendar-facts" } ;
 

diff --git a/basis/calendar/calendar-tests.factor b/basis/calendar/calendar-tests.factor
index 5990f83769c5401170066165e57aa7f335bf0b5c..0cac2525195236272e3b81702cad528fc7a9a7fa 100644 (file)
--- a/basis/calendar/calendar-tests.factor
+++ b/basis/calendar/calendar-tests.factor
@@ -443,3 +443,22 @@ IN: calendar
         { second 59+999/1000 }
     }
 } [ 2023 4 13 <date-gmt> start-of-week dup end-of-week ] unit-test
+
+{
+    T{ timestamp
+        { year 2024 }
+        { month 2 }
+        { day 15 }
+        { hour 7 }
+    }
+} [ 2024 02 15 <date-gmt> 48.87 2.67 sunrise >gmt ] unit-test
+
+{
+    T{ timestamp
+        { year 2024 }
+        { month 2 }
+        { day 15 }
+        { hour 17 }
+        { minute 7 }
+    }
+} [ 2024 02 15 <date-gmt> 48.87 2.67 sunset >gmt ] unit-test

diff --git a/basis/calendar/calendar.factor b/basis/calendar/calendar.factor
index 1682001e6b6b3a7ae2bbe913a480ec13c603153f..0b28bdbf18eb484c448f3e11036169a2df667b63 100644 (file)
--- a/basis/calendar/calendar.factor
+++ b/basis/calendar/calendar.factor
@@ -1,9 +1,9 @@
 ! Copyright (C) 2007 Doug Coleman.
 ! See https://factorcode.org/license.txt for BSD license.
 USING: accessors arrays classes.tuple combinators
-combinators.short-circuit kernel literals math math.functions
-math.intervals math.order math.statistics sequences slots.syntax
-system vocabs vocabs.loader ;
+combinators.short-circuit kernel literals math math.constants
+math.functions math.intervals math.order math.statistics
+sequences slots.syntax system vocabs vocabs.loader ;
 FROM: ranges => [a..b) ;
 IN: calendar
 
@@ -42,6 +42,7 @@ TUPLE: timestamp
     { gmt-offset duration } ;
 
 <PRIVATE
+
 <<
 CONSTANT: day-counts { 0 31 28 31 30 31 30 31 31 30 31 30 31 }
 >>
@@ -830,6 +831,66 @@ M: integer weeks-in-week-year
 M: timestamp weeks-in-week-year
     { [ january 1 >>day thursday? ] [ december 31 >>day thursday? ] } 1|| 53 52 ? ;
 
+! https://gml.noaa.gov/grad/solcalc/solareqns.PDF
+
+<PRIVATE
+
+: fractional-year ( timestamp -- radians )
+    [ days-in-year 2pi swap / ]
+    [ day-of-year 1 - ]
+    [ hour>> 12 - 24 / + * ] tri ;
+
+:: declination ( timestamp -- radians )
+    timestamp fractional-year :> γ
+    0.006918
+    0.399912 γ cos * -
+    0.070257 γ sin * +
+    0.006758 γ 2 * cos * -
+    0.000907 γ 2 * sin * +
+    0.002697 γ 3 * cos * -
+    0.00148  γ 3 * sin * + ;
+
+:: hour-angle ( timestamp latitude -- degrees )
+    timestamp declination :> decl
+    latitude deg>rad :> lat
+    90.833 deg>rad cos
+    lat cos decl cos * /
+    lat tan decl tan * -
+    acos rad>deg ;
+
+:: equation-of-time ( timestamp -- minutes )
+    timestamp fractional-year :> γ
+    0.000075
+    0.001868 γ cos * +
+    0.032077 γ sin * -
+    0.014615 γ 2 * cos * -
+    0.040849 γ 2 * sin * -
+    229.18 * ;
+
+: preserve-gmt-offset ( timestamp quot -- timestamp' )
+    '[ >utc @ ] [ gmt-offset>> convert-timezone ] bi ; inline
+
+: (sunrise/sunset) ( timestamp latitude longitude quot -- new-timestamp )
+    '[
+        [ noon ]
+        [ _ hour-angle _ swap @ 4 * ]
+        [ equation-of-time ] tri + round >integer minutes time-
+    ] preserve-gmt-offset ; inline
+
+PRIVATE>
+
+: sunrise ( timestamp latitude longitude -- new-timestamp )
+    [ + ] (sunrise/sunset) ;
+
+: sunset ( timestamp latitude longitude -- new-timestamp )
+    [ - ] (sunrise/sunset) ;
+
+: solar-noon ( timestamp longitude -- new-timestamp )
+    '[
+        [ noon _ 4 * ] [ equation-of-time ] bi + minutes time-
+        [ round >integer ] change-second
+    ] preserve-gmt-offset ;
+
 {
     { [ os unix? ] [ "calendar.unix" ] }
     { [ os windows? ] [ "calendar.windows" ] }