Enhance solar position calculations by adding standard time and fractional hour functions, and update related algorithms and tests

Author: langestefan

src/Basic/Basic.jl

@@ -13,11 +13,12 @@ export BasicAlgorithm
 using PVSimBase.GeoLocation: Location
 using DocStringExtensions
 using TimeZones
-using DynamicQuantities: Quantity
+using DynamicQuantities
 using Interfaces
 
 using SolarPosition.PositionInterface: SolarPositionAlgorithm, SolarPositionInterface,
-                                       declination, equation_of_time, offset_hours
+                                       declination, equation_of_time, offset_hours,
+                                       standard_time, fractional_hour
 import ..PositionInterface
 
 """
@@ -28,18 +29,47 @@ struct BasicAlgorithm <: SolarPositionAlgorithm
 end
 
 function PositionInterface.sunpos(algorithm::BasicAlgorithm, timestamp::ZonedDateTime)
-    longitude = algorithm.location.longitude
-    latitude = algorithm.location.latitude
 
     # declination [deg]
     β = declination(timestamp)
 
-    # Local standard Time Merdidian [deg]
-    LSTM = 15 * offset_hours(timestamp)
+    # Longitude of the standard Time Merdidian [deg]
+    λ_LSTM = 15 * offset_hours(timestamp)us"deg"
 
-    # equation of time 
-    E_qt = equation_of_time(timestamp)
-    return (Quantity(0.0, deg = 1), Quantity(0.0, deg = 1))
+    # Local longitude  / latitude [deg]
+    λ_local = algorithm.location.longitude
+    ϕ_local = algorithm.location.latitude
+
+    # Local time, without daylights saving time [h]
+    T_local = fractional_hour(standard_time(timestamp)) * us"h"
+
+    # equation of time from [deg] to [h], 360 deg = 24 h
+    E_qt = equation_of_time(timestamp) * 24us"h" / 360us"deg"
+    println("Timestamp: ", timestamp)
+    println("E_qt: ", E_qt)
+    println("T_local: ", T_local)
+    println("λ_LSTM: ", λ_LSTM)
+    println("λ_local: ", λ_local)
+    println("ϕ_local: ", ϕ_local)
+    println("β: ", β)
+
+    # solar time [h]
+    T_solar = T_local + E_qt + (λ_LSTM - λ_local) / 15us"deg/h"
+    println("T_solar: ", T_solar)
+
+    # hour angle [deg]
+    H = 15us"deg/h" * (12us"h" - T_solar)
+    println("H: ", H)
+
+    # solar elevation angle [deg]
+    α = asind(sind(β.value) * sind(ϕ_local.value) +
+              cosd(β.value) * cosd(ϕ_local.value) * cosd(H.value))
+
+    # solar azimuth angle [deg]
+    A = acosd((sind(β.value) * cosd(ϕ_local.value) -
+               cosd(β.value) * sind(ϕ_local.value) * cosd(H.value)) / cosd(α))
+    return (
+        Quantity(α, SymbolicDimensions, deg = 1), Quantity(A, SymbolicDimensions, deg = 1))
 end
 
 algorithm = BasicAlgorithm(Location(latitude = 0.0, longitude = 0.0))

src/PositionInterface/PositionInterface.jl

@@ -1,7 +1,7 @@
 module PositionInterface
 
 using Interfaces
-using DynamicQuantities: Quantity
+using DynamicQuantities
 using Dates
 using PVSimBase
 using TimeZones

src/PositionInterface/angles.jl

@@ -15,7 +15,7 @@ The output is a `Quantity` with the unit of degrees.
 function declination(d::Number)
     (0 <= d <= 365) || throw(ArgumentError("day of the year must be in the range [0, 365]"))
     δ = -23.45cosd((360 / 365) * (d + 10))
-    return Quantity(δ, deg = 1)
+    return Quantity(δ, SymbolicDimensions, deg = 1)
 end
 
 """

src/PositionInterface/time.jl

@@ -1,4 +1,4 @@
-export equation_of_time, offset_hours
+export equation_of_time, offset_hours, standard_time, fractional_hour
 
 """
     equation_of_time(t::Union{Number, DateTime})
@@ -36,11 +36,11 @@ function equation_of_time(jd::Number)
     # Ecliptic latitude of the Sun.
     λ_ecliptic = mod(λ_m + 1.914666471sin_Ms + 0.019994643sin_2Ms, 360) |> deg2rad
 
-    # Compute the equation of time [rad].
+    # Compute the equation of time [deg].
     eot = -1.914666471sin_Ms - 0.019994643sin_2Ms + 2.466sin(2λ_ecliptic) -
           0.0053sin(4λ_ecliptic)
 
-    return Quantity(eot, deg = 1)
+    return Quantity(eot, SymbolicDimensions, deg = 1)
 end
 
 """
@@ -57,4 +57,24 @@ equation_of_time(t::ZonedDateTime) = equation_of_time(DateTime(t))
 Return the offset in hours of the timezone `tz` with respect to UTC-0.
 """
 offset_hours(tz::FixedTimeZone) = Hour(tz.offset.std).value
-offset_hours(t::ZonedDateTime) = offset_hours(FixedTimeZone(t))
+offset_hours(t::ZonedDateTime) = offset_hours(FixedTimeZone(t))
+
+"""
+    standard_time(t::ZonedDateTime)
+
+Get the standard time of a given `ZonedDateTime` object.
+"""
+standard_time(t::ZonedDateTime) = t.utc_datetime + t.zone.offset.std
+
+"""
+    fractional_hour(t::DateTime)
+
+Get the current time as a fraction x/24 of the day.
+"""
+# fractional_hour(t::ZonedDateTime) = hour(t) + minute(t) / 60 + second(t) / 3600
+# dt = t - floor(t, Hour(24))
+function fractional_hour(t::DateTime)
+    dt = t - floor(t, Hour(24))
+    return Float64(dt.value) / (60 * 60 * 1000)
+end
+fractional_hour(t::ZonedDateTime) = fractional_hour(DateTime(t))

src/SolarPositionAlgorithms/SolarPositionAlgorithms.jl

@@ -9,7 +9,7 @@
 """
 
 using Dates: DateTime
-using DynamicQuantities: Quantity
+using DynamicQuantities
 using SolarPosition.PositionInterface
 
 @testset "Function declination" begin

test/test-angles.jl

@@ -1,16 +1,33 @@
-using TimeZones: now, TimeZone
+using TimeZones
+using Dates
 using Interfaces
 
 using PVSimBase.GeoLocation: Location
 using SolarPosition.Basic
 using SolarPosition.PositionInterface
 
-@testset "Basic algorithm" begin
-    test_location = Location(latitude = 0.0, longitude = 0.0)
-    alg = BasicAlgorithm(test_location)
-    @test alg isa BasicAlgorithm
-    @test alg.location == test_location
-    pos = sunpos(alg, now(TimeZone("UTC")))
+# @testset "Basic algorithm" begin
+#     test_location = Location(latitude = 0.0, longitude = 0.0)
+#     alg = BasicAlgorithm(test_location)
+#     @test alg isa BasicAlgorithm
+#     @test alg.location == test_location
+#     pos = sunpos(alg, now(TimeZone("UTC")))
+# end
+
+@testset "Algorithms" begin
+    # latitude, longitude, timestamp, elevation, azimuth
+    cases = [
+        (39.743, -105.178,
+        ZonedDateTime(2020, 2, 10, 0, 0, TimeZone("UTC-7")), -39.44us"deg", 75.84us"deg")
+    ]
+
+    for (lat, lon, ts, el, az) in cases
+        loc = Location(latitude = lat, longitude = lon)
+        alg = BasicAlgorithm(loc)
+        pos = sunpos(alg, ts)
+        @test pos[1]≈el atol=1.0us"deg"
+        @test pos[2]≈az atol=1.0us"deg"
+    end
 end
 
 @testset "PositionInterfaces" begin

test/test-interfaces.jl

@@ -9,7 +9,7 @@
 """
 
 using Dates: DateTime
-using DynamicQuantities: Quantity
+using DynamicQuantities
 using SolarPosition.PositionInterface
 
 @testset "Function equation_of_time" begin