From 8848b063f7737fa58d96f6b17b45fbca5fa20c77 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Tue, 11 Nov 2025 23:23:48 +0100
Subject: [PATCH 01/18] Initial version of MTK extension

---
 Project.toml                           |   5 +
 examples/Project.toml                  |   6 +
 examples/mtk-integration.jl            | 314 +++++++++++++++++++++++++
 ext/SolarPositionModelingToolkitExt.jl |  52 ++++
 src/Positioning/Positioning.jl         |   5 +-
 src/SolarPosition.jl                   |  19 +-
 test/Project.toml                      |   4 +
 test/runtests.jl                       |  13 +-
 test/test-mtk.jl                       | 164 +++++++++++++
 9 files changed, 572 insertions(+), 10 deletions(-)
 create mode 100644 examples/Project.toml
 create mode 100644 examples/mtk-integration.jl
 create mode 100644 ext/SolarPositionModelingToolkitExt.jl
 create mode 100644 test/test-mtk.jl

diff --git a/Project.toml b/Project.toml
index 498055b..4deb3e9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,15 +12,20 @@ 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"
+Symbolics = "6, 7"
 StructArrays = "0.6, 0.7"
 Tables = "1"
 Test = "1"
diff --git a/examples/Project.toml b/examples/Project.toml
new file mode 100644
index 0000000..8123984
--- /dev/null
+++ b/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"
diff --git a/examples/mtk-integration.jl b/examples/mtk-integration.jl
new file mode 100644
index 0000000..1744639
--- /dev/null
+++ b/examples/mtk-integration.jl
@@ -0,0 +1,314 @@
+# # # ModelingToolkit Integration Example
+# #
+# # This example demonstrates how to use the SolarPosition.jl ModelingToolkit extension
+# # to integrate solar position calculations into larger modeling workflows.
+# #
+# # The ModelingToolkit extension allows you to:
+# # 1. Create reusable solar position components
+# # 2. Compose solar models with other physical systems
+# # 3. Leverage MTK's symbolic simplification and code generation
+# # 4. Build complex solar energy system models
+
+# using ModelingToolkit
+# using ModelingToolkit: t_nounits as t, D_nounits as D
+# using SolarPosition
+# using Dates
+# using OrdinaryDiffEq
+# using Plots
+
+# # ## Example 1: Basic Solar Component
+# #
+# # Create a solar position component for a specific location
+
+# @named sun_sf = SolarPositionBlock(
+#     latitude = 37.7749,      # San Francisco
+#     longitude = -122.4194,
+#     altitude = 100.0,
+#     t0 = DateTime(2024, 6, 21, 0, 0, 0),  # Summer solstice midnight
+#     time_unit = Hour(1),      # t=1.0 means 1 hour from t0
+# )
+
+# # The solar component has variables for azimuth, elevation, and zenith angles
+# println("Solar component unknowns: ", unknowns(sun_sf))
+# println("Solar component parameters: ", parameters(sun_sf))
+# println("Solar component equations: ", equations(sun_sf))
+
+# # ## Example 2: Solar Panel Power Model
+# #
+# # Build a simple solar panel model that uses the solar position
+
+# function simple_solar_panel(; name)
+#     @parameters begin
+#         area = 10.0           # Panel area in m²
+#         efficiency = 0.20     # Panel efficiency (20%)
+#         tilt_angle = 30.0     # Panel tilt from horizontal (degrees)
+#         azimuth_angle = 180.0 # Panel azimuth (180° = South facing)
+#     end
+
+#     @variables begin
+#         power(t)              # Power output in W
+#         dni(t) = 800.0        # Direct Normal Irradiance in W/m²
+#         incidence_angle(t)    # Angle between sun and panel normal
+#         effective_irr(t)      # Effective irradiance on panel
+#     end
+
+#     # Note: This is a simplified model for demonstration
+#     # In a real model, you would connect sun.elevation and sun.azimuth
+#     # to calculate the actual incidence angle
+#     eqs = [
+#         # Simplified: assume some base irradiance pattern
+#         dni ~ 800.0 + 200.0 * sin(t * π / 12),  # Daily variation
+#         # Effective irradiance (simplified)
+#         effective_irr ~ dni * max(0.0, cos(deg2rad(45.0))),
+#         # Power output
+#         power ~ area * efficiency * effective_irr,
+#     ]
+
+#     return System(eqs, t; name = name)
+# end
+
+# @named panel = simple_solar_panel()
+
+# # Compose the solar component with the panel model
+# @named solar_system = compose(System(Equation[], t; name = :solar_system), sun_sf, panel)
+
+# println("\nComposed system unknowns: ", length(unknowns(solar_system)))
+# println("Composed system equations: ", length(equations(solar_system)))
+
+# # ## Example 3: Working Simulation with Solar Panel
+# #
+# # Create a complete working example that can be compiled and simulated
+
+# function solar_panel_with_sun(; name, sun_elevation)
+#     @parameters begin
+#         area = 10.0           # Panel area in m²
+#         efficiency = 0.20     # Panel efficiency (20%)
+#         base_irradiance = 1000.0  # Base solar irradiance W/m²
+#     end
+
+#     @variables begin
+#         power(t)              # Power output in W
+#         irradiance(t)         # Effective irradiance on panel
+#     end
+
+#     eqs = [
+#         # Irradiance depends on sun elevation
+#         irradiance ~ base_irradiance * max(0.0, sin(deg2rad(sun_elevation))),
+#         # Power output
+#         power ~ area * efficiency * irradiance,
+#     ]
+
+#     return System(eqs, t; name = name)
+# end
+
+# # Pre-compute solar position for a specific time
+# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
+# ref_time = DateTime(2024, 6, 21, 20, 0, 0)  # Summer solstice, local solar noon (UTC)
+# ref_pos = solar_position(obs, ref_time, PSA())
+
+# println("\nReference solar position at summer solstice noon:")
+# println("  Elevation: $(ref_pos.elevation)°")
+# println("  Azimuth: $(ref_pos.azimuth)°")
+
+# # Create solar panel with the computed elevation
+# @named panel_working = solar_panel_with_sun(sun_elevation = ref_pos.elevation)
+
+# # Compile and simulate
+# panel_sys = mtkcompile(panel_working)
+# prob = ODEProblem(panel_sys, Dict(), (0.0, 10.0))
+# sol = solve(prob, Tsit5())
+
+# println("\nSolar panel simulation:")
+# println("  Power output: $(sol[panel_working.power][end]) W")
+# println("  Irradiance: $(sol[panel_working.irradiance][end]) W/m²")
+
+# # ## Example 4: Building Thermal Model with Solar Gain
+# #
+# # Model a simple building with solar heat gain through windows
+
+# function building_room(; name)
+#     @parameters begin
+#         mass = 1000.0           # Thermal mass in kg
+#         cp = 1000.0             # Specific heat capacity J/(kg·K)
+#         U_wall = 0.5            # Wall U-value W/(m²·K)
+#         wall_area = 50.0        # Wall area m²
+#         T_outside = 20.0        # Outside temperature °C
+#     end
+
+#     @variables begin
+#         T(t) = 20.0             # Room temperature °C
+#         Q_loss(t)               # Heat loss through walls W
+#         Q_solar(t) = 0.0        # Solar heat gain W (placeholder)
+#     end
+
+#     eqs = [
+#         # Heat loss through walls
+#         Q_loss ~ U_wall * wall_area * (T - T_outside),
+#         # Energy balance (simplified)
+#         D(T) ~ (Q_solar - Q_loss) / (mass * cp),
+#     ]
+
+#     return System(eqs, t; name = name)
+# end
+
+# function solar_window(; name)
+#     @parameters begin
+#         window_area = 5.0       # Window area m²
+#         transmittance = 0.7     # Glass transmittance
+#         normal_azimuth = 180.0  # South-facing
+#     end
+
+#     @variables begin
+#         Q_solar(t)              # Solar heat gain W
+#         irradiance(t)           # Solar irradiance W/m²
+#     end
+
+#     # Simplified solar gain calculation
+#     eqs = [
+#         # Base irradiance pattern (in real model, use sun.elevation)
+#         irradiance ~ 500.0 * (1 + sin(t * π / 12)),
+#         Q_solar ~ window_area * transmittance * irradiance,
+#     ]
+
+#     return System(eqs, t; name = name)
+# end
+
+# @named room = building_room()
+# @named window = solar_window()
+# @named sun_building = SolarPositionBlock(latitude = 40.7128, longitude = -74.0060)  # NYC
+
+# # Connect the window solar gain to the room
+# @named building = compose(
+#     System([room.Q_solar ~ window.Q_solar], t; name = :building),
+#     room,
+#     window,
+#     sun_building,
+# )
+
+# println("\nBuilding model unknowns: ", length(unknowns(building)))
+# println("Building model parameters: ", length(parameters(building)))
+
+# # ## Example 5: Time-Varying Solar Simulation
+# #
+# # Simulate a simple system over a day with changing solar position
+
+# function simple_thermal_mass(; name, Q_solar_func)
+#     @parameters begin
+#         mass = 100.0        # Thermal mass in kg
+#         cp = 1000.0         # Specific heat J/(kg·K)
+#         T_amb = 20.0        # Ambient temperature °C
+#         h = 5.0             # Heat transfer coefficient W/(m²·K)
+#         area = 1.0          # Surface area m²
+#     end
+
+#     @variables begin
+#         T(t) = 20.0         # Temperature °C
+#         Q_solar(t)          # Solar heat input W
+#     end
+
+#     eqs = [Q_solar ~ Q_solar_func, D(T) ~ (Q_solar - h * area * (T - T_amb)) / (mass * cp)]
+
+#     return System(eqs, t; name = name)
+# end
+
+# # Create a time-varying solar input (simplified as sinusoidal)
+# # In practice, you'd use pre-computed solar positions with interpolation
+# solar_input_pattern = 500.0 * max(0.0, sin(π * t / 12))  # Peaks at noon (t=6)
+
+# @named thermal = simple_thermal_mass(Q_solar_func = solar_input_pattern)
+# thermal_sys = mtkcompile(thermal)
+
+# # Simulate over 24 hours
+# prob_thermal = ODEProblem(thermal_sys, Dict(), (0.0, 24.0))
+# sol_thermal = solve(prob_thermal, Tsit5())
+
+# println("\nThermal mass simulation with solar input:")
+# println("  Initial temperature: $(sol_thermal[thermal.T][1]) °C")
+# println("  Final temperature: $(sol_thermal[thermal.T][end]) °C")
+# println("  Peak temperature: $(maximum(sol_thermal[thermal.T])) °C")
+
+# # Plot temperature evolution
+# p_thermal = plot(
+#     sol_thermal,
+#     idxs = [thermal.T],
+#     xlabel = "Time (hours)",
+#     ylabel = "Temperature (°C)",
+#     label = "Temperature",
+#     title = "Thermal Mass with Solar Heating",
+#     linewidth = 2,
+#     legend = :best,
+# )
+
+# display(p_thermal)
+
+# # ## Example 6: Using Real Solar Position Data
+# #
+# # For accurate solar position calculations, you can compute positions offline
+# # and use them as time-dependent parameters or callbacks
+
+# # Create observer and time range
+# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
+# start_time = DateTime(2024, 6, 21, 0, 0, 0)  # Summer solstice
+# times = [start_time + Hour(h) for h = 0:23]
+
+# # Compute solar positions
+# positions = solar_position(obs, times, PSA())
+
+# # Extract data
+# hours = 0:23
+# azimuths = positions.azimuth
+# elevations = positions.elevation
+# zeniths = positions.zenith
+
+# # Plot solar position over the day
+# p3 = plot(
+#     hours,
+#     elevations,
+#     xlabel = "Hour of Day",
+#     ylabel = "Solar Elevation (°)",
+#     label = "Elevation Angle",
+#     title = "Solar Position - Summer Solstice, San Francisco",
+#     linewidth = 2,
+#     marker = :circle,
+#     legend = :best,
+# )
+
+# p4 = plot(
+#     hours,
+#     azimuths,
+#     xlabel = "Hour of Day",
+#     ylabel = "Solar Azimuth (°)",
+#     label = "Azimuth Angle",
+#     title = "Solar Azimuth Over the Day",
+#     linewidth = 2,
+#     marker = :circle,
+#     legend = :best,
+# )
+
+# plot(p3, p4, layout = (2, 1), size = (800, 600))
+
+# # ## Notes for Advanced Usage
+# #
+# # 1. **Time Mapping**: In a real application, you need to map simulation time `t`
+# #    to actual DateTime values. This can be done with callbacks or custom functions.
+# #
+# # 2. **Callbacks**: For dynamic solar position updates during simulation, implement
+# #    a DiscreteCallback that updates solar position parameters at each time step.
+# #
+# # 3. **Performance**: For long simulations, consider pre-computing solar positions
+# #    and using interpolation rather than calculating at each solver step.
+# #
+# # 4. **Refraction**: Include atmospheric refraction for more accurate apparent
+# #    solar positions, especially important near sunrise/sunset.
+# #
+# # 5. **Integration**: The MTK extension enables integration with:
+# #    - PV system models
+# #    - Solar thermal collectors
+# #    - Building energy models (EnergyPlus-like)
+# #    - Daylighting calculations
+# #    - Solar tracking systems
+# #    - Agricultural/greenhouse models
+
+# println("\n✓ ModelingToolkit integration examples completed!")
+# println("  The extension enables composable solar energy system modeling.")
+# println("  See the documentation for more advanced usage patterns.")
diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
new file mode 100644
index 0000000..7e6021b
--- /dev/null
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -0,0 +1,52 @@
+module SolarPositionModelingToolkitExt
+
+using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
+using ModelingToolkit: @parameters, @variables, System, @register_symbolic, t_nounits
+using Symbolics
+using Symbolics: term
+using Dates: DateTime, Millisecond
+
+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
+get_elevation(pos) = pos.elevation
+get_zenith(pos) = pos.zenith
+
+@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
+@register_symbolic solar_position(
+    observer::Observer,
+    time::DateTime,
+    algorithm::SolarAlgorithm,
+    refraction::RefractionAlgorithm,
+)
+@register_symbolic get_azimuth(pos)
+@register_symbolic get_elevation(pos)
+@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_nounits) [output = true]
+    @variables elevation(t_nounits) [output = true]
+    @variables zenith(t_nounits) [output = true]
+
+    time_expr = term(seconds_to_datetime, t_nounits, 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_nounits; name = name)
+end
+
+end # module SolarPositionModelingToolkitExt
diff --git a/src/Positioning/Positioning.jl b/src/Positioning/Positioning.jl
index a789892..48b2225 100644
--- a/src/Positioning/Positioning.jl
+++ b/src/Positioning/Positioning.jl
@@ -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])"
diff --git a/src/SolarPosition.jl b/src/SolarPosition.jl
index 20ca2cd..3483cef 100644
--- a/src/SolarPosition.jl
+++ b/src/SolarPosition.jl
@@ -5,14 +5,24 @@ include("Positioning/Positioning.jl")
 
 using .Positioning:
     Observer, PSA, NOAA, Walraven, USNO, SPA, solar_position, solar_position!
-using .Positioning: SolPos, ApparentSolPos, SPASolPos
+using .Positioning:
+    SolPos,
+    ApparentSolPos,
+    SPASolPos,
+    SolarAlgorithm,
+    AbstractApparentSolPos,
+    AbstractSolPos
 using .Refraction: RefractionAlgorithm, NoRefraction
 using .Refraction: HUGHES, ARCHER, BENNETT, MICHALSKY, SG2, SPARefraction
 
-export solar_position, solar_position!, Observer, PSA, NOAA, Walraven, USNO, SPA
+export solar_position, solar_position!, SolarAlgorithm, Observer
+export PSA, NOAA, Walraven, USNO, SPA
+
 export RefractionAlgorithm, NoRefraction
 export HUGHES, ARCHER, BENNETT, MICHALSKY, SG2, SPARefraction
+
 export SolPos, ApparentSolPos, SPASolPos
+export AbstractSolPos, AbstractApparentSolPos
 
 # to make the makie extension work
 export sunpathplot
@@ -25,4 +35,9 @@ function sunpathplot! end
 function sunpathpolarplot end
 function sunpathpolarplot! end
 
+# to make the ModelingToolkit extension work
+export SolarPositionBlock
+
+function SolarPositionBlock end
+
 end # module
diff --git a/test/Project.toml b/test/Project.toml
index 736cacb..4d886c4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,10 +1,14 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SolarPosition = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
diff --git a/test/runtests.jl b/test/runtests.jl
index 40b4e54..6f93d83 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,10 +1,11 @@
 using Test
-using SolarPosition
-using SolarPosition.Positioning:
-    Observer, NOAA, PSA, Walraven, USNO, solar_position, NoRefraction
-using SolarPosition.Refraction: HUGHES
-using Dates, TimeZones
-using DataFrames
+# using SolarPosition
+# using SolarPosition.Positioning:
+#     Observer, NOAA, PSA, Walraven, USNO, solar_position, NoRefraction
+# using SolarPosition.Refraction: HUGHES
+# using Dates, TimeZones
+# using DataFrames
+# using ModelingToolkit
 
 include("setup.jl")
 include("linting.jl")
diff --git a/test/test-mtk.jl b/test/test-mtk.jl
new file mode 100644
index 0000000..f0b7f56
--- /dev/null
+++ b/test/test-mtk.jl
@@ -0,0 +1,164 @@
+using Test
+using SolarPosition:
+    Observer, solar_position, PSA, NoRefraction, SolarAlgorithm, SolarPositionBlock, BENNETT
+using ModelingToolkit:
+    @named, @variables, @parameters, get_variables, unknowns, System, mtkcompile, observed
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using Dates: DateTime
+using OrdinaryDiffEq
+using CairoMakie
+
+@testset "ModelingToolkit Extension" begin
+    obs = Observer(37.7749, -122.4194, 100.0)
+    t0 = DateTime(2024, 6, 21, 0, 0, 0)  # Summer solstice
+
+    @testset "SolarPositionBlock Creation No Refraction" begin
+        @named sun = SolarPositionBlock()
+        @test sun isa System
+        @test sun.azimuth isa Any
+        @test sun.elevation isa Any
+        @test sun.zenith isa Any
+    end
+
+    @testset "Component Composition" begin
+        @named sun = SolarPositionBlock()
+
+        @parameters begin
+            area = 10.0
+            efficiency = 0.2
+        end
+
+        @variables begin
+            irradiance(t) = 0.0
+            power(t) = 0.0
+        end
+
+        eqs = [
+            irradiance ~ 1000 * max(0, sind(sun.elevation) * cosd(sun.azimuth)),
+            power ~ area * efficiency * irradiance,
+        ]
+
+        @named model = System(eqs, t; systems = [sun])
+        sys = mtkcompile(model)
+
+        @test sys isa System
+        # test if we can access variables without errors
+        @test sys.irradiance isa Any
+        @test sys.power isa Any
+    end
+
+    @testset "ODEProblem Solution and Plotting" begin
+        @named sun = SolarPositionBlock()
+
+        # Compile the system
+        sys = mtkcompile(sun)
+
+        # Create parameter dictionary using the parameters from the compiled system
+        pmap = [
+            sys.observer => obs,
+            sys.t0 => t0,
+            sys.algorithm => PSA(),
+            sys.refraction => NoRefraction(),
+        ]
+
+        # Create ODEProblem for 24 hours (86400 seconds)
+        tspan = (0.0, 86400.0)
+        prob = ODEProblem(sys, pmap, tspan)
+
+        # Solve with fixed time step - save every 60 seconds
+        sol = solve(prob; saveat = 60.0)
+
+        @test sol isa Any
+        @test length(sol.t) > 0
+        @test length(sol.t) == 1441  # 0 to 86400 in steps of 60 = 1441 points
+
+        # Validate that the solution contains actual values
+        azimuth_vals = sol[sys.azimuth]
+        elevation_vals = sol[sys.elevation]
+        zenith_vals = sol[sys.zenith]
+
+        @test length(azimuth_vals) > 0
+        @test length(elevation_vals) > 0
+        @test length(zenith_vals) > 0
+
+        # Verify that we have values at different times
+        @test length(azimuth_vals) == length(sol.t)
+        @test length(elevation_vals) == length(sol.t)
+        @test length(zenith_vals) == length(sol.t)
+
+        # Check that zenith ≈ 90 - elevation
+        @test all(isapprox.(zenith_vals, 90 .- elevation_vals; atol = 1e-6))
+
+        # For summer solstice in San Francisco, check that:
+        # 1. The sun rises (elevation goes from negative to positive)
+        # 2. The sun sets (elevation goes from positive to negative)
+        # 3. Maximum elevation is reasonable (should be around 70-75° for San Francisco summer solstice)
+        max_elevation = maximum(elevation_vals)
+        min_elevation = minimum(elevation_vals)
+
+        @test max_elevation > 60.0  # Sun should get high in sky on summer solstice
+        @test min_elevation < 0.0   # Sun should be below horizon at some point
+        @test max_elevation - min_elevation > 100.0  # Should have significant variation over 24h
+
+        # Check azimuth sweeps across a wide range
+        azimuth_range = maximum(azimuth_vals) - minimum(azimuth_vals)
+        @test azimuth_range > 180.0  # Sun should move significantly across the sky
+
+        println("Solution validation:")
+        println("  Time range: ", sol.t[1], " to ", sol.t[end], " seconds")
+        println("  Number of points: ", length(sol.t))
+        println("  Elevation range: ", min_elevation, "° to ", max_elevation, "°")
+        println(
+            "  Azimuth range: ",
+            minimum(azimuth_vals),
+            "° to ",
+            maximum(azimuth_vals),
+            "°",
+        )
+
+        # Create plots
+        fig = Figure(; size = (1200, 800))
+
+        # Plot azimuth
+        ax1 = Axis(
+            fig[1, 1];
+            xlabel = "Time (hours)",
+            ylabel = "Azimuth (°)",
+            title = "Solar Azimuth",
+        )
+        lines!(ax1, sol.t ./ 3600, sol[sys.azimuth])
+
+        # Plot elevation
+        ax2 = Axis(
+            fig[1, 2];
+            xlabel = "Time (hours)",
+            ylabel = "Elevation (°)",
+            title = "Solar Elevation",
+        )
+        lines!(ax2, sol.t ./ 3600, sol[sys.elevation])
+
+        # Plot zenith
+        ax3 = Axis(
+            fig[2, 1];
+            xlabel = "Time (hours)",
+            ylabel = "Zenith (°)",
+            title = "Solar Zenith",
+        )
+        lines!(ax3, sol.t ./ 3600, sol[sys.zenith])
+
+        # Sky plot (polar plot of azimuth vs elevation)
+        ax4 = Axis(
+            fig[2, 2];
+            xlabel = "Azimuth (°)",
+            ylabel = "Elevation (°)",
+            title = "Sun Path",
+        )
+        lines!(ax4, sol[sys.azimuth], sol[sys.elevation])
+
+        # Test that the figure was created
+        @test fig isa Figure
+
+        # Optionally save the plot (commented out to avoid file creation in tests)
+        save("test_solar_position_plot.png", fig)
+    end
+end

From d2e25e412b8f89b59a66c29e8912b51c7de4e7f9 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Wed, 12 Nov 2025 22:03:39 +0100
Subject: [PATCH 02/18] Update imports in MTK extension

---
 Project.toml                           |  4 ++--
 ext/SolarPositionModelingToolkitExt.jl | 20 +++++++++-----------
 test/Project.toml                      |  1 -
 3 files changed, 11 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4deb3e9..470ee99 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,12 +24,12 @@ Aqua = "0.8"
 Dates = "1"
 DocStringExtensions = "0.8, 0.9"
 Makie = "0.24"
-ModelingToolkit = "10"
-Symbolics = "6, 7"
+ModelingToolkit = "10.26.1"
 StructArrays = "0.6, 0.7"
 Tables = "1"
 Test = "1"
 TimeZones = "1.22.0"
+Symbolics = "6,7"
 julia = "1.10"
 
 [extras]
diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
index 7e6021b..96dbacb 100644
--- a/ext/SolarPositionModelingToolkitExt.jl
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -1,34 +1,32 @@
 module SolarPositionModelingToolkitExt
 
 using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
-using ModelingToolkit: @parameters, @variables, System, @register_symbolic, t_nounits
-using Symbolics
-using Symbolics: term
+using ModelingToolkit: @parameters, @variables, System, t_nounits
 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
+# helper functions to extract fields from solar position
 get_azimuth(pos) = pos.azimuth
 get_elevation(pos) = pos.elevation
 get_zenith(pos) = pos.zenith
 
-@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
-@register_symbolic solar_position(
+Symbolics.@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
+Symbolics.@register_symbolic solar_position(
     observer::Observer,
     time::DateTime,
     algorithm::SolarAlgorithm,
     refraction::RefractionAlgorithm,
 )
-@register_symbolic get_azimuth(pos)
-@register_symbolic get_elevation(pos)
-@register_symbolic get_zenith(pos)
+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]
@@ -37,7 +35,7 @@ function SolarPositionBlock(; name)
     @variables elevation(t_nounits) [output = true]
     @variables zenith(t_nounits) [output = true]
 
-    time_expr = term(seconds_to_datetime, t_nounits, t0; type = DateTime)
+    time_expr = Symbolics.term(seconds_to_datetime, t_nounits, t0; type = DateTime)
     pos = solar_position(observer, time_expr, algorithm, refraction)
 
     eqs = [
diff --git a/test/Project.toml b/test/Project.toml
index 4d886c4..2db91c1 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,7 +8,6 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SolarPosition = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"

From 8ac601e88fc67e872b06a1253054a697160708fa Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Thu, 20 Nov 2025 03:40:55 +0100
Subject: [PATCH 03/18] Add docs for ModelingToolkit.jl extension

---
 Project.toml                           |   2 +-
 docs/Project.toml                      |   2 +
 docs/make.jl                           |   6 +-
 docs/src/examples/modelingtoolkit.md   | 218 +++++++++++++++++++++++++
 docs/src/index.md                      |   8 +
 ext/SolarPositionModelingToolkitExt.jl |   1 +
 src/SolarPosition.jl                   |  61 +++++++
 7 files changed, 296 insertions(+), 2 deletions(-)
 create mode 100644 docs/src/examples/modelingtoolkit.md

diff --git a/Project.toml b/Project.toml
index 470ee99..b6ba11d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,7 +24,7 @@ Aqua = "0.8"
 Dates = "1"
 DocStringExtensions = "0.8, 0.9"
 Makie = "0.24"
-ModelingToolkit = "10.26.1"
+ModelingToolkit = "10.3.0 - 10.26.1"
 StructArrays = "0.6, 0.7"
 Tables = "1"
 Test = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 2c62fa5..ce1956c 100644
--- a/docs/Project.toml
+++ b/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"
 
diff --git a/docs/make.jl b/docs/make.jl
index d88a76a..c69b19f 100644
--- a/docs/make.jl
+++ b/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",
diff --git a/docs/src/examples/modelingtoolkit.md b/docs/src/examples/modelingtoolkit.md
new file mode 100644
index 0000000..842f513
--- /dev/null
+++ b/docs/src/examples/modelingtoolkit.md
@@ -0,0 +1,218 @@
+# 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
+
+# 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, Tsit5(); 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
diff --git a/docs/src/index.md b/docs/src/index.md
index 535b580..c66dccd 100644
--- a/docs/src/index.md
+++ b/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
diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
index 96dbacb..4a9153f 100644
--- a/ext/SolarPositionModelingToolkitExt.jl
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -26,6 +26,7 @@ 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]
diff --git a/src/SolarPosition.jl b/src/SolarPosition.jl
index 3483cef..6ec5d55 100644
--- a/src/SolarPosition.jl
+++ b/src/SolarPosition.jl
@@ -3,6 +3,8 @@ module SolarPosition
 include("Refraction/Refraction.jl")
 include("Positioning/Positioning.jl")
 
+using DocStringExtensions: TYPEDSIGNATURES
+
 using .Positioning:
     Observer, PSA, NOAA, Walraven, USNO, SPA, solar_position, solar_position!
 using .Positioning:
@@ -38,6 +40,65 @@ function sunpathpolarplot! end
 # to make the ModelingToolkit extension work
 export SolarPositionBlock
 
+"""
+    $(TYPEDSIGNATURES)
+
+Return a [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl) component
+that computes solar position as a function of time and can be integrated into symbolic
+modeling workflows.
+
+The [`SolarPositionBlock`](@ref) is a [`System`](https://docs.sciml.ai/ModelingToolkit/stable/API/System/#ModelingToolkit.System)
+which exposes `azimuth`, `elevation`, and `zenith` as output variables computed from the
+simulation time `t` (in seconds) relative to a reference time `t0`.
+
+!!! note
+    This function requires [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl)
+    to be loaded. The extension is automatically loaded when both [SolarPosition.jl](https://github.com/JuliaAstro/SolarPosition.jl)
+    and [`ModelingToolkit.jl`](https://github.com/SciML/ModelingToolkit.jl) are available.
+
+# Parameters
+
+- `observer::Observer`: location (latitude, longitude, altitude). See [`Observer`](@ref)
+- `t0::DateTime`: Reference time (time when simulation time `t = 0`)
+- `algorithm::SolarAlgorithm`: Solar positioning algorithm (default: [`PSA`](@ref))
+- `refraction::RefractionAlgorithm`: Atmospheric refraction correction (default: [`NoRefraction`](@ref))
+
+# Variables (Outputs)
+
+- `azimuth(t)`: Solar azimuth angle in degrees (0° = North, 90° = East, 180° = South, 270° = West)
+- `elevation(t)`: Solar elevation angle in degrees (angle above horizon, positive when sun is visible)
+- `zenith(t)`: Solar zenith angle in degrees (angle from vertical, complementary to elevation: `zenith = 90° - elevation`)
+
+# Time Convention
+
+The simulation time `t` (accessed via `t_nounits`) is in **seconds** from the reference time `t0`. For example:
+- `t = 0` corresponds to `t0`
+- `t = 3600` corresponds to `t0 + 1 hour`
+- `t = 86400` corresponds to `t0 + 24 hours`
+
+# Example
+
+```julia
+using SolarPosition, ModelingToolkit
+using ModelingToolkit: t_nounits as t
+using Dates
+
+@named sun = SolarPositionBlock()
+obs = Observer(51.5, -0.18, 15.0)
+t0 = DateTime(2024, 6, 21, 12, 0, 0)
+
+sys = mtkcompile(sun)
+pmap = [
+    sys.observer => obs,
+    sys.t0 => t0,
+    sys.algorithm => PSA(),
+    sys.refraction => NoRefraction(),
+]
+
+prob = ODEProblem(sys, pmap, (0.0, 86400.0))
+sol = solve(prob; saveat = 3600.0)
+```
+"""
 function SolarPositionBlock end
 
 end # module

From 38110dc89d8c7327cc884c58d2e493776bbfa312 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <37669773+langestefan@users.noreply.github.com>
Date: Fri, 21 Nov 2025 14:35:38 +0100
Subject: [PATCH 04/18] Update test/test-mtk.jl

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 test/test-mtk.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index f0b7f56..e0efc8d 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -159,6 +159,6 @@ using CairoMakie
         @test fig isa Figure
 
         # Optionally save the plot (commented out to avoid file creation in tests)
-        save("test_solar_position_plot.png", fig)
+        # save("test_solar_position_plot.png", fig)
     end
 end

From 0d356204b564d7a0642ae441412715d17f5b4761 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <37669773+langestefan@users.noreply.github.com>
Date: Fri, 21 Nov 2025 14:43:23 +0100
Subject: [PATCH 05/18] Update examples/mtk-integration.jl

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 examples/mtk-integration.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/examples/mtk-integration.jl b/examples/mtk-integration.jl
index 1744639..5f7eb94 100644
--- a/examples/mtk-integration.jl
+++ b/examples/mtk-integration.jl
@@ -20,13 +20,13 @@
 # #
 # # Create a solar position component for a specific location
 
-# @named sun_sf = SolarPositionBlock(
-#     latitude = 37.7749,      # San Francisco
-#     longitude = -122.4194,
-#     altitude = 100.0,
-#     t0 = DateTime(2024, 6, 21, 0, 0, 0),  # Summer solstice midnight
-#     time_unit = Hour(1),      # t=1.0 means 1 hour from t0
-# )
+# @named sun_sf = SolarPositionBlock()
+# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
+# t0 = DateTime(2024, 6, 21, 0, 0, 0)        # Summer solstice midnight
+# sys = mtkcompile(sun_sf)
+# pmap = [sys.observer => obs, sys.t0 => t0]
+#
+# # Now you can use `pmap` when simulating or evaluating the system
 
 # # The solar component has variables for azimuth, elevation, and zenith angles
 # println("Solar component unknowns: ", unknowns(sun_sf))

From 2353a3953050c1f17b036d820efec80ba380afba Mon Sep 17 00:00:00 2001
From: Stefan de Lange <37669773+langestefan@users.noreply.github.com>
Date: Fri, 21 Nov 2025 14:43:54 +0100
Subject: [PATCH 06/18] Update examples/mtk-integration.jl

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 examples/mtk-integration.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/mtk-integration.jl b/examples/mtk-integration.jl
index 5f7eb94..b170a81 100644
--- a/examples/mtk-integration.jl
+++ b/examples/mtk-integration.jl
@@ -175,7 +175,7 @@
 
 # @named room = building_room()
 # @named window = solar_window()
-# @named sun_building = SolarPositionBlock(latitude = 40.7128, longitude = -74.0060)  # NYC
+# @named sun_building = SolarPositionBlock()  # Set latitude/longitude via parameter map after compilation
 
 # # Connect the window solar gain to the room
 # @named building = compose(

From 4fa1382065db24b931f31414001891c44c1ba5ba Mon Sep 17 00:00:00 2001
From: Stefan de Lange <37669773+langestefan@users.noreply.github.com>
Date: Fri, 21 Nov 2025 14:44:11 +0100
Subject: [PATCH 07/18] Update test/test-mtk.jl

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 test/test-mtk.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index e0efc8d..b89e86b 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -1,8 +1,8 @@
 using Test
 using SolarPosition:
-    Observer, solar_position, PSA, NoRefraction, SolarAlgorithm, SolarPositionBlock, BENNETT
+    Observer, PSA, NoRefraction, SolarPositionBlock
 using ModelingToolkit:
-    @named, @variables, @parameters, get_variables, unknowns, System, mtkcompile, observed
+    @named, @variables, @parameters, unknowns, System, mtkcompile
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Dates: DateTime
 using OrdinaryDiffEq

From 6225851be82bcccabeb8d8e90df66d8617e481c2 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 22 Nov 2025 15:13:26 +0000
Subject: [PATCH 08/18] Initial plan


From 4acb563941540288f5b0488d89bad33c54ab447f Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Sat, 22 Nov 2025 15:17:45 +0000
Subject: [PATCH 09/18] Run Julia formatter to fix linting

Co-authored-by: langestefan <37669773+langestefan@users.noreply.github.com>
---
 test/test-mtk.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index b89e86b..4aa07c2 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -1,8 +1,6 @@
 using Test
-using SolarPosition:
-    Observer, PSA, NoRefraction, SolarPositionBlock
-using ModelingToolkit:
-    @named, @variables, @parameters, unknowns, System, mtkcompile
+using SolarPosition: Observer, PSA, NoRefraction, SolarPositionBlock
+using ModelingToolkit: @named, @variables, @parameters, unknowns, System, mtkcompile
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Dates: DateTime
 using OrdinaryDiffEq

From b1b00a3ccc8dcf3d5fe47999a54da9048ceb2178 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 20:00:47 +0100
Subject: [PATCH 10/18] Add apparent position to MTK component and validate
 positions

---
 ext/SolarPositionModelingToolkitExt.jl | 28 ++++++++---
 test/test-mtk.jl                       | 70 +++++++++++++++++++++++++-
 2 files changed, 89 insertions(+), 9 deletions(-)

diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
index 4a9153f..b3e0e64 100644
--- a/ext/SolarPositionModelingToolkitExt.jl
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -1,7 +1,9 @@
 module SolarPositionModelingToolkitExt
 
 using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
+using SolarPosition: SolPos, ApparentSolPos, SPASolPos
 using ModelingToolkit: @parameters, @variables, System, t_nounits
+using ModelingToolkit: t_nounits as t
 using Dates: DateTime, Millisecond
 import Symbolics
 
@@ -12,8 +14,18 @@ seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1
 
 # helper functions to extract fields from solar position
 get_azimuth(pos) = pos.azimuth
-get_elevation(pos) = pos.elevation
-get_zenith(pos) = pos.zenith
+
+# for SolPos: use elevation
+get_elevation(pos::SolPos) = pos.elevation
+get_zenith(pos::SolPos) = pos.zenith
+
+# for ApparentSolPos: use apparent_elevation and apparent_zenith
+get_elevation(pos::ApparentSolPos) = pos.apparent_elevation
+get_zenith(pos::ApparentSolPos) = pos.apparent_zenith
+
+# for SPASolPos: use apparent_elevation and apparent_zenith
+get_elevation(pos::SPASolPos) = pos.apparent_elevation
+get_zenith(pos::SPASolPos) = pos.apparent_zenith
 
 Symbolics.@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
 Symbolics.@register_symbolic solar_position(
@@ -22,21 +34,21 @@ Symbolics.@register_symbolic solar_position(
     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_nounits) [output = true]
-    @variables elevation(t_nounits) [output = true]
-    @variables zenith(t_nounits) [output = true]
+    @variables azimuth(t) [output = true]
+    @variables elevation(t) [output = true]
+    @variables zenith(t) [output = true]
 
-    time_expr = Symbolics.term(seconds_to_datetime, t_nounits, t0; type = DateTime)
+    time_expr = Symbolics.term(seconds_to_datetime, t, t0; type = DateTime)
     pos = solar_position(observer, time_expr, algorithm, refraction)
 
     eqs = [
@@ -45,7 +57,7 @@ function SolarPositionBlock(; name)
         zenith ~ get_zenith(pos),
     ]
 
-    return System(eqs, t_nounits; name = name)
+    return System(eqs, t; name = name)
 end
 
 end # module SolarPositionModelingToolkitExt
diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index 4aa07c2..bef3d0b 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -1,11 +1,20 @@
 using Test
-using SolarPosition: Observer, PSA, NoRefraction, SolarPositionBlock
+using SolarPosition:
+    Observer, PSA, NOAA, Walraven, USNO, SPA, NoRefraction, SolarPositionBlock
+using SolarPosition: HUGHES, BENNETT, ARCHER, MICHALSKY, SG2
 using ModelingToolkit: @named, @variables, @parameters, unknowns, System, mtkcompile
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Dates: DateTime
 using OrdinaryDiffEq
 using CairoMakie
 
+# Include helper functions from other test files for expected values
+include("test-psa.jl")
+include("test-noaa.jl")
+include("test-walraven.jl")
+include("test-usno.jl")
+include("test-spa.jl")
+
 @testset "ModelingToolkit Extension" begin
     obs = Observer(37.7749, -122.4194, 100.0)
     t0 = DateTime(2024, 6, 21, 0, 0, 0)  # Summer solstice
@@ -159,4 +168,63 @@ using CairoMakie
         # Optionally save the plot (commented out to avoid file creation in tests)
         # save("test_solar_position_plot.png", fig)
     end
+
+    @testset "Algorithm Accuracy" begin
+        # Test that MTK extension produces correct results for each algorithm
+        # Test all conditions from each algorithm's test file
+        conds = test_conditions()
+
+        @testset "$alg_name" for (alg_name, alg, exp_func, refr) in [
+            ("PSA", PSA(2020), expected_2020, NoRefraction()),
+            ("Walraven", Walraven(), expected_walraven, NoRefraction()),
+            ("USNO", USNO(), expected_usno, NoRefraction()),
+            ("NOAA", NOAA(), expected_noaa, HUGHES(101325.0, 10.0)),
+            ("SPA", SPA(), expected_spa, NoRefraction()),
+        ]
+            # Get expected values for all test cases
+            df_expected = exp_func()
+
+            for (i, ((zdt, lat, lon, alt), row)) in
+                enumerate(zip(eachrow(conds), eachrow(df_expected)))
+                # Convert ZonedDateTime to UTC DateTime for ModelingToolkit
+                # astimezone converts to UTC, then DateTime extracts the DateTime part
+                dt = DateTime(astimezone(zdt, tz"UTC"))
+
+                # Create observer
+                if ismissing(alt)
+                    obs = Observer(lat, lon)
+                else
+                    obs = Observer(lat, lon, altitude = alt)
+                end
+
+                @named sun = SolarPositionBlock()
+                sys = mtkcompile(sun)
+
+                pmap = [
+                    sys.observer => obs,
+                    sys.t0 => dt,
+                    sys.algorithm => alg,
+                    sys.refraction => refr,
+                ]
+
+                prob = ODEProblem(sys, pmap, (0.0, 1.0))
+                sol = solve(prob; saveat = [0.0])
+
+                # For algorithms with refraction or SPA, use apparent_elevation/apparent_zenith
+                # Otherwise use elevation/zenith
+                if haskey(row, :apparent_elevation)
+                    @test isapprox(
+                        sol[sys.elevation][1],
+                        row.apparent_elevation,
+                        atol = 1e-6,
+                    )
+                    @test isapprox(sol[sys.zenith][1], row.apparent_zenith, atol = 1e-6)
+                else
+                    @test isapprox(sol[sys.elevation][1], row.elevation, atol = 1e-6)
+                    @test isapprox(sol[sys.zenith][1], row.zenith, atol = 1e-6)
+                end
+                @test isapprox(sol[sys.azimuth][1], row.azimuth, atol = 1e-6)
+            end
+        end
+    end
 end

From b3e4ba66a6780c22e1b62fe6d51ed79f28437cb0 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 20:01:17 +0100
Subject: [PATCH 11/18] Add maxlog = 1 to @warn

---
 src/Positioning/Positioning.jl | 4 ++--
 src/Positioning/deltat.jl      | 2 +-
 src/Positioning/spa.jl         | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/Positioning/Positioning.jl b/src/Positioning/Positioning.jl
index 48b2225..c53e5d5 100644
--- a/src/Positioning/Positioning.jl
+++ b/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)
diff --git a/src/Positioning/deltat.jl b/src/Positioning/deltat.jl
index 56fcea2..7cbd478 100644
--- a/src/Positioning/deltat.jl
+++ b/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
diff --git a/src/Positioning/spa.jl b/src/Positioning/spa.jl
index d393a85..fa22aaa 100644
--- a/src/Positioning/spa.jl
+++ b/src/Positioning/spa.jl
@@ -84,10 +84,10 @@ struct SPAObserver{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)

From b580e188fe4aec585d5aecb869ee93aaaaa3c059 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 20:43:07 +0100
Subject: [PATCH 12/18] Put all expected values in expected-values.jl

---
 docs/src/examples/modelingtoolkit.md |   2 +
 test/expected-values.jl              | 209 +++++++++++++++++++++++++++
 test/setup.jl                        |   2 +
 test/test-mtk.jl                     |   7 -
 test/test-noaa.jl                    |  28 ----
 test/test-psa.jl                     |  56 -------
 test/test-spa.jl                     |  35 -----
 test/test-usno.jl                    |  56 -------
 test/test-walraven.jl                |  28 ----
 9 files changed, 213 insertions(+), 210 deletions(-)
 create mode 100644 test/expected-values.jl

diff --git a/docs/src/examples/modelingtoolkit.md b/docs/src/examples/modelingtoolkit.md
index 842f513..5497a4f 100644
--- a/docs/src/examples/modelingtoolkit.md
+++ b/docs/src/examples/modelingtoolkit.md
@@ -27,7 +27,9 @@ using ModelingToolkit
 using ModelingToolkit: t_nounits as t
 using Dates
 using OrdinaryDiffEq
+```
 
+```@example mtk
 # Create a solar position block
 @named sun = SolarPositionBlock()
 
diff --git a/test/expected-values.jl b/test/expected-values.jl
new file mode 100644
index 0000000..bfef360
--- /dev/null
+++ b/test/expected-values.jl
@@ -0,0 +1,209 @@
+"""Expected values for algorithm tests (shared helper functions)"""
+
+# PSA expected values
+function expected_2020()
+    columns = [:elevation, :zenith, :azimuth]
+
+    values = [
+        [32.21595946, 57.78404054, 204.9167547],
+        [32.20349382, 57.79650618, 204.96304003],
+        [34.92071906, 55.07928094, 169.37535422],
+        [18.63439275, 71.36560725, 234.1884117],
+        [35.7535716, 54.2464284, 197.67796274],
+        [-9.5321134, 99.5321134, 201.18736969],
+        [66.85741664, 23.14258336, 245.08558596],
+        [9.52730057, 80.47269943, 338.81263031],
+        [50.10853074, 39.89146926, 326.23536385],
+        [35.35979872, 54.64020128, 175.38781378],
+        [-53.24299848, 143.24299848, 18.64664538],
+        [-53.24299848, 143.24299848, 18.64664538],
+        [32.47215358, 57.52784642, 204.9305477],
+        [32.43510874, 57.56489126, 204.98638231],
+        [-23.41028925, 113.41028925, 79.54884709],
+        [1.10556787, 88.89443213, 104.54003062],
+        [32.21595946, 57.78404054, 204.9167547],
+        [32.21595946, 57.78404054, 204.9167547],
+        [32.21595946, 57.78404054, 204.9167547],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+function expected_2001()
+    columns = [:elevation, :zenith, :azimuth]
+
+    values = [
+        [32.21434385, 57.78565615, 204.91957868],
+        [32.20187689, 57.79812311, 204.96586258],
+        [34.92027644, 55.07972356, 169.3788305],
+        [18.63211707, 71.36788293, 234.19028111],
+        [35.75446859, 54.24553141, 197.67713395],
+        [-9.53293173, 99.53293173, 201.19017401],
+        [66.85455175, 23.14544825, 245.08643499],
+        [9.52811892, 80.47188108, 338.80982599],
+        [50.10817913, 39.89182087, 326.23090017],
+        [35.35914111, 54.64085889, 175.39125721],
+        [-53.24316093, 143.24316093, 18.65145716],
+        [-53.24316093, 143.24316093, 18.65145716],
+        [32.46428802, 57.53571198, 204.93415982],
+        [32.43981197, 57.56018803, 204.98955089],
+        [-23.40890746, 113.40890746, 79.55160745],
+        [1.10690714, 88.89309286, 104.54258663],
+        [32.21434385, 57.78565615, 204.91957868],
+        [32.21434385, 57.78565615, 204.91957868],
+        [32.21434385, 57.78565615, 204.91957868],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+# NOAA expected values
+function expected_noaa()
+    columns = [:elevation, :apparent_elevation, :zenith, :apparent_zenith, :azimuth]
+
+    values = [
+        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
+        [32.20468378, 32.23022967, 57.79531622, 57.76977033, 204.96860803],
+        [34.92392895, 34.94698595, 55.07607105, 55.05301405, 169.38094302],
+        [18.63438988, 18.68174893, 71.36561012, 71.31825107, 234.19290241],
+        [35.75618186, 35.77854313, 54.24381814, 54.22145687, 197.67357003],
+        [-9.52911395, -9.49473778, 99.52911395, 99.49473778, 201.19219097],
+        [66.85423515, 66.8611327, 23.14576485, 23.1388673, 245.09172279],
+        [9.52911581, 9.62132639, 80.47088419, 80.3786736, 338.80780907],
+        [50.10765752, 50.12113671, 39.89234248, 39.87886329, 326.22893168],
+        [35.36265374, 35.38534047, 54.63734626, 54.61465953, 175.39359304],
+        [-53.23987161, -53.23556094, 143.23987161, 143.23556094, 18.65415239],
+        [-53.23987161, -53.23556094, 143.23987161, 143.23556094, 18.65415239],
+        [32.46248831, 32.48778263, 57.53751169, 57.51221737, 204.95376627],
+        [32.44117331, 32.4664883, 57.55882669, 57.5335117, 204.97518601],
+        [-23.40444445, -23.39111232, 113.40444445, 113.39111232, 79.55187937],
+        [1.11161226, 1.46445929, 88.88838774, 88.53554071, 104.54291476],
+        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
+        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
+        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+# Walraven expected values
+function expected_walraven()
+    columns = [:elevation, :zenith, :azimuth]
+
+    values = [
+        [32.21577184, 57.78422816, 204.9145408],
+        [32.20330752, 57.79669248, 204.96082497],
+        [34.91988403, 55.08011597, 169.37445826],
+        [18.63514306, 71.36485694, 234.18579686],
+        [35.75942569, 54.24057431, 197.67543014],
+        [-9.53241956, 99.53241956, 201.18624892],
+        [66.85832643, 23.14167357, 245.07813388],
+        [9.53241956, 80.46758044, 338.81375108],
+        [50.11302395, 39.88697605, 326.23525905],
+        [35.35901685, 54.64098315, 175.38665251],
+        [-53.24443195, 143.24443195, 18.64588668],
+        [-53.24443195, 143.24443195, 18.64588668],
+        [33.19945919, 56.80054081, 205.08690939],
+        [32.078221, 57.921779, 204.90430273],
+        [-23.41074972, 113.41074972, 79.55008786],
+        [1.10528968, 88.89471032, 104.541125],
+        [32.21577184, 57.78422816, 204.9145408],
+        [32.21577184, 57.78422816, 204.9145408],
+        [32.21577184, 57.78422816, 204.9145408],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+# USNO expected values
+function expected_usno()
+    columns = [:elevation, :zenith, :azimuth]
+
+    values = [
+        [32.21913058, 57.78086942, 204.91396212],
+        [32.20666654, 57.79333346, 204.9602485],
+        [34.92271632, 55.07728368, 169.37217155],
+        [18.63848631, 71.36151369, 234.18640074],
+        [35.75823472, 54.24176528, 197.67107057],
+        [-9.52936557, 99.52936557, 201.18474698],
+        [66.86088833, 23.13911167, 245.08379957],
+        [9.52936557, 80.47063443, 338.81525302],
+        [50.11081313, 39.88918687, 326.23927513],
+        [35.36198031, 54.63801969, 175.38462327],
+        [-53.24179879, 143.24179879, 18.6423076],
+        [-53.24179879, 143.24179879, 18.6423076],
+        [32.46463884, 57.53536116, 204.9389006],
+        [32.44304455, 57.55695545, 204.98724112],
+        [-23.40963192, 113.40963192, 79.54658304],
+        [1.10645791, 88.89354209, 104.53792899],
+        [32.21913058, 57.78086942, 204.91396212],
+        [32.21913058, 57.78086942, 204.91396212],
+        [32.21913058, 57.78086942, 204.91396212],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+function expected_usno_option_2()
+    columns = [:elevation, :zenith, :azimuth]
+
+    values = [
+        [32.21913452, 57.78086548, 204.91394742],
+        [32.20667049, 57.79332951, 204.96023379],
+        [34.92271497, 55.07728503, 169.37215926],
+        [18.63849557, 71.36150443, 234.18638706],
+        [35.7582376, 54.2417624, 197.67105461],
+        [-9.52936557, 99.52936557, 201.18473375],
+        [66.86090033, 23.13909967, 245.08378652],
+        [9.52936557, 80.47063443, 338.81526625],
+        [50.11081833, 39.88918167, 326.2392938],
+        [35.36197956, 54.63802044, 175.38460729],
+        [-53.24180178, 143.24180178, 18.64228637],
+        [-53.24180178, 143.24180178, 18.64228637],
+        [32.46465875, 57.53534125, 204.93882613],
+        [32.44303542, 57.55696458, 204.98727521],
+        [-23.40963198, 113.40963198, 79.54658298],
+        [1.10645552, 88.89354448, 104.53792651],
+        [32.21913452, 57.78086548, 204.91394742],
+        [32.21913452, 57.78086548, 204.91394742],
+        [32.21913452, 57.78086548, 204.91394742],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
+
+# SPA expected values
+function expected_spa()
+    columns = [
+        :elevation,
+        :apparent_elevation,
+        :zenith,
+        :apparent_zenith,
+        :azimuth,
+        :equation_of_time,
+    ]
+
+    values = [
+        [32.21708674, 32.24367654, 57.78291326, 57.75632346, 204.9188164, 14.75816027],
+        [32.20462016, 32.23122263, 57.79537984, 57.76877737, 204.96510208, 14.75818295],
+        [34.92249857, 34.94652348, 55.07750143, 55.05347652, 169.3767211, 14.74179647],
+        [18.63493009, 18.68392203, 71.36506991, 71.31607797, 234.19049983, 14.77445424],
+        [35.75662074, 35.77992239, 54.24337926, 54.22007761, 197.67009996, -12.4319873],
+        [-9.53051391, -9.53051391, 99.53051391, 99.53051391, 201.18870641, 14.75816027],
+        [66.85682918, 66.86401769, 23.14317082, 23.13598231, 245.09065324, 14.75816027],
+        [9.52569499, 9.61953906, 80.47430501, 80.38046094, 338.81129359, 14.75816027],
+        [50.10653615, 50.12059933, 39.89346385, 39.87940067, 326.23446721, 14.75816027],
+        [35.36147343, 35.38511409, 54.63852657, 54.61488591, 175.38931352, 14.75816027],
+        [-53.24113451, -53.24113451, 143.24113451, 143.24113451, 18.64817122, 14.75816027],
+        [-53.24113451, -53.24113451, 143.24113451, 143.24113451, 18.64817122, 14.75816027],
+        [32.4622024, 32.48854464, 57.5377976, 57.51145536, 204.94926574, 14.55795391],
+        [32.4407131, 32.46707691, 57.5592869, 57.53292309, 204.96985758, 14.64779166],
+        [-23.40808954, -23.40808954, 113.40808954, 113.40808954, 79.54868063, 14.68397497],
+        [1.10773228, 1.45847482, 88.89226772, 88.54152518, 104.53989829, 14.70334338],
+        [32.21708674, 32.24367654, 57.78291326, 57.75632346, 204.9188164, 14.75816027],
+        [32.21708677, 32.24367657, 57.78291323, 57.75632343, 204.9188164, 14.75816027],
+        [32.21708544, 32.24367524, 57.78291456, 57.75632476, 204.9188164, 14.75816027],
+    ]
+
+    return DataFrame(reduce(hcat, values)', columns)
+end
diff --git a/test/setup.jl b/test/setup.jl
index d82f64c..e6536c6 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -4,6 +4,8 @@ using DataFrames
 using TimeZones
 using Dates
 
+include("expected-values.jl")
+
 function test_conditions()
     inputs = DataFrame(
         time = [
diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index bef3d0b..da48f3c 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -8,13 +8,6 @@ using Dates: DateTime
 using OrdinaryDiffEq
 using CairoMakie
 
-# Include helper functions from other test files for expected values
-include("test-psa.jl")
-include("test-noaa.jl")
-include("test-walraven.jl")
-include("test-usno.jl")
-include("test-spa.jl")
-
 @testset "ModelingToolkit Extension" begin
     obs = Observer(37.7749, -122.4194, 100.0)
     t0 = DateTime(2024, 6, 21, 0, 0, 0)  # Summer solstice
diff --git a/test/test-noaa.jl b/test/test-noaa.jl
index ae16bca..4dcfe85 100644
--- a/test/test-noaa.jl
+++ b/test/test-noaa.jl
@@ -1,33 +1,5 @@
 """Unit tests for NOAA.jl"""
 
-function expected_noaa()
-    columns = [:elevation, :apparent_elevation, :zenith, :apparent_zenith, :azimuth]
-
-    values = [
-        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
-        [32.20468378, 32.23022967, 57.79531622, 57.76977033, 204.96860803],
-        [34.92392895, 34.94698595, 55.07607105, 55.05301405, 169.38094302],
-        [18.63438988, 18.68174893, 71.36561012, 71.31825107, 234.19290241],
-        [35.75618186, 35.77854313, 54.24381814, 54.22145687, 197.67357003],
-        [-9.52911395, -9.49473778, 99.52911395, 99.49473778, 201.19219097],
-        [66.85423515, 66.8611327, 23.14576485, 23.1388673, 245.09172279],
-        [9.52911581, 9.62132639, 80.47088419, 80.3786736, 338.80780907],
-        [50.10765752, 50.12113671, 39.89234248, 39.87886329, 326.22893168],
-        [35.36265374, 35.38534047, 54.63734626, 54.61465953, 175.39359304],
-        [-53.23987161, -53.23556094, 143.23987161, 143.23556094, 18.65415239],
-        [-53.23987161, -53.23556094, 143.23987161, 143.23556094, 18.65415239],
-        [32.46248831, 32.48778263, 57.53751169, 57.51221737, 204.95376627],
-        [32.44117331, 32.4664883, 57.55882669, 57.5335117, 204.97518601],
-        [-23.40444445, -23.39111232, 113.40444445, 113.39111232, 79.55187937],
-        [1.11161226, 1.46445929, 88.88838774, 88.53554071, 104.54291476],
-        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
-        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
-        [32.21715174, 32.24268539, 57.78284826, 57.75731461, 204.92232395],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
 @testset "NOAA" begin
     df_expected = expected_noaa()
     conds = test_conditions()
diff --git a/test/test-psa.jl b/test/test-psa.jl
index 7116c31..805080a 100644
--- a/test/test-psa.jl
+++ b/test/test-psa.jl
@@ -1,61 +1,5 @@
 """Unit tests for PSA.jl"""
 
-function expected_2020()
-    columns = [:elevation, :zenith, :azimuth]
-
-    values = [
-        [32.21595946, 57.78404054, 204.9167547],
-        [32.20349382, 57.79650618, 204.96304003],
-        [34.92071906, 55.07928094, 169.37535422],
-        [18.63439275, 71.36560725, 234.1884117],
-        [35.7535716, 54.2464284, 197.67796274],
-        [-9.5321134, 99.5321134, 201.18736969],
-        [66.85741664, 23.14258336, 245.08558596],
-        [9.52730057, 80.47269943, 338.81263031],
-        [50.10853074, 39.89146926, 326.23536385],
-        [35.35979872, 54.64020128, 175.38781378],
-        [-53.24299848, 143.24299848, 18.64664538],
-        [-53.24299848, 143.24299848, 18.64664538],
-        [32.47215358, 57.52784642, 204.9305477],
-        [32.43510874, 57.56489126, 204.98638231],
-        [-23.41028925, 113.41028925, 79.54884709],
-        [1.10556787, 88.89443213, 104.54003062],
-        [32.21595946, 57.78404054, 204.9167547],
-        [32.21595946, 57.78404054, 204.9167547],
-        [32.21595946, 57.78404054, 204.9167547],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
-function expected_2001()
-    columns = [:elevation, :zenith, :azimuth]
-
-    values = [
-        [32.21434385, 57.78565615, 204.91957868],
-        [32.20187689, 57.79812311, 204.96586258],
-        [34.92027644, 55.07972356, 169.3788305],
-        [18.63211707, 71.36788293, 234.19028111],
-        [35.75446859, 54.24553141, 197.67713395],
-        [-9.53293173, 99.53293173, 201.19017401],
-        [66.85455175, 23.14544825, 245.08643499],
-        [9.52811892, 80.47188108, 338.80982599],
-        [50.10817913, 39.89182087, 326.23090017],
-        [35.35914111, 54.64085889, 175.39125721],
-        [-53.24316093, 143.24316093, 18.65145716],
-        [-53.24316093, 143.24316093, 18.65145716],
-        [32.46428802, 57.53571198, 204.93415982],
-        [32.43981197, 57.56018803, 204.98955089],
-        [-23.40890746, 113.40890746, 79.55160745],
-        [1.10690714, 88.89309286, 104.54258663],
-        [32.21434385, 57.78565615, 204.91957868],
-        [32.21434385, 57.78565615, 204.91957868],
-        [32.21434385, 57.78565615, 204.91957868],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
 @testset "PSA" begin
     coeffs = Dict(2020 => expected_2020, 2001 => expected_2001)
 
diff --git a/test/test-spa.jl b/test/test-spa.jl
index 3f6289e..8cb6a77 100644
--- a/test/test-spa.jl
+++ b/test/test-spa.jl
@@ -1,40 +1,5 @@
 """Unit tests for SPA algorithm"""
 
-function expected_spa()
-    columns = [
-        :elevation,
-        :apparent_elevation,
-        :zenith,
-        :apparent_zenith,
-        :azimuth,
-        :equation_of_time,
-    ]
-
-    values = [
-        [32.21708674, 32.24367654, 57.78291326, 57.75632346, 204.9188164, 14.75816027],
-        [32.20462016, 32.23122263, 57.79537984, 57.76877737, 204.96510208, 14.75818295],
-        [34.92249857, 34.94652348, 55.07750143, 55.05347652, 169.3767211, 14.74179647],
-        [18.63493009, 18.68392203, 71.36506991, 71.31607797, 234.19049983, 14.77445424],
-        [35.75662074, 35.77992239, 54.24337926, 54.22007761, 197.67009996, -12.4319873],
-        [-9.53051391, -9.53051391, 99.53051391, 99.53051391, 201.18870641, 14.75816027],
-        [66.85682918, 66.86401769, 23.14317082, 23.13598231, 245.09065324, 14.75816027],
-        [9.52569499, 9.61953906, 80.47430501, 80.38046094, 338.81129359, 14.75816027],
-        [50.10653615, 50.12059933, 39.89346385, 39.87940067, 326.23446721, 14.75816027],
-        [35.36147343, 35.38511409, 54.63852657, 54.61488591, 175.38931352, 14.75816027],
-        [-53.24113451, -53.24113451, 143.24113451, 143.24113451, 18.64817122, 14.75816027],
-        [-53.24113451, -53.24113451, 143.24113451, 143.24113451, 18.64817122, 14.75816027],
-        [32.4622024, 32.48854464, 57.5377976, 57.51145536, 204.94926574, 14.55795391],
-        [32.4407131, 32.46707691, 57.5592869, 57.53292309, 204.96985758, 14.64779166],
-        [-23.40808954, -23.40808954, 113.40808954, 113.40808954, 79.54868063, 14.68397497],
-        [1.10773228, 1.45847482, 88.89226772, 88.54152518, 104.53989829, 14.70334338],
-        [32.21708674, 32.24367654, 57.78291326, 57.75632346, 204.9188164, 14.75816027],
-        [32.21708677, 32.24367657, 57.78291323, 57.75632343, 204.9188164, 14.75816027],
-        [32.21708544, 32.24367524, 57.78291456, 57.75632476, 204.9188164, 14.75816027],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
 @testset "SPA" begin
     df_expected = expected_spa()
     conds = test_conditions()
diff --git a/test/test-usno.jl b/test/test-usno.jl
index 1b76003..732b920 100644
--- a/test/test-usno.jl
+++ b/test/test-usno.jl
@@ -1,61 +1,5 @@
 """Unit tests for USNO algorithm"""
 
-function expected_usno()
-    columns = [:elevation, :zenith, :azimuth]
-
-    values = [
-        [32.21913058, 57.78086942, 204.91396212],
-        [32.20666654, 57.79333346, 204.9602485],
-        [34.92271632, 55.07728368, 169.37217155],
-        [18.63848631, 71.36151369, 234.18640074],
-        [35.75823472, 54.24176528, 197.67107057],
-        [-9.52936557, 99.52936557, 201.18474698],
-        [66.86088833, 23.13911167, 245.08379957],
-        [9.52936557, 80.47063443, 338.81525302],
-        [50.11081313, 39.88918687, 326.23927513],
-        [35.36198031, 54.63801969, 175.38462327],
-        [-53.24179879, 143.24179879, 18.6423076],
-        [-53.24179879, 143.24179879, 18.6423076],
-        [32.46463884, 57.53536116, 204.9389006],
-        [32.44304455, 57.55695545, 204.98724112],
-        [-23.40963192, 113.40963192, 79.54658304],
-        [1.10645791, 88.89354209, 104.53792899],
-        [32.21913058, 57.78086942, 204.91396212],
-        [32.21913058, 57.78086942, 204.91396212],
-        [32.21913058, 57.78086942, 204.91396212],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
-function expected_usno_option_2()
-    columns = [:elevation, :zenith, :azimuth]
-
-    values = [
-        [32.21913452, 57.78086548, 204.91394742],
-        [32.20667049, 57.79332951, 204.96023379],
-        [34.92271497, 55.07728503, 169.37215926],
-        [18.63849557, 71.36150443, 234.18638706],
-        [35.7582376, 54.2417624, 197.67105461],
-        [-9.52936557, 99.52936557, 201.18473375],
-        [66.86090033, 23.13909967, 245.08378652],
-        [9.52936557, 80.47063443, 338.81526625],
-        [50.11081833, 39.88918167, 326.2392938],
-        [35.36197956, 54.63802044, 175.38460729],
-        [-53.24180178, 143.24180178, 18.64228637],
-        [-53.24180178, 143.24180178, 18.64228637],
-        [32.46465875, 57.53534125, 204.93882613],
-        [32.44303542, 57.55696458, 204.98727521],
-        [-23.40963198, 113.40963198, 79.54658298],
-        [1.10645552, 88.89354448, 104.53792651],
-        [32.21913452, 57.78086548, 204.91394742],
-        [32.21913452, 57.78086548, 204.91394742],
-        [32.21913452, 57.78086548, 204.91394742],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
 @testset "USNO" begin
     @testset "With default delta_t (gmst_option=1)" begin
         df_expected = expected_usno()
diff --git a/test/test-walraven.jl b/test/test-walraven.jl
index 5949bb5..cb6aa69 100644
--- a/test/test-walraven.jl
+++ b/test/test-walraven.jl
@@ -1,33 +1,5 @@
 """Unit tests for Walraven algorithm"""
 
-function expected_walraven()
-    columns = [:elevation, :zenith, :azimuth]
-
-    values = [
-        [32.21577184, 57.78422816, 204.9145408],
-        [32.20330752, 57.79669248, 204.96082497],
-        [34.91988403, 55.08011597, 169.37445826],
-        [18.63514306, 71.36485694, 234.18579686],
-        [35.75942569, 54.24057431, 197.67543014],
-        [-9.53241956, 99.53241956, 201.18624892],
-        [66.85832643, 23.14167357, 245.07813388],
-        [9.53241956, 80.46758044, 338.81375108],
-        [50.11302395, 39.88697605, 326.23525905],
-        [35.35901685, 54.64098315, 175.38665251],
-        [-53.24443195, 143.24443195, 18.64588668],
-        [-53.24443195, 143.24443195, 18.64588668],
-        [33.19945919, 56.80054081, 205.08690939],
-        [32.078221, 57.921779, 204.90430273],
-        [-23.41074972, 113.41074972, 79.55008786],
-        [1.10528968, 88.89471032, 104.541125],
-        [32.21577184, 57.78422816, 204.9145408],
-        [32.21577184, 57.78422816, 204.9145408],
-        [32.21577184, 57.78422816, 204.9145408],
-    ]
-
-    return DataFrame(reduce(hcat, values)', columns)
-end
-
 @testset "Walraven" begin
     df_expected = expected_walraven()
     conds = test_conditions()

From 2cbc6b29d02f01960058e2c2cc428831a0827ff0 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 20:58:19 +0100
Subject: [PATCH 13/18] Fix typo

---
 docs/src/examples/modelingtoolkit.md | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/docs/src/examples/modelingtoolkit.md b/docs/src/examples/modelingtoolkit.md
index 5497a4f..4ab7973 100644
--- a/docs/src/examples/modelingtoolkit.md
+++ b/docs/src/examples/modelingtoolkit.md
@@ -216,5 +216,4 @@ engine, so its internals cannot be symbolically simplified.
 
 - [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
+- [ModelingToolkit.jl Documentation](https://docs.sciml.ai/ModelingToolkit/stable/) - MTK framework documentation

From d0ccd85a183ffc3d460ad3a441a7fcc3a73267e3 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 21:07:23 +0100
Subject: [PATCH 14/18] Use abstract solpos types in MTK ext

---
 ext/SolarPositionModelingToolkitExt.jl | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
index b3e0e64..736ccf8 100644
--- a/ext/SolarPositionModelingToolkitExt.jl
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -1,7 +1,7 @@
 module SolarPositionModelingToolkitExt
 
 using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
-using SolarPosition: SolPos, ApparentSolPos, SPASolPos
+using SolarPosition: SolPos, ApparentSolPos, SPASolPos, AbstractApparentSolPos
 using ModelingToolkit: @parameters, @variables, System, t_nounits
 using ModelingToolkit: t_nounits as t
 using Dates: DateTime, Millisecond
@@ -15,17 +15,13 @@ seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1
 # helper functions to extract fields from solar position
 get_azimuth(pos) = pos.azimuth
 
-# for SolPos: use elevation
+# for SolPos: use elevation and zenith
 get_elevation(pos::SolPos) = pos.elevation
 get_zenith(pos::SolPos) = pos.zenith
 
-# for ApparentSolPos: use apparent_elevation and apparent_zenith
-get_elevation(pos::ApparentSolPos) = pos.apparent_elevation
-get_zenith(pos::ApparentSolPos) = pos.apparent_zenith
-
-# for SPASolPos: use apparent_elevation and apparent_zenith
-get_elevation(pos::SPASolPos) = pos.apparent_elevation
-get_zenith(pos::SPASolPos) = pos.apparent_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(

From 2d8ce0fe991866a8b78883b039627f9567c7c87b Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 21:08:46 +0100
Subject: [PATCH 15/18] Remove whitespace

---
 docs/src/examples/modelingtoolkit.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/examples/modelingtoolkit.md b/docs/src/examples/modelingtoolkit.md
index 4ab7973..95cde60 100644
--- a/docs/src/examples/modelingtoolkit.md
+++ b/docs/src/examples/modelingtoolkit.md
@@ -62,7 +62,7 @@ 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
+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

From d001a0c14be6abaf29c0a4365d41da4e24be3e3c Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 21:13:10 +0100
Subject: [PATCH 16/18] Remove unecessary enumerate

---
 test/test-mtk.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/test-mtk.jl b/test/test-mtk.jl
index da48f3c..b4f07dc 100644
--- a/test/test-mtk.jl
+++ b/test/test-mtk.jl
@@ -177,8 +177,7 @@ using CairoMakie
             # Get expected values for all test cases
             df_expected = exp_func()
 
-            for (i, ((zdt, lat, lon, alt), row)) in
-                enumerate(zip(eachrow(conds), eachrow(df_expected)))
+            for ((zdt, lat, lon, alt), row) in zip(eachrow(conds), eachrow(df_expected))
                 # Convert ZonedDateTime to UTC DateTime for ModelingToolkit
                 # astimezone converts to UTC, then DateTime extracts the DateTime part
                 dt = DateTime(astimezone(zdt, tz"UTC"))

From 35877b1d8dbd8d700bde9f125bfc921b2ba02ed2 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <37669773+langestefan@users.noreply.github.com>
Date: Sun, 23 Nov 2025 21:44:10 +0100
Subject: [PATCH 17/18] Update ext/SolarPositionModelingToolkitExt.jl

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 ext/SolarPositionModelingToolkitExt.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/SolarPositionModelingToolkitExt.jl b/ext/SolarPositionModelingToolkitExt.jl
index 736ccf8..5111736 100644
--- a/ext/SolarPositionModelingToolkitExt.jl
+++ b/ext/SolarPositionModelingToolkitExt.jl
@@ -2,7 +2,7 @@ module SolarPositionModelingToolkitExt
 
 using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
 using SolarPosition: SolPos, ApparentSolPos, SPASolPos, AbstractApparentSolPos
-using ModelingToolkit: @parameters, @variables, System, t_nounits
+using ModelingToolkit: @parameters, @variables, System
 using ModelingToolkit: t_nounits as t
 using Dates: DateTime, Millisecond
 import Symbolics

From 5e1824bd7aa033291110266e54acf9a6cae7ab49 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <langestefan@msn.com>
Date: Sun, 23 Nov 2025 23:58:49 +0100
Subject: [PATCH 18/18] Cleanup

---
 docs/src/examples/modelingtoolkit.md |   2 +-
 examples/mtk-integration.jl          | 314 ---------------------------
 src/SolarPosition.jl                 |   3 +-
 3 files changed, 3 insertions(+), 316 deletions(-)
 delete mode 100644 examples/mtk-integration.jl

diff --git a/docs/src/examples/modelingtoolkit.md b/docs/src/examples/modelingtoolkit.md
index 95cde60..cccf29b 100644
--- a/docs/src/examples/modelingtoolkit.md
+++ b/docs/src/examples/modelingtoolkit.md
@@ -184,7 +184,7 @@ pmap = [
 ]
 
 prob = ODEProblem(sys_building, pmap, (0.0, 86400.0))
-sol = solve(prob, Tsit5(); saveat = 600.0)
+sol = solve(prob; saveat = 600.0)
 
 # Plot temperature evolution
 fig = Figure(size = (1200, 400))
diff --git a/examples/mtk-integration.jl b/examples/mtk-integration.jl
deleted file mode 100644
index b170a81..0000000
--- a/examples/mtk-integration.jl
+++ /dev/null
@@ -1,314 +0,0 @@
-# # # ModelingToolkit Integration Example
-# #
-# # This example demonstrates how to use the SolarPosition.jl ModelingToolkit extension
-# # to integrate solar position calculations into larger modeling workflows.
-# #
-# # The ModelingToolkit extension allows you to:
-# # 1. Create reusable solar position components
-# # 2. Compose solar models with other physical systems
-# # 3. Leverage MTK's symbolic simplification and code generation
-# # 4. Build complex solar energy system models
-
-# using ModelingToolkit
-# using ModelingToolkit: t_nounits as t, D_nounits as D
-# using SolarPosition
-# using Dates
-# using OrdinaryDiffEq
-# using Plots
-
-# # ## Example 1: Basic Solar Component
-# #
-# # Create a solar position component for a specific location
-
-# @named sun_sf = SolarPositionBlock()
-# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
-# t0 = DateTime(2024, 6, 21, 0, 0, 0)        # Summer solstice midnight
-# sys = mtkcompile(sun_sf)
-# pmap = [sys.observer => obs, sys.t0 => t0]
-#
-# # Now you can use `pmap` when simulating or evaluating the system
-
-# # The solar component has variables for azimuth, elevation, and zenith angles
-# println("Solar component unknowns: ", unknowns(sun_sf))
-# println("Solar component parameters: ", parameters(sun_sf))
-# println("Solar component equations: ", equations(sun_sf))
-
-# # ## Example 2: Solar Panel Power Model
-# #
-# # Build a simple solar panel model that uses the solar position
-
-# function simple_solar_panel(; name)
-#     @parameters begin
-#         area = 10.0           # Panel area in m²
-#         efficiency = 0.20     # Panel efficiency (20%)
-#         tilt_angle = 30.0     # Panel tilt from horizontal (degrees)
-#         azimuth_angle = 180.0 # Panel azimuth (180° = South facing)
-#     end
-
-#     @variables begin
-#         power(t)              # Power output in W
-#         dni(t) = 800.0        # Direct Normal Irradiance in W/m²
-#         incidence_angle(t)    # Angle between sun and panel normal
-#         effective_irr(t)      # Effective irradiance on panel
-#     end
-
-#     # Note: This is a simplified model for demonstration
-#     # In a real model, you would connect sun.elevation and sun.azimuth
-#     # to calculate the actual incidence angle
-#     eqs = [
-#         # Simplified: assume some base irradiance pattern
-#         dni ~ 800.0 + 200.0 * sin(t * π / 12),  # Daily variation
-#         # Effective irradiance (simplified)
-#         effective_irr ~ dni * max(0.0, cos(deg2rad(45.0))),
-#         # Power output
-#         power ~ area * efficiency * effective_irr,
-#     ]
-
-#     return System(eqs, t; name = name)
-# end
-
-# @named panel = simple_solar_panel()
-
-# # Compose the solar component with the panel model
-# @named solar_system = compose(System(Equation[], t; name = :solar_system), sun_sf, panel)
-
-# println("\nComposed system unknowns: ", length(unknowns(solar_system)))
-# println("Composed system equations: ", length(equations(solar_system)))
-
-# # ## Example 3: Working Simulation with Solar Panel
-# #
-# # Create a complete working example that can be compiled and simulated
-
-# function solar_panel_with_sun(; name, sun_elevation)
-#     @parameters begin
-#         area = 10.0           # Panel area in m²
-#         efficiency = 0.20     # Panel efficiency (20%)
-#         base_irradiance = 1000.0  # Base solar irradiance W/m²
-#     end
-
-#     @variables begin
-#         power(t)              # Power output in W
-#         irradiance(t)         # Effective irradiance on panel
-#     end
-
-#     eqs = [
-#         # Irradiance depends on sun elevation
-#         irradiance ~ base_irradiance * max(0.0, sin(deg2rad(sun_elevation))),
-#         # Power output
-#         power ~ area * efficiency * irradiance,
-#     ]
-
-#     return System(eqs, t; name = name)
-# end
-
-# # Pre-compute solar position for a specific time
-# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
-# ref_time = DateTime(2024, 6, 21, 20, 0, 0)  # Summer solstice, local solar noon (UTC)
-# ref_pos = solar_position(obs, ref_time, PSA())
-
-# println("\nReference solar position at summer solstice noon:")
-# println("  Elevation: $(ref_pos.elevation)°")
-# println("  Azimuth: $(ref_pos.azimuth)°")
-
-# # Create solar panel with the computed elevation
-# @named panel_working = solar_panel_with_sun(sun_elevation = ref_pos.elevation)
-
-# # Compile and simulate
-# panel_sys = mtkcompile(panel_working)
-# prob = ODEProblem(panel_sys, Dict(), (0.0, 10.0))
-# sol = solve(prob, Tsit5())
-
-# println("\nSolar panel simulation:")
-# println("  Power output: $(sol[panel_working.power][end]) W")
-# println("  Irradiance: $(sol[panel_working.irradiance][end]) W/m²")
-
-# # ## Example 4: Building Thermal Model with Solar Gain
-# #
-# # Model a simple building with solar heat gain through windows
-
-# function building_room(; name)
-#     @parameters begin
-#         mass = 1000.0           # Thermal mass in kg
-#         cp = 1000.0             # Specific heat capacity J/(kg·K)
-#         U_wall = 0.5            # Wall U-value W/(m²·K)
-#         wall_area = 50.0        # Wall area m²
-#         T_outside = 20.0        # Outside temperature °C
-#     end
-
-#     @variables begin
-#         T(t) = 20.0             # Room temperature °C
-#         Q_loss(t)               # Heat loss through walls W
-#         Q_solar(t) = 0.0        # Solar heat gain W (placeholder)
-#     end
-
-#     eqs = [
-#         # Heat loss through walls
-#         Q_loss ~ U_wall * wall_area * (T - T_outside),
-#         # Energy balance (simplified)
-#         D(T) ~ (Q_solar - Q_loss) / (mass * cp),
-#     ]
-
-#     return System(eqs, t; name = name)
-# end
-
-# function solar_window(; name)
-#     @parameters begin
-#         window_area = 5.0       # Window area m²
-#         transmittance = 0.7     # Glass transmittance
-#         normal_azimuth = 180.0  # South-facing
-#     end
-
-#     @variables begin
-#         Q_solar(t)              # Solar heat gain W
-#         irradiance(t)           # Solar irradiance W/m²
-#     end
-
-#     # Simplified solar gain calculation
-#     eqs = [
-#         # Base irradiance pattern (in real model, use sun.elevation)
-#         irradiance ~ 500.0 * (1 + sin(t * π / 12)),
-#         Q_solar ~ window_area * transmittance * irradiance,
-#     ]
-
-#     return System(eqs, t; name = name)
-# end
-
-# @named room = building_room()
-# @named window = solar_window()
-# @named sun_building = SolarPositionBlock()  # Set latitude/longitude via parameter map after compilation
-
-# # Connect the window solar gain to the room
-# @named building = compose(
-#     System([room.Q_solar ~ window.Q_solar], t; name = :building),
-#     room,
-#     window,
-#     sun_building,
-# )
-
-# println("\nBuilding model unknowns: ", length(unknowns(building)))
-# println("Building model parameters: ", length(parameters(building)))
-
-# # ## Example 5: Time-Varying Solar Simulation
-# #
-# # Simulate a simple system over a day with changing solar position
-
-# function simple_thermal_mass(; name, Q_solar_func)
-#     @parameters begin
-#         mass = 100.0        # Thermal mass in kg
-#         cp = 1000.0         # Specific heat J/(kg·K)
-#         T_amb = 20.0        # Ambient temperature °C
-#         h = 5.0             # Heat transfer coefficient W/(m²·K)
-#         area = 1.0          # Surface area m²
-#     end
-
-#     @variables begin
-#         T(t) = 20.0         # Temperature °C
-#         Q_solar(t)          # Solar heat input W
-#     end
-
-#     eqs = [Q_solar ~ Q_solar_func, D(T) ~ (Q_solar - h * area * (T - T_amb)) / (mass * cp)]
-
-#     return System(eqs, t; name = name)
-# end
-
-# # Create a time-varying solar input (simplified as sinusoidal)
-# # In practice, you'd use pre-computed solar positions with interpolation
-# solar_input_pattern = 500.0 * max(0.0, sin(π * t / 12))  # Peaks at noon (t=6)
-
-# @named thermal = simple_thermal_mass(Q_solar_func = solar_input_pattern)
-# thermal_sys = mtkcompile(thermal)
-
-# # Simulate over 24 hours
-# prob_thermal = ODEProblem(thermal_sys, Dict(), (0.0, 24.0))
-# sol_thermal = solve(prob_thermal, Tsit5())
-
-# println("\nThermal mass simulation with solar input:")
-# println("  Initial temperature: $(sol_thermal[thermal.T][1]) °C")
-# println("  Final temperature: $(sol_thermal[thermal.T][end]) °C")
-# println("  Peak temperature: $(maximum(sol_thermal[thermal.T])) °C")
-
-# # Plot temperature evolution
-# p_thermal = plot(
-#     sol_thermal,
-#     idxs = [thermal.T],
-#     xlabel = "Time (hours)",
-#     ylabel = "Temperature (°C)",
-#     label = "Temperature",
-#     title = "Thermal Mass with Solar Heating",
-#     linewidth = 2,
-#     legend = :best,
-# )
-
-# display(p_thermal)
-
-# # ## Example 6: Using Real Solar Position Data
-# #
-# # For accurate solar position calculations, you can compute positions offline
-# # and use them as time-dependent parameters or callbacks
-
-# # Create observer and time range
-# obs = Observer(37.7749, -122.4194, 100.0)  # San Francisco
-# start_time = DateTime(2024, 6, 21, 0, 0, 0)  # Summer solstice
-# times = [start_time + Hour(h) for h = 0:23]
-
-# # Compute solar positions
-# positions = solar_position(obs, times, PSA())
-
-# # Extract data
-# hours = 0:23
-# azimuths = positions.azimuth
-# elevations = positions.elevation
-# zeniths = positions.zenith
-
-# # Plot solar position over the day
-# p3 = plot(
-#     hours,
-#     elevations,
-#     xlabel = "Hour of Day",
-#     ylabel = "Solar Elevation (°)",
-#     label = "Elevation Angle",
-#     title = "Solar Position - Summer Solstice, San Francisco",
-#     linewidth = 2,
-#     marker = :circle,
-#     legend = :best,
-# )
-
-# p4 = plot(
-#     hours,
-#     azimuths,
-#     xlabel = "Hour of Day",
-#     ylabel = "Solar Azimuth (°)",
-#     label = "Azimuth Angle",
-#     title = "Solar Azimuth Over the Day",
-#     linewidth = 2,
-#     marker = :circle,
-#     legend = :best,
-# )
-
-# plot(p3, p4, layout = (2, 1), size = (800, 600))
-
-# # ## Notes for Advanced Usage
-# #
-# # 1. **Time Mapping**: In a real application, you need to map simulation time `t`
-# #    to actual DateTime values. This can be done with callbacks or custom functions.
-# #
-# # 2. **Callbacks**: For dynamic solar position updates during simulation, implement
-# #    a DiscreteCallback that updates solar position parameters at each time step.
-# #
-# # 3. **Performance**: For long simulations, consider pre-computing solar positions
-# #    and using interpolation rather than calculating at each solver step.
-# #
-# # 4. **Refraction**: Include atmospheric refraction for more accurate apparent
-# #    solar positions, especially important near sunrise/sunset.
-# #
-# # 5. **Integration**: The MTK extension enables integration with:
-# #    - PV system models
-# #    - Solar thermal collectors
-# #    - Building energy models (EnergyPlus-like)
-# #    - Daylighting calculations
-# #    - Solar tracking systems
-# #    - Agricultural/greenhouse models
-
-# println("\n✓ ModelingToolkit integration examples completed!")
-# println("  The extension enables composable solar energy system modeling.")
-# println("  See the documentation for more advanced usage patterns.")
diff --git a/src/SolarPosition.jl b/src/SolarPosition.jl
index 6ec5d55..52f04fe 100644
--- a/src/SolarPosition.jl
+++ b/src/SolarPosition.jl
@@ -80,8 +80,9 @@ The simulation time `t` (accessed via `t_nounits`) is in **seconds** from the re
 
 ```julia
 using SolarPosition, ModelingToolkit
-using ModelingToolkit: t_nounits as t
+using ModelingToolkit: t_nounits as t, @named, mtkcompile
 using Dates
+using OrdinaryDiffEq: ODEProblem, solve
 
 @named sun = SolarPositionBlock()
 obs = Observer(51.5, -0.18, 15.0)