Lines 97.82% 90 / 92
Functions and Methods 88.88% 16 / 18
Classes and Traits 0.00% 0 / 1
Name Lines Functions and Methods CRAP Classes and Traits
NoaaSolarCalculator 97.82% 90 / 92 88.88% 16 / 18 30 0.00% 0 / 1
 findEventOn 75.00% 3 / 4 0.00% 0 / 1 2.06
 findSolarNoonOn 100.00% 5 / 5 100.00% 1 / 1 2
 findCrossingOn 100.00% 11 / 11 100.00% 1 / 1 4
 computePositionAt 100.00% 5 / 5 100.00% 1 / 1 1
 isSunUpAt 100.00% 3 / 3 100.00% 1 / 1 1
 computeGeometricPositionAt 100.00% 15 / 15 100.00% 1 / 1 2
 findUtcMidnightOfTransit 100.00% 4 / 4 100.00% 1 / 1 1
 computeTransitAt 100.00% 2 / 2 100.00% 1 / 1 1
 computeHourAngleDegrees 100.00% 4 / 4 100.00% 1 / 1 2
 computeCoordinatesAt 100.00% 21 / 21 100.00% 1 / 1 1
 computeGeometricElevationOf 100.00% 2 / 2 100.00% 1 / 1 1
 computeRefractionDegrees 85.71% 6 / 7 0.00% 0 / 1 5.07
 computeHorizonDipDegrees 100.00% 3 / 3 100.00% 1 / 1 2
 convertToJulianDay 100.00% 1 / 1 100.00% 1 / 1 1
 convertToInstant 100.00% 2 / 2 100.00% 1 / 1 1
 clampLatitude 100.00% 1 / 1 100.00% 1 / 1 1
 clampUnit 100.00% 1 / 1 100.00% 1 / 1 1
 wrap 100.00% 1 / 1 100.00% 1 / 1 1
1<?php
2
3declare(strict_types=1);
4
5namespace Stewart\Sun;
6
7use DateTimeImmutable;
8use Stewart\Contracts\Sun\GeoLocation;
9use Stewart\Contracts\Sun\SunDirection;
10use Stewart\Contracts\Sun\SunElevationCrossing;
11use Stewart\Contracts\Sun\SunEvent;
12use Stewart\Contracts\Sun\SunPosition;
13use Stewart\Contracts\Time\Instant;
14
15final readonly class NoaaSolarCalculator
16{
17    private const float JULIAN_DAY_AT_UNIX_EPOCH = 2440587.5;
18    private const float JULIAN_DAY_AT_J2000 = 2451545.0;
19    private const float DAYS_PER_JULIAN_CENTURY = 36525.0;
20    private const float SECONDS_PER_DAY = 86400.0;
21    private const float MINUTES_PER_DAY = 1440.0;
22    private const float MINUTES_PER_LONGITUDE_DEGREE = 4.0;
23    private const float NOON_MINUTES = 720.0;
24    private const float EARTH_RADIUS_METERS = 6356900.0;
25    // Hour angles are undefined at the poles themselves.
26    private const float LATITUDE_LIMIT_DEGREES = 89.8;
27    private const int REFINEMENT_PASSES = 3;
28
29    public function findEventOn(DateTimeImmutable $localDate, GeoLocation $location, SunEvent $event): ?Instant
30    {
31        $crossing = $event->findElevationCrossing();
32
33        if ($crossing === null) {
34            return $this->findSolarNoonOn($localDate, $location);
35        }
36
37        return $this->findCrossingOn($localDate, $location, $crossing);
38    }
39
40    public function findSolarNoonOn(DateTimeImmutable $localDate, GeoLocation $location): Instant
41    {
42        $utcMidnight = $this->findUtcMidnightOfTransit($localDate, $location);
43        $transit = $utcMidnight + 0.5;
44
45        for ($pass = 0; $pass < self::REFINEMENT_PASSES; ++$pass) {
46            $transit = $this->computeTransitAt($utcMidnight, $location, $this->computeCoordinatesAt($transit));
47        }
48
49        return $this->convertToInstant($transit);
50    }
51
52    public function findCrossingOn(DateTimeImmutable $localDate, GeoLocation $location, SunElevationCrossing $crossing): ?Instant
53    {
54        $utcMidnight = $this->findUtcMidnightOfTransit($localDate, $location);
55        $targetElevation = $this->computeGeometricElevationOf($crossing, $location);
56        $sign = $crossing->direction === SunDirection::Rising ? -1.0 : 1.0;
57        $estimate = $utcMidnight + 0.5;
58
59        for ($pass = 0; $pass < self::REFINEMENT_PASSES; ++$pass) {
60            $coordinates = $this->computeCoordinatesAt($estimate);
61            $hourAngle = $this->computeHourAngleDegrees($location, $coordinates, $targetElevation);
62
63            if ($hourAngle === null) {
64                return null;
65            }
66
67            $estimate = $this->computeTransitAt($utcMidnight, $location, $coordinates) + $sign * $hourAngle / 360.0;
68        }
69
70        return $this->convertToInstant($estimate);
71    }
72
73    public function computePositionAt(Instant $moment, GeoLocation $location): SunPosition
74    {
75        $geometric = $this->computeGeometricPositionAt($moment, $location);
76
77        return new SunPosition(
78            $geometric->azimuthDegrees,
79            $geometric->elevationDegrees + $this->computeRefractionDegrees($geometric->elevationDegrees),
80        );
81    }
82
83    public function isSunUpAt(Instant $moment, GeoLocation $location): bool
84    {
85        $sunrise = SunEvent::Sunrise->findElevationCrossing();
86        \assert($sunrise !== null);
87
88        return $this->computeGeometricPositionAt($moment, $location)->elevationDegrees > $this->computeGeometricElevationOf($sunrise, $location);
89    }
90
91    private function computeGeometricPositionAt(Instant $moment, GeoLocation $location): SunPosition
92    {
93        $julianDay = $this->convertToJulianDay($moment);
94        $coordinates = $this->computeCoordinatesAt($julianDay);
95        $latitude = deg2rad($this->clampLatitude($location->latitude));
96        $declination = deg2rad($coordinates->declinationDegrees);
97
98        $utcMinutes = ($julianDay + 0.5 - floor($julianDay + 0.5)) * self::MINUTES_PER_DAY;
99        $trueSolarMinutes = $this->wrap($utcMinutes + $coordinates->equationOfTimeMinutes + self::MINUTES_PER_LONGITUDE_DEGREE * $location->longitude, self::MINUTES_PER_DAY);
100        $hourAngle = $trueSolarMinutes / self::MINUTES_PER_LONGITUDE_DEGREE - 180.0;
101
102        $cosZenith = sin($latitude) * sin($declination) + cos($latitude) * cos($declination) * cos(deg2rad($hourAngle));
103        $zenith = acos($this->clampUnit($cosZenith));
104        $cosAzimuth = (sin($latitude) * cos($zenith) - sin($declination)) / (cos($latitude) * sin($zenith));
105        $azimuth = rad2deg(acos($this->clampUnit($cosAzimuth)));
106
107        return new SunPosition(
108            $hourAngle > 0 ? $this->wrap($azimuth + 180.0, 360.0) : $this->wrap(540.0 - $azimuth, 360.0),
109            90.0 - rad2deg($zenith),
110        );
111    }
112
113    private function findUtcMidnightOfTransit(DateTimeImmutable $localDate, GeoLocation $location): float
114    {
115        $localNoon = $this->convertToJulianDay(Instant::fromDateTime($localDate->setTime(12, 0)));
116        $utcMidnight = floor($localNoon + 0.5) - 0.5;
117        $transit = $this->computeTransitAt($utcMidnight, $location, $this->computeCoordinatesAt($localNoon));
118
119        return $utcMidnight + round($localNoon - $transit);
120    }
121
122    private function computeTransitAt(float $utcMidnight, GeoLocation $location, SolarCoordinates $coordinates): float
123    {
124        $minutes = self::NOON_MINUTES - self::MINUTES_PER_LONGITUDE_DEGREE * $location->longitude - $coordinates->equationOfTimeMinutes;
125
126        return $utcMidnight + $minutes / self::MINUTES_PER_DAY;
127    }
128
129    private function computeHourAngleDegrees(GeoLocation $location, SolarCoordinates $coordinates, float $elevationDegrees): ?float
130    {
131        $latitude = deg2rad($this->clampLatitude($location->latitude));
132        $declination = deg2rad($coordinates->declinationDegrees);
133
134        $cosHourAngle = (sin(deg2rad($elevationDegrees)) - sin($latitude) * sin($declination)) / (cos($latitude) * cos($declination));
135
136        return abs($cosHourAngle) > 1.0 ? null : rad2deg(acos($cosHourAngle));
137    }
138
139    private function computeCoordinatesAt(float $julianDay): SolarCoordinates
140    {
141        $century = ($julianDay - self::JULIAN_DAY_AT_J2000) / self::DAYS_PER_JULIAN_CENTURY;
142
143        $meanLongitude = deg2rad($this->wrap(280.46646 + $century * (36000.76983 + $century * 0.0003032), 360.0));
144        $meanAnomaly = deg2rad(357.52911 + $century * (35999.05029 - 0.0001537 * $century));
145        $eccentricity = 0.016708634 - $century * (0.000042037 + 0.0000001267 * $century);
146
147        $equationOfCenter = sin($meanAnomaly) * (1.914602 - $century * (0.004817 + 0.000014 * $century))
148            + sin(2 * $meanAnomaly) * (0.019993 - 0.000101 * $century)
149            + sin(3 * $meanAnomaly) * 0.000289;
150        $nutation = deg2rad(125.04 - 1934.136 * $century);
151        $apparentLongitude = deg2rad(rad2deg($meanLongitude) + $equationOfCenter - 0.00569 - 0.00478 * sin($nutation));
152
153        $meanObliquity = 23 + (26 + (21.448 - $century * (46.815 + $century * (0.00059 - $century * 0.001813))) / 60) / 60;
154        $obliquity = deg2rad($meanObliquity + 0.00256 * cos($nutation));
155
156        $obliquityFactor = tan($obliquity / 2) ** 2;
157        $equationOfTime = $obliquityFactor * sin(2 * $meanLongitude)
158            - 2 * $eccentricity * sin($meanAnomaly)
159            + 4 * $eccentricity * $obliquityFactor * sin($meanAnomaly) * cos(2 * $meanLongitude)
160            - 0.5 * $obliquityFactor * $obliquityFactor * sin(4 * $meanLongitude)
161            - 1.25 * $eccentricity * $eccentricity * sin(2 * $meanAnomaly);
162
163        return new SolarCoordinates(
164            rad2deg(asin(sin($obliquity) * sin($apparentLongitude))),
165            self::MINUTES_PER_LONGITUDE_DEGREE * rad2deg($equationOfTime),
166        );
167    }
168
169    private function computeGeometricElevationOf(SunElevationCrossing $crossing, GeoLocation $location): float
170    {
171        $seenElevation = $crossing->elevationDegrees - $this->computeHorizonDipDegrees($location);
172
173        return $seenElevation - $this->computeRefractionDegrees($seenElevation);
174    }
175
176    // Astral's piecewise refraction model, so times agree with Home Assistant's sun integration.
177    private function computeRefractionDegrees(float $elevationDegrees): float
178    {
179        if ($elevationDegrees >= 85.0) {
180            return 0.0;
181        }
182
183        $tangent = tan(deg2rad($elevationDegrees));
184
185        $arcSeconds = match (true) {
186            $elevationDegrees > 5.0 => 58.1 / $tangent - 0.07 / $tangent ** 3 + 0.000086 / $tangent ** 5,
187            $elevationDegrees > -0.575 => 1735.0 + $elevationDegrees * (-518.2 + $elevationDegrees * (103.4 + $elevationDegrees * (-12.79 + $elevationDegrees * 0.711))),
188            default => -20.774 / $tangent,
189        };
190
191        return $arcSeconds / 3600.0;
192    }
193
194    private function computeHorizonDipDegrees(GeoLocation $location): float
195    {
196        if ($location->elevationMeters <= 0.0) {
197            return 0.0;
198        }
199
200        return rad2deg(acos(self::EARTH_RADIUS_METERS / (self::EARTH_RADIUS_METERS + $location->elevationMeters)));
201    }
202
203    private function convertToJulianDay(Instant $moment): float
204    {
205        return $moment->toEpochMicroseconds() / 1_000_000 / self::SECONDS_PER_DAY + self::JULIAN_DAY_AT_UNIX_EPOCH;
206    }
207
208    private function convertToInstant(float $julianDay): Instant
209    {
210        $epochSeconds = (int) round(($julianDay - self::JULIAN_DAY_AT_UNIX_EPOCH) * self::SECONDS_PER_DAY);
211
212        return Instant::fromEpochMicroseconds($epochSeconds * 1_000_000);
213    }
214
215    private function clampLatitude(float $latitude): float
216    {
217        return max(-self::LATITUDE_LIMIT_DEGREES, min(self::LATITUDE_LIMIT_DEGREES, $latitude));
218    }
219
220    private function clampUnit(float $value): float
221    {
222        return max(-1.0, min(1.0, $value));
223    }
224
225    private function wrap(float $value, float $period): float
226    {
227        return $value - $period * floor($value / $period);
228    }
229}