Merge pull request #16 from JuliaAstro/mtk-extension

Add ModelingToolkit.jl extension

Author: langestefan

Project.toml

@@ -12,19 +12,24 @@ TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [extensions]
 SolarPositionMakieExt = "Makie"
+SolarPositionModelingToolkitExt = ["ModelingToolkit", "Symbolics"]
 
 [compat]
 Aqua = "0.8"
 Dates = "1"
 DocStringExtensions = "0.8, 0.9"
 Makie = "0.24"
+ModelingToolkit = "10.3.0 - 10.26.1"
 StructArrays = "0.6, 0.7"
 Tables = "1"
 Test = "1"
 TimeZones = "1.22.0"
+Symbolics = "6,7"
 julia = "1.10"
 
 [extras]

docs/Project.toml

@@ -6,6 +6,8 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SolarPosition = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 

docs/make.jl

@@ -34,7 +34,11 @@ makedocs(;
     plugins = [bib],
     pages = [
         "index.md",
-        "Examples" => ["examples/getting-started.md", "examples/plotting.md"],
+        "Examples" => [
+            "examples/getting-started.md",
+            "examples/plotting.md",
+            "examples/modelingtoolkit.md",
+        ],
         "reference.md",
         "positioning.md",
         "refraction.md",

docs/src/examples/modelingtoolkit.md

@@ -0,0 +1,219 @@
+# Building models with ModelingToolkit.jl
+
+SolarPosition.jl provides a [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl)
+extension that enables integration of solar position calculations into symbolic modeling
+workflows. This allows you to compose solar position components with other physical
+systems for applications like solar energy modeling, building thermal analysis, and
+solar tracking systems.
+
+## Installation
+
+The ModelingToolkit extension is loaded automatically when both [`SolarPosition.jl`](https://github.com/JuliaAstro/SolarPosition.jl) and [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl)
+are loaded:
+
+```julia
+using SolarPosition
+using ModelingToolkit
+```
+
+## Quick Start
+
+The extension provides the [`SolarPositionBlock`](@ref) component, which outputs solar
+azimuth, elevation, and zenith angles as time-varying quantities.
+
+```@example mtk
+using SolarPosition
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t
+using Dates
+using OrdinaryDiffEq
+```
+
+```@example mtk
+# Create a solar position block
+@named sun = SolarPositionBlock()
+
+# Define observer location and reference time
+obs = Observer(51.50274937708521, -0.17782150375214803, 15.0)  # Natural History Museum
+t0 = DateTime(2024, 6, 21, 12, 0, 0)  # Summer solstice noon
+
+# Compile the system
+sys = mtkcompile(sun)
+
+# Set parameters using the compiled system's parameter references
+pmap = [
+    sys.observer => obs,
+    sys.t0 => t0,
+    sys.algorithm => PSA(),
+    sys.refraction => NoRefraction(),
+]
+
+# Solve over 24 hours (time in seconds)
+tspan = (0.0, 86400.0)
+prob = ODEProblem(sys, pmap, tspan)
+sol = solve(prob; saveat = 3600.0)  # Save every hour
+
+# Show some results
+println("Solar position at noon (t=12 hours):")
+println("  Azimuth: ", round(sol[sys.azimuth][1], digits=2), "°")
+println("  Elevation: ", round(sol[sys.elevation][1], digits=2), "°")
+println("  Zenith: ", round(sol[sys.zenith][1], digits=2), "°")
+```
+
+## SolarPositionBlock
+
+The [`SolarPositionBlock`](@ref) is a [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl) component that computes solar position angles based on time, observer location, and
+chosen positioning and refraction algorithms.
+
+```@docs
+SolarPositionBlock
+```
+
+## Composing with Other Systems
+
+The real power of the ModelingToolkit extension comes from composing solar position with other physical systems.
+
+### Example: Solar Panel Power Model
+
+```@example mtk
+using CairoMakie: Figure, Axis, lines!
+
+# Create solar position block
+@named sun = SolarPositionBlock()
+
+# Create a simple solar panel model
+@parameters begin
+    area = 10.0           # Panel area (m²)
+    efficiency = 0.2      # Panel efficiency (20%)
+    dni_peak = 1000.0     # Peak direct normal irradiance (W/m²)
+end
+
+@variables begin
+    irradiance(t) = 0.0   # Effective irradiance on panel (W/m²)
+    power(t) = 0.0        # Power output (W)
+end
+
+# Simplified model: irradiance depends on sun elevation
+# In reality, you'd account for panel orientation, azimuth, etc.
+eqs = [
+    irradiance ~ dni_peak * max(0, sind(sun.elevation)),
+    power ~ area * efficiency * irradiance,
+]
+
+# Compose the complete system
+@named model = System(eqs, t; systems = [sun])
+sys_model = mtkcompile(model)
+
+# Set up and solve
+obs = Observer(37.7749, -122.4194, 100.0)
+t0 = DateTime(2024, 6, 21, 0, 0, 0)
+
+pmap = [
+    sys_model.sun.observer => obs,
+    sys_model.sun.t0 => t0,
+    sys_model.sun.algorithm => PSA(),
+    sys_model.sun.refraction => NoRefraction(),
+]
+
+prob = ODEProblem(sys_model, pmap, (0.0, 86400.0))
+sol = solve(prob; saveat = 600.0)  # Save every 10 minutes
+
+# Plot results
+fig = Figure(size = (1000, 400))
+
+ax1 = Axis(fig[1, 1]; xlabel = "Time (hours)", ylabel = "Elevation (°)", title = "Solar Elevation")
+lines!(ax1, sol.t ./ 3600, sol[sys_model.sun.elevation])
+
+ax2 = Axis(fig[1, 2]; xlabel = "Time (hours)", ylabel = "Power (W)", title = "Solar Panel Power")
+lines!(ax2, sol.t ./ 3600, sol[sys_model.power])
+
+fig
+```
+
+### Example: Building Thermal Model with Solar Gain
+
+```@example mtk
+using CairoMakie: Figure, Axis, lines!
+using ModelingToolkit: D_nounits as D
+
+# Solar position component
+@named sun = SolarPositionBlock()
+
+# Building thermal model with solar gain
+@parameters begin
+    mass = 1000.0         # Thermal mass (kg)
+    cp = 1000.0           # Specific heat capacity (J/(kg·K))
+    U = 0.5               # Overall heat transfer coefficient (W/(m²·K))
+    wall_area = 50.0      # Wall area (m²)
+    window_area = 5.0     # Window area (m²)
+    window_trans = 0.7    # Window transmittance
+    T_outside = 20.0      # Outside temperature (°C)
+    dni_peak = 800.0      # Peak solar irradiance (W/m²)
+end
+
+@variables begin
+    T(t) = 20.0           # Room temperature (°C)
+    Q_loss(t)             # Heat loss through walls (W)
+    Q_solar(t)            # Solar heat gain (W)
+    irradiance(t)         # Solar irradiance (W/m²)
+end
+
+eqs = [
+    # Solar irradiance based on sun elevation
+    irradiance ~ dni_peak * max(0, sind(sun.elevation)),
+    # Solar heat gain through windows
+    Q_solar ~ window_area * window_trans * irradiance,
+    # Heat loss through walls
+    Q_loss ~ U * wall_area * (T - T_outside),
+    # Energy balance
+    D(T) ~ (Q_solar - Q_loss) / (mass * cp),
+]
+
+@named building = System(eqs, t; systems = [sun])
+sys_building = mtkcompile(building)
+
+# Simulate
+obs = Observer(40.7128, -74.0060, 100.0)  # New York City
+t0 = DateTime(2024, 6, 21, 0, 0, 0)
+
+pmap = [
+    sys_building.sun.observer => obs,
+    sys_building.sun.t0 => t0,
+    sys_building.sun.algorithm => PSA(),
+    sys_building.sun.refraction => NoRefraction(),
+]
+
+prob = ODEProblem(sys_building, pmap, (0.0, 86400.0))
+sol = solve(prob; saveat = 600.0)
+
+# Plot temperature evolution
+fig = Figure(size = (1200, 400))
+
+ax1 = Axis(fig[1, 1]; xlabel = "Time (hours)", ylabel = "Temperature (°C)", title = "Room Temperature")
+lines!(ax1, sol.t ./ 3600, sol[sys_building.T])
+
+ax2 = Axis(fig[1, 2]; xlabel = "Time (hours)", ylabel = "Solar Gain (W)", title = "Solar Heat Gain")
+lines!(ax2, sol.t ./ 3600, sol[sys_building.Q_solar])
+
+ax3 = Axis(fig[1, 3]; xlabel = "Time (hours)", ylabel = "Elevation (°)", title = "Sun Elevation")
+lines!(ax3, sol.t ./ 3600, sol[sys_building.sun.elevation])
+
+fig
+```
+
+## Implementation Details
+
+The extension works by registering the [`solar_position`](@ref) function and helper functions as
+symbolic operations in ModelingToolkit. The actual solar position calculation happens
+during ODE solving, with the simulation time `t` being converted to a [`DateTime`](https://docs.julialang.org/en/v1/stdlib/Dates/#Dates.DateTime) relative to the reference time `t0`.
+
+## Limitations
+
+The solar position calculation is treated as a black-box function by MTK's symbolic
+engine, so its internals cannot be symbolically simplified.
+
+## See Also
+
+- [Solar Positioning](@ref solar-positioning-algorithms) - Available positioning algorithms
+- [Refraction Correction](@ref refraction-correction) - Atmospheric refraction methods
+- [ModelingToolkit.jl Documentation](https://docs.sciml.ai/ModelingToolkit/stable/) - MTK framework documentation

docs/src/index.md

@@ -24,6 +24,14 @@ azimuth angles, which are essential for various applications where the sun is im
 - Climate studies
 - Astronomy
 
+## Extensions
+
+SolarPosition.jl provides package extensions for advanced use cases:
+
+- **ModelingToolkit Extension**: Integrate solar position calculations into symbolic modeling workflows. Create composable solar energy system models with ModelingToolkit.jl. See the [ModelingToolkit Extension](examples/modelingtoolkit.md) guide for details.
+
+- **Makie Extension**: Plotting recipes for solar position visualization.
+
 ## Acknowledgement
 
 This package is based on the work done by readers in the field of solar photovoltaics

examples/Project.toml

@@ -0,0 +1,6 @@
+[deps]
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SolarPosition = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"

ext/SolarPositionModelingToolkitExt.jl

@@ -0,0 +1,59 @@
+module SolarPositionModelingToolkitExt
+
+using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
+using SolarPosition: SolPos, ApparentSolPos, SPASolPos, AbstractApparentSolPos
+using ModelingToolkit: @parameters, @variables, System
+using ModelingToolkit: t_nounits as t
+using Dates: DateTime, Millisecond
+import Symbolics
+
+import SolarPosition: SolarPositionBlock, solar_position
+
+
+seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1e3))
+
+# helper functions to extract fields from solar position
+get_azimuth(pos) = pos.azimuth
+
+# for SolPos: use elevation and zenith
+get_elevation(pos::SolPos) = pos.elevation
+get_zenith(pos::SolPos) = pos.zenith
+
+# for ApparentSolPos and SPASolPos: use apparent_elevation and apparent_zenith
+get_elevation(pos::AbstractApparentSolPos) = pos.apparent_elevation
+get_zenith(pos::AbstractApparentSolPos) = pos.apparent_zenith
+
+Symbolics.@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
+Symbolics.@register_symbolic solar_position(
+    observer::Observer,
+    time::DateTime,
+    algorithm::SolarAlgorithm,
+    refraction::RefractionAlgorithm,
+)
+
+Symbolics.@register_symbolic get_azimuth(pos)
+Symbolics.@register_symbolic get_elevation(pos)
+Symbolics.@register_symbolic get_zenith(pos)
+
+function SolarPositionBlock(; name)
+    @parameters t0::DateTime [tunable = false] observer::Observer [tunable = false]
+    @parameters algorithm::SolarAlgorithm = PSA() [tunable = false]
+    @parameters refraction::RefractionAlgorithm = NoRefraction() [tunable = false]
+
+    @variables azimuth(t) [output = true]
+    @variables elevation(t) [output = true]
+    @variables zenith(t) [output = true]
+
+    time_expr = Symbolics.term(seconds_to_datetime, t, t0; type = DateTime)
+    pos = solar_position(observer, time_expr, algorithm, refraction)
+
+    eqs = [
+        azimuth ~ get_azimuth(pos),
+        elevation ~ get_elevation(pos),
+        zenith ~ get_zenith(pos),
+    ]
+
+    return System(eqs, t; name = name)
+end
+
+end # module SolarPositionModelingToolkitExt

src/Positioning/Positioning.jl

@@ -74,10 +74,10 @@ struct Observer{T<:AbstractFloat}
         # apply pole corrections to avoid numerical issues
         if lat == 90.0
             lat -= 1e-6
-            @warn "Latitude was 90°. Adjusted to $(lat)° to avoid singularities."
+            @warn "Latitude is 90°. Adjusted to $(lat)° to avoid singularities." maxlog = 1
         elseif lat == -90.0
             lat += 1e-6
-            @warn "Latitude was -90°. Adjusted to $(lat)° to avoid singularities."
+            @warn "Latitude is -90°. Adjusted to $(lat)° to avoid singularities." maxlog = 1
         end
 
         lat_rad = deg2rad(lat)
@@ -97,6 +97,7 @@ Base.show(io::IO, obs::Observer) = print(
 )
 
 abstract type AbstractSolPos end
+abstract type AbstractApparentSolPos <: AbstractSolPos end
 
 """
     $(TYPEDEF)
@@ -126,7 +127,7 @@ Also includes apparent elevation and zenith angles.
 # Fields
 $(TYPEDFIELDS)
 """
-struct ApparentSolPos{T} <: AbstractSolPos where {T<:AbstractFloat}
+struct ApparentSolPos{T} <: AbstractApparentSolPos where {T<:AbstractFloat}
     "Azimuth (degrees, 0=N, +clockwise, range [-180, 180])"
     azimuth::T
     "Elevation (degrees, range [-90, 90])"
@@ -159,7 +160,7 @@ Solar position result from SPA algorithm including equation of time.
 # Fields
 $(TYPEDFIELDS)
 """
-struct SPASolPos{T} <: AbstractSolPos where {T<:AbstractFloat}
+struct SPASolPos{T} <: AbstractApparentSolPos where {T<:AbstractFloat}
     "Azimuth (degrees, 0=N, +clockwise, range [-180, 180])"
     azimuth::T
     "Elevation (degrees, range [-90, 90])"

src/Positioning/deltat.jl

@@ -106,7 +106,7 @@ The polynomial expressions for ΔT are from [NASADeltaT](@cite), based on the wo
 """
 function calculate_deltat(year::Real, month::Real)
     if year < -1999 || year > 3000
-        @warn "ΔT is undefined for years before -1999 or after 3000."
+        @warn "ΔT is undefined for years before -1999 or after 3000." maxlog = 1
     end
 
     y = year + (month - 0.5) / 12