Merge pull request #29 from statnett/fix-hour-angle

Correct hour angle for longitude

docs/api/equations/solar_angles.rst

@@ -0,0 +1,5 @@
+Solar angles
+-------------
+
+.. automodule:: linerate.equations.solar_angles
+    :members:

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 = {

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(

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)

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():