Skip to content

Commit c39a4de

Browse files
authored
Merge pull request #57 from JuliaAstro/sun-set-rise
Add sunrise / sunset calculation using SPA algorithm
2 parents 486727a + c712527 commit c39a4de

36 files changed

+1725
-404
lines changed

Project.toml

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "SolarPosition"
22
uuid = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"
3-
version = "0.3.2"
3+
version = "0.4.0"
44
authors = ["Stefan de Lange"]
55

66
[workspace]
@@ -9,6 +9,7 @@ projects = ["test", "docs", "benchmark", "examples"]
99
[deps]
1010
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
1111
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
12+
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
1213
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
1314
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
1415
TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
@@ -31,6 +32,7 @@ DocStringExtensions = "0.8, 0.9"
3132
Makie = "0.24"
3233
ModelingToolkit = "11"
3334
OhMyThreads = "0.8"
35+
Reexport = "1"
3436
StructArrays = "0.7"
3537
Symbolics = "6,7"
3638
Tables = "1"

docs/make.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ makedocs(;
3030
format = Documenter.HTML(;
3131
canonical = "https://juliaastro.org/SolarPosition/stable/",
3232
size_threshold = 2^20, # 1 MB
33+
assets = String["assets/citations.css"],
3334
),
3435
plugins = [bib],
3536
pages = [
@@ -45,6 +46,7 @@ makedocs(;
4546
"reference.md",
4647
"positioning.md",
4748
"refraction.md",
49+
"utilities.md",
4850
"deltat.md",
4951
"literature.md",
5052
"contributing.md",

docs/src/assets/citations.css

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
.citation dl {
2+
display: grid;
3+
grid-template-columns: max-content auto; }
4+
.citation dt {
5+
grid-column-start: 1; }
6+
.citation dd {
7+
grid-column-start: 2;
8+
margin-bottom: 0.75em; }
9+
.citation ul {
10+
padding: 0 0 2.25em 0;
11+
margin: 0;
12+
list-style: none !important;}
13+
.citation ul li {
14+
text-indent: -2.25em;
15+
margin: 0.33em 0.5em 0.5em 2.25em;}
16+
.citation ol li {
17+
padding-left:0.75em;}

docs/src/guides/new-algorithm.md

Lines changed: 2 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -158,28 +158,10 @@ Edit `src/Positioning/Positioning.jl` to include your new file:
158158

159159
```julia
160160
# Near the bottom of the file, with other includes
161-
include("utils.jl")
162-
include("deltat.jl")
163-
include("psa.jl")
164-
include("noaa.jl")
165-
include("walraven.jl")
166-
include("usno.jl")
167-
include("spa.jl")
168161
include("simple.jl") # Add your new file
169162

170-
# Add to the export list
171-
export Observer,
172-
PSA,
173-
NOAA,
174-
Walraven,
175-
USNO,
176-
SPA,
177-
SimpleAlgorithm, # Add your algorithm
178-
solar_position,
179-
solar_position!,
180-
SolPos,
181-
ApparentSolPos,
182-
SPASolPos
163+
# Add your algorithm to the export list
164+
export SimpleAlgorithm
183165
```
184166

185167
### 4.2 Export from Main Module

docs/src/reference.md

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,6 @@ SolarPosition.Positioning.solar_position!
2727
SolarPosition.Positioning.Observer
2828
SolarPosition.Positioning.SolPos
2929
SolarPosition.Positioning.ApparentSolPos
30-
SolarPosition.Positioning.SPASolPos
3130
SolarPosition.Positioning.SPAObserver
3231
```
3332

docs/src/utilities.md

Lines changed: 298 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,298 @@
1+
# [Utilities](@id utilities)
2+
3+
`SolarPosition.jl` provides several utility functions for common solar position-related
4+
calculations, such as determining solar transit, sunrise, and sunset times for a given
5+
location and date.
6+
7+
For now, only the analytical [SPA](@ref spa-algorithm) algorithm is supported for these
8+
utility functions. This algorithm was developed by [Jean Meeus](https://en.wikipedia.org/wiki/Jean_Meeus)
9+
in his book *Astronomical Algorithms* [MEEUS91](@cite) and is widely used in the
10+
astronomical community for solar position calculations, particularly because it is
11+
relatively accurate but computationally efficient and simple to implement.
12+
13+
## Sunrise, Sunset, and Solar Noon
14+
15+
The module exports the following convenience functions:
16+
17+
* [`transit_sunrise_sunset`](@ref) — solar noon, sunrise, and sunset for a given day
18+
* [`next_sunrise`](@ref), [`previous_sunrise`](@ref) — sunrise times
19+
* [`next_sunset`](@ref), [`previous_sunset`](@ref) — sunset times
20+
* [`next_solar_noon`](@ref), [`previous_solar_noon`](@ref) — solar noon times
21+
22+
!!! info
23+
Solar noon is defined as the time when the sun reaches its highest elevation for the
24+
day. It is also referred to as the solar transit, the point the sun crosses the
25+
local meridian. The next/previous sunrise and sunset functions allow you to find the
26+
next or previous occurrence of these events relative to a given input time.
27+
28+
## Example: Sunrise and Sunset
29+
30+
As an example, let's calculate the solar noon, sunrise, and sunset times for the Van
31+
Gogh museum in Amsterdam on June 21, 2023:
32+
33+
```@example utilities
34+
using SolarPosition, Dates, TimeZones
35+
36+
tz_amsterdam = TimeZone("Europe/Amsterdam", TimeZones.Class(:LEGACY))
37+
obs = Observer(52.35888, 4.88185, 100.0)
38+
39+
# Summer solstice
40+
zdt = ZonedDateTime(2023, 6, 21, 12, 0, tz_amsterdam)
41+
42+
events = transit_sunrise_sunset(obs, zdt)
43+
44+
println("Solar noon: ", events.transit)
45+
println("Sunrise: ", events.sunrise)
46+
println("Sunset: ", events.sunset)
47+
```
48+
49+
We can confirm these results by consulting [timeanddate.com](https://www.timeanddate.com/sun/@52.35888,4.88185?month=6&year=2023) for our location and date.
50+
51+
Another option is to use the [`next_sunrise`](@ref) and [`next_sunset`](@ref) functions
52+
to return the sunrise and sunset times directly:
53+
54+
```@example utilities
55+
using SolarPosition, Dates, TimeZones
56+
57+
next_sunrise_time = next_sunrise(obs, zdt)
58+
next_sunset_time = next_sunset(obs, zdt)
59+
println("Next Sunrise: ", next_sunrise_time)
60+
println("Next Sunset: ", next_sunset_time)
61+
```
62+
63+
## Plotting the Solar Altitude
64+
65+
To visualize the solar altitude throughout the day, we can use the `solar_position`
66+
function to compute the solar positions at regular intervals and plot the results. We
67+
will make use of [`next_sunrise`](@ref) and [`next_sunset`](@ref) to mark the sunrise
68+
and sunset times on the plot.
69+
70+
!!! details "Visualization"
71+
```@example utilities
72+
using CairoMakie
73+
74+
# Define time range for the entire day (every 5 minutes)
75+
start_time = ZonedDateTime(2023, 6, 21, 0, 0, tz_amsterdam)
76+
end_time = ZonedDateTime(2023, 6, 21, 23, 59, tz_amsterdam)
77+
times = collect(start_time:Minute(5):end_time)
78+
79+
# Compute solar positions for all times
80+
positions = solar_position(obs, times)
81+
82+
# Get key events
83+
events = transit_sunrise_sunset(obs, zdt)
84+
sunrise_elev = solar_position(obs, events.sunrise).elevation
85+
sunset_elev = solar_position(obs, events.sunset).elevation
86+
transit_elev = solar_position(obs, events.transit).elevation
87+
88+
# Convert times to hours for plotting
89+
times_hours = hour.(times) .+ minute.(times) ./ 60
90+
91+
# Create the plot with styling
92+
fig = Figure(backgroundcolor=:transparent, textcolor="#f5ab35", size=(800, 400))
93+
ax = Axis(fig[1, 1],
94+
xlabel="Time (hours)",
95+
ylabel="Solar Altitude (°)",
96+
title="Solar Altitude - Amsterdam, June 21, 2023",
97+
backgroundcolor=:transparent,
98+
xticks=0:3:24)
99+
100+
# Plot the solar altitude curve
101+
lines!(ax, times_hours, positions.elevation,
102+
linewidth=2, color="#f5ab35")
103+
104+
# Add vertical markers and labels for events
105+
sunrise_hour = hour(events.sunrise) + minute(events.sunrise) / 60
106+
transit_hour = hour(events.transit) + minute(events.transit) / 60
107+
sunset_hour = hour(events.sunset) + minute(events.sunset) / 60
108+
109+
vlines!(ax, sunrise_hour, linestyle=:dash, color=:gold, linewidth=1.5)
110+
text!(ax, sunrise_hour, sunrise_elev + 5,
111+
text=Dates.format(events.sunrise, "HH:MM"),
112+
align=(:center, :bottom), color=:gold, fontsize=12)
113+
114+
vlines!(ax, transit_hour, linestyle=:dash, color=:red, linewidth=1.5)
115+
text!(ax, transit_hour, transit_elev + 5,
116+
text=Dates.format(events.transit, "HH:MM"),
117+
align=(:center, :bottom), color=:red, fontsize=12)
118+
119+
vlines!(ax, sunset_hour, linestyle=:dash, color=:purple, linewidth=1.5)
120+
text!(ax, sunset_hour, sunset_elev + 5,
121+
text=Dates.format(events.sunset, "HH:MM"),
122+
align=(:center, :bottom), color=:purple, fontsize=12)
123+
124+
# Add horizontal line at horizon
125+
hlines!(ax, 0, linestyle=:dash, color=:gray, linewidth=1)
126+
```
127+
128+
```@example utilities
129+
fig # hide
130+
```
131+
132+
As you can see, the sunrise and sunset events occur slightly below the horizon line (0°
133+
elevation). This is due to atmospheric [refraction](@ref refraction-correction) effects,
134+
which cause the sun to appear slightly higher in the sky when it is near the horizon.
135+
136+
## Sun Graph
137+
138+
We can plot the sunrise, sunset, and solar noon times on a sun graph to visualize the
139+
number of daylight hours throughout the day for our location in an entire year.
140+
141+
!!! details "Visualization"
142+
```@example utilities
143+
144+
# Generate dates for the entire year 2023 (every day)
145+
year_start = ZonedDateTime(2023, 1, 1, 12, 0, tz_amsterdam)
146+
year_dates = [year_start + Day(i) for i in 0:364]
147+
148+
# Calculate sunrise, sunset, and solar noon for each day
149+
sunrise_times = Float64[]
150+
sunset_times = Float64[]
151+
solar_noon_times = Float64[]
152+
153+
for date in year_dates
154+
events = transit_sunrise_sunset(obs, date)
155+
156+
# Convert to hours since midnight
157+
push!(sunrise_times, hour(events.sunrise) + minute(events.sunrise) / 60)
158+
push!(sunset_times, hour(events.sunset) + minute(events.sunset) / 60)
159+
push!(solar_noon_times, hour(events.transit) + minute(events.transit) / 60)
160+
end
161+
162+
# Calculate daylight hours for each day
163+
daylight_hours = sunset_times .- sunrise_times
164+
165+
# Find solstices (longest and shortest days)
166+
summer_solstice_idx = argmax(daylight_hours)
167+
winter_solstice_idx = argmin(daylight_hours)
168+
169+
# Create day of year array for x-axis
170+
day_of_year = 1:365
171+
172+
# Create the sun graph
173+
fig = Figure(backgroundcolor=:transparent, textcolor="#f5ab35", size=(900, 500))
174+
ax = Axis(fig[1, 1],
175+
xlabel="Month",
176+
ylabel="Time (hours)",
177+
title="Sun Graph 2023 - Amsterdam (52.36°N, 4.88°E)",
178+
backgroundcolor=:transparent,
179+
yticks=0:2:24,
180+
xticks=(
181+
[1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335],
182+
["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
183+
))
184+
185+
xlims!(ax, 1, 365)
186+
187+
# Fill night time (top part: midnight to sunrise)
188+
band!(ax, day_of_year, sunrise_times, fill(24.0, length(day_of_year)),
189+
color=(:darkslategray, 0.8))
190+
191+
# Fill night time (bottom part: sunset to midnight)
192+
band!(ax, day_of_year, fill(0.0, length(day_of_year)), sunset_times,
193+
color=(:darkslategray, 0.8))
194+
195+
# Fill daylight time
196+
band!(ax, day_of_year, sunrise_times, sunset_times,
197+
color=(:lightblue, 0.6))
198+
199+
# Plot solar noon line
200+
lines!(ax, day_of_year, solar_noon_times,
201+
color=:red, linewidth=2, label="Solar Noon")
202+
203+
# Mark solstices
204+
vlines!(ax, summer_solstice_idx, linestyle=:dash, color=:orange, linewidth=2)
205+
text!(ax, summer_solstice_idx, 12,
206+
text="Summer\nSolstice",
207+
align=(:center, :bottom), color=:orange, fontsize=10, rotation=π/2)
208+
209+
vlines!(ax, winter_solstice_idx, linestyle=:dash, color=:steelblue, linewidth=2)
210+
text!(ax, winter_solstice_idx, 12,
211+
text="Winter\nSolstice",
212+
align=(:center, :bottom), color=:steelblue, fontsize=10, rotation=π/2)
213+
214+
ylims!(ax, 0, 24)
215+
```
216+
217+
```@example utilities
218+
fig # hide
219+
```
220+
221+
We also marked the summer and winter solstices, which correspond to the longest and
222+
shortest days of the year, respectively.
223+
224+
!!! note
225+
Note the two discontinuities in March and October. These are due to the start and
226+
end of Daylight Saving Time (DST). The DST period starts on the last Sunday of March
227+
and ends on the last Sunday of October. Clocks are set one hour ahead in March,
228+
meaning sunrise and sunset times are later by one hour. This effect is reversed in
229+
October when clocks are set back one hour. This effectively turns the UTC offset
230+
from +1 hour to +2 hours during the DST period.
231+
232+
## Date, DateTime and ZonedDateTime
233+
234+
The utility functions accept three different time input types:
235+
236+
* **`ZonedDateTime`** — Recommended for timezone-aware calculations. The functions will
237+
return results in the same timezone as the input.
238+
* **`DateTime`** — Assumed to be in UTC. Results will be returned as `DateTime` in UTC.
239+
* **`Date`** — Assumed to be in UTC. Results will be returned as `DateTime` in UTC.
240+
241+
!!! tip
242+
When working with specific geographic locations, it's best to use `ZonedDateTime` to
243+
ensure results are in the local timezone and to correctly handle Daylight Saving Time
244+
transitions.
245+
246+
Here's an example showing the different input types:
247+
248+
```@example utilities
249+
using SolarPosition, Dates, TimeZones
250+
251+
obs = Observer(52.35888, 4.88185, 100.0) # Amsterdam
252+
tz_amsterdam = TimeZone("Europe/Amsterdam", TimeZones.Class(:LEGACY))
253+
254+
# Using ZonedDateTime (recommended - timezone aware)
255+
zdt = ZonedDateTime(2023, 6, 21, 12, 0, tz_amsterdam)
256+
events_zdt = transit_sunrise_sunset(obs, zdt)
257+
println("ZonedDateTime input:")
258+
println(" Sunrise: ", events_zdt.sunrise)
259+
260+
# Using DateTime (assumed UTC at 00:00)
261+
dt = DateTime(2023, 6, 21, 12, 0)
262+
events_dt = transit_sunrise_sunset(obs, dt)
263+
println("\nDateTime input (UTC):")
264+
println(" Sunrise: ", events_dt.sunrise)
265+
266+
# Using Date (assumed UTC at 00:00)
267+
d = Date(2023, 6, 21)
268+
events_d = transit_sunrise_sunset(obs, d)
269+
println("\nDate input (UTC 00:00):")
270+
println(" Sunrise: ", events_d.sunrise)
271+
```
272+
273+
Note that `DateTime` and `Date` inputs produce results in UTC, while `ZonedDateTime`
274+
preserves the input timezone. For Amsterdam in summer, the local time is UTC+2 (CEST),
275+
which explains the 2-hour difference in the sunrise times shown above.
276+
277+
## Forward looking functions
278+
279+
```@docs
280+
SolarPosition.Utilities.next_sunrise
281+
SolarPosition.Utilities.next_sunset
282+
SolarPosition.Utilities.next_solar_noon
283+
```
284+
285+
## Backward looking functions
286+
287+
```@docs
288+
SolarPosition.Utilities.previous_sunrise
289+
SolarPosition.Utilities.previous_sunset
290+
SolarPosition.Utilities.previous_solar_noon
291+
```
292+
293+
## Docs
294+
295+
```@docs
296+
SolarPosition.Utilities.transit_sunrise_sunset
297+
SolarPosition.Utilities.TransitSunriseSunset
298+
```

0 commit comments

Comments
 (0)