Adding function sun_position_el (ported from SatelliteToolbox.jl)

Project.toml

@@ -10,12 +10,14 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SatelliteToolboxBase = "9e17983a-0463-41a7-9a16-1682db6d8b66"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 
 [compat]
 Reexport = "1"
 SatelliteToolboxBase = "0.2.1"
 StaticArrays = "1"
 julia = "1.6"
+SatelliteToolboxTransformations = "0.1.3"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/SatelliteToolboxCelestialBodies.jl

@@ -5,6 +5,7 @@ using StaticArrays
 
 @reexport using Dates
 @reexport using SatelliteToolboxBase
+@reexport using SatelliteToolboxTransformations
 
 ############################################################################################
 #                                         Includes

src/sun.jl

@@ -19,9 +19,12 @@
 #
 #   [4] The Astronomical Almanac for the year 2006.
 #
+#   [5] Blanco, Manuel Jesus, Milidonis, Kypros, Bonanos, Aristides. Updating the PSA sun
+#       position algorithm. Solar Energy, vol.212, Elsevier BV,2020-12.
+#
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 
-export sun_position_mod, sun_velocity_mod
+export sun_position_mod, sun_velocity_mod, sun_position_el
 
 """
     sun_position_mod(jd_tdb::Number) -> SVector{3, Float64}
@@ -59,10 +62,10 @@ function sun_position_mod(jd_tdb::Number)
     sin_2Ms = 2 * sin_Ms * cos_Ms
     cos_2Ms = cos_Ms * cos_Ms - sin_Ms * sin_Ms
 
-    # Mean longitude of the Sun [deg].
+    # Mean long of the Sun [deg].
     λ_m = 280.460 + 36_000.771t_tdb
 
-    # Ecliptic latitude of the Sun [deg].
+    # Ecliptic lat of the Sun [deg].
     λ_e = λ_m + 1.914_666_471sin_Ms + 0.019_994_643sin_2Ms
 
     # Obliquity of the ecliptic [deg].
@@ -133,10 +136,10 @@ function sun_velocity_mod(jd_tdb::Number)
     sin_2Ms = 2 * sin_Ms * cos_Ms
     cos_2Ms = cos_Ms * cos_Ms - sin_Ms * sin_Ms
 
-    # Mean longitude of the Sun [deg].
+    # Mean long of the Sun [deg].
     λ_m = 280.460 + 36_000.771t_tdb
 
-    # Ecliptic latitude of the Sun [deg].
+    # Ecliptic lat of the Sun [deg].
     λ_e = λ_m + 1.914_666_471sin_Ms + 0.019_994_643sin_2Ms
 
     # Obliquity of the ecliptic [deg].
@@ -173,3 +176,103 @@ function sun_velocity_mod(jd_tdb::Number)
 
     return vsun_mod
 end
+
+"""
+    sun_position_el(
+        jd::Number,
+        lat::Real=0.0,
+        long::Real=0.0,
+        flag::Char='l',
+    )
+
+Compute the Sun position represented in the Local Horizon Reference System at the
+Julian Day `jd`, lat `lat`, long `long`, Atmospheric pressure `Pressure`,
+Ambient Temperature `Temperature`, Output Flag `flag` and Algorithm `algorithm`. It utilises
+the PSA+ algorithm, (\\[5] (accessed on 2023-07-09)), to compute the topocentric sun position.
+
+# Inputs
+
+    - jd: Julian Day;
+    - lat: Latitude of the observer, in degrees, WSG84;
+    - long: Longitude of the observer, in degrees, WSG84;
+
+# Outputs
+
+## Equatorial system:
+    - α, Right Ascension in degrees;
+    - δ, Declination in degrees;
+
+## Local Coordinates:
+    - ω, Hour angle in degrees;
+    - θ, Zenith in degrees;
+    - γ, Azimuth in degrees;
+
+## Sun vector in (East, North, Zenith):
+    - SunVec, [Nx3] sun vector in (east, north, zenith);
+
+# References
+
+- **[5]**: Blanco, Manuel Jesus, Milidonis, Kypros, Bonanos, Aristides. Updating the PSA sun
+            position algorithm. Solar Energy, vol.212, Elsevier BV,2020-12.
+
+"""
+function sun_position_el(
+    jd::Real,
+    lat::Real=0.0,
+    long::Real=0.0,
+)
+    # Get time data from Julian Date `jd`
+    elapsedJD = jd - JD_J2000
+    _, _, _, Hour, Minute, Second = jd_to_date(jd)
+    DecimalHours = Hour + Minute/60.0 + Second/3600.0
+
+    # PSA+ Algorithm
+    ## Ecliptic Coordinates
+    Ω =  2.267127827e+00 - 9.300339267e-04*elapsedJD #
+    ML = 4.895036035e+00 + 1.720279602e-02*elapsedJD # Mean long
+    MA = 6.239468336e+00 + 1.720200135e-02*elapsedJD # Mean Anomaly
+
+    # Ecliptic long
+    λ₀ = (
+        ML + 3.338320972e-02*sin( MA )
+        + 3.497596876e-04 * sin( 2*MA ) - 1.544353226e-04
+        - 8.689729360e-06*sin( Ω )
+    )
+    # Ecliptic Obliquity
+    ϵ₀ = 4.090904909e-01 - 6.213605399e-09*elapsedJD + 4.418094944e-05*cos(Ω)
+
+    ## Celestial coordinates
+    # Right ascension & declination
+    dY1 = cos( ϵ₀ ) .* sin( λ₀ )
+    dX1 = cos( λ₀ )
+    α = mod2pi(atan( dY1, dX1))
+    δ = asin( sin( ϵ₀ ) .* sin( λ₀ ) )
+
+    ## Topocentric coordinates
+    # Greenwich & Local sidereal time
+    GMST = 6.697096103e+00 + 6.570984737e-02*elapsedJD + DecimalHours
+    LMST = deg2rad( GMST*15 + long )
+    # Hour angle
+    ω = mod2pi(LMST - α)
+
+    # Local coordinates
+    # Zenith
+    θ = acos(cosd(lat)*cos( ω ).*cos( δ ) + sin( δ )*sind(lat))
+    dY2 = -sin(ω)
+    dX2 = tan( δ) * cosd( lat ) - sind(lat)*cos(ω)
+    # Azimuth
+    γ = mod2pi(atan(dY2, dX2))
+
+    # Parallax correction
+    # ToDo: Add radius of earth in Base constants? Then import and replace it below?
+    θ += 6371.01e+03/ASTRONOMICAL_UNIT * sin(θ)
+
+    # East North Zenith Frame
+    SunVec = [
+        sin(γ).*sin(θ),
+        cos(γ).*sin(θ),
+        cos(θ)
+    ]
+
+    return (rad2deg(α), rad2deg(δ), rad2deg(ω), rad2deg(θ), rad2deg(γ), SunVec)
+end

test/sun.jl

@@ -13,6 +13,9 @@
 #   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications. 4th ed.
 #       Microcosm Press, Hawthorn, CA, USA.
 #
+#   [2] Blanco, Manuel Jesus, Milidonis, Kypros, Bonanos, Aristides. Updating the PSA sun
+#       position algorithm. Solar Energy, vol.212, Elsevier BV,2020-12.
+#
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 
 # File: ./src/sun.jl
@@ -71,3 +74,46 @@ end
         @test norm(v - v_n) / norm(v) * 100 < 0.055
     end
 end
+
+############################################################################################
+#                                       Test Results
+############################################################################################
+#
+# Finding the local sun position in local frame using the following data:
+#   JD = 2.458985282701389e6
+#   Latitude = 37.1 °N
+#   Longitude = -2.36 °W
+#
+# Calling the function as:
+#   sun_position_el(2.458985282701389e6, 37.1, -2.36)
+#
+# Must result in :
+#   RA ~ 53.04443267321162 °
+#   DEC ~ 19.108142467612954 °
+#   HA ~ 100.31967018417757 °
+#   ZEN ~ 86.42495381753338 °
+#   AZM ~ 291.337579668121 °
+#   SUN ~ [-0.9296, 0.3632, 0.0624]
+#
+# Accuracy: 50 arcsecs ~ 0.0138889 ° [2, p. 341]
+#
+############################################################################################
+
+@testset "Sun Position Local" begin
+    JD = 2.458985282701389e6
+    Latitude = 37.1
+    Longitude = -2.36
+
+    # Tolerances
+    tol_ang = 1.38889e-2
+    tol_vec = 1e-2
+
+    (ra, dec, ha, zen, az, sun) = sun_position_el(JD, Latitude, Longitude)
+
+    @test ra ≈ 53.04443267321162 atol = tol_ang
+    @test dec ≈ 19.108142467612954 atol = tol_ang
+    @test ha ≈ 100.31967018417757 atol = tol_ang
+    @test zen ≈ 86.42495381753338 atol = tol_ang
+    @test az ≈ 291.337579668121 atol = tol_ang
+    @test sun ≈ [-0.9296, 0.3632, 0.0624] atol = tol_vec
+end