From edf3837cf34206fa1f4cf3165e1ae0d14192af38 Mon Sep 17 00:00:00 2001 From: Halvor Lund Date: Thu, 30 May 2024 14:35:23 +0200 Subject: [PATCH] Correct hour angle for longitude Previously, hour angle was calculated at longitude 0, which was incorrect. --- examples/plot_solar_heating_comparison.py | 2 +- linerate/equations/solar_angles.py | 14 ++++++++++---- linerate/model.py | 4 ++-- tests/equations/test_solar_angles.py | 18 +++++++++++++----- 4 files changed, 26 insertions(+), 12 deletions(-) diff --git a/examples/plot_solar_heating_comparison.py b/examples/plot_solar_heating_comparison.py index 8d78650..3680781 100644 --- a/examples/plot_solar_heating_comparison.py +++ b/examples/plot_solar_heating_comparison.py @@ -21,7 +21,7 @@ # ^^^^^^^^^^^^^^^^^^^^^ time = np.datetime64("2022-06-01T13:00") -omega = solar_angles.compute_hour_angle_relative_to_noon(time) +omega = solar_angles.compute_hour_angle_relative_to_noon(time, 0) delta = solar_angles.compute_solar_declination(time) vals_with_range = { diff --git a/linerate/equations/solar_angles.py b/linerate/equations/solar_angles.py index 63b64a6..15ec729 100644 --- a/linerate/equations/solar_angles.py +++ b/linerate/equations/solar_angles.py @@ -24,7 +24,7 @@ def _get_minute_of_hour(when: Date) -> Unitless: return (when.astype(MinuteResolutionType) - when.astype(HourResolutionType)).astype(float) -def compute_hour_angle_relative_to_noon(when: Date) -> Radian: +def compute_hour_angle_relative_to_noon(when: Date, longitude: Degrees) -> Radian: r"""Compute the hour angle. Described in the text on p. 18 of :cite:p:`ieee738`. The hour angle is the number of hours @@ -35,6 +35,9 @@ def compute_hour_angle_relative_to_noon(when: Date) -> Radian: which is the same as 15^\circ. The hour angle is used when calculating the solar altitude. + This function does not take into account the difference between apparent/actual + and mean solar time, which means that the result may be up to 15 minutes from the + correct hour angle. Parameters ---------- @@ -46,10 +49,13 @@ def compute_hour_angle_relative_to_noon(when: Date) -> Radian: Union[float, float64, ndarray[Any, dtype[float64]]] :math:'\omega~\left[\text{radian}\right]`. The hour angle relative to noon. """ - hour = _get_hour_of_day(when) - minute = _get_minute_of_hour(when) + utc_hour = _get_hour_of_day(when) + utc_minute = _get_minute_of_hour(when) pi = np.pi - return (-12 + hour + minute / 60) * (pi / 12) # pi/12 is 15 degrees + # We add longitude/15 since 15 degrees of longitude increases solar hour by 1 + return np.mod((-12 + utc_hour + utc_minute / 60 + longitude / 15), 24) * ( + pi / 12 + ) # pi/12 is 15 degrees def compute_solar_declination( diff --git a/linerate/model.py b/linerate/model.py index 29bcf31..3977c18 100644 --- a/linerate/model.py +++ b/linerate/model.py @@ -339,7 +339,7 @@ def compute_solar_heating( N_s = self.weather.clearness_ratio D = self.span.conductor.conductor_diameter - omega = solar_angles.compute_hour_angle_relative_to_noon(self.time) + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) delta = solar_angles.compute_solar_declination(self.time) sin_H_s = solar_angles.compute_sin_solar_altitude(phi, delta, omega) chi = solar_angles.compute_solar_azimuth_variable(phi, delta, omega) @@ -491,7 +491,7 @@ def compute_solar_heating( y = self.span.conductor_altitude # H_e in IEEE D = self.span.conductor.conductor_diameter # D_0 in IEEE - omega = solar_angles.compute_hour_angle_relative_to_noon(self.time) + omega = solar_angles.compute_hour_angle_relative_to_noon(self.time, self.span.longitude) delta = solar_angles.compute_solar_declination(self.time) sin_H_c = solar_angles.compute_sin_solar_altitude(phi, delta, omega) Q_s = ieee738.solar_heating.compute_total_heat_flux_density(sin_H_c, True) diff --git a/tests/equations/test_solar_angles.py b/tests/equations/test_solar_angles.py index 92e7b3b..89c2a05 100644 --- a/tests/equations/test_solar_angles.py +++ b/tests/equations/test_solar_angles.py @@ -41,7 +41,14 @@ def test_get_minute_of_hour_with_example(): def test_hour_angle_relative_to_noon_with_example(): when = np.datetime64("2022-06-01T12:00") omega = 0 - assert solar_angles.compute_hour_angle_relative_to_noon(when) == approx(omega) + assert solar_angles.compute_hour_angle_relative_to_noon(when, 0) == approx(omega) + + +def test_hour_angle_relative_to_noon_norway(): + when = np.datetime64("2022-12-01T12:00+01:00") + omega = 0 + longitude = 15 # Longitude at timezone +01:00 + assert solar_angles.compute_hour_angle_relative_to_noon(when, longitude) == approx(omega) def test_solar_azimuth_variable_with_example(): @@ -119,12 +126,13 @@ def test_solar_declination_scales_correctly_with_day_of_year(day): when=st.datetimes( min_value=datetime.datetime(2022, 1, 1, 1, 0), max_value=datetime.datetime(2022, 12, 31, 23, 0), - ) + ), + longitude=st.floats(min_value=-180, max_value=180), ) -def test_solar_declination_scales_with_dates_and_times(when): - omega = (-12 + when.hour + when.minute / 60) * np.pi / 12 +def test_solar_declination_scales_with_dates_and_times(when, longitude): + omega = ((-12 + when.hour + when.minute / 60 + longitude / 15) % 24) * np.pi / 12 when = np.datetime64(when) - assert omega == approx(solar_angles.compute_hour_angle_relative_to_noon(when)) + assert omega == approx(solar_angles.compute_hour_angle_relative_to_noon(when, longitude)) def test_solar_declination_with_examples():