r.sun: Incorrect beam radiation for some steep inclined surfaces

raster/r.sun/main.c

@@ -1690,6 +1690,8 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     double locTimeOffset;
     double latitude, longitude;
     double coslat = 0.0;
+    bool shouldBeBestAM = false;
+    bool isBestAM = false;
 
     struct SunGeometryConstDay sunGeom;
     struct SunGeometryVarDay sunVarGeom;
@@ -1817,7 +1819,8 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
     latitude, longitude, sin_phi_l, latid_l, sin_u, cos_u, sin_v, cos_v, lum,  \
     gridGeom, elevin, aspin, slopein, civiltime, linkein, albedo, latin,       \
     coefbh, coefdh, incidout, longin, horizon, beam_rad, insol_time, diff_rad, \
-    refl_rad, glob_rad, mapset, per, decimals, str_step)
+    refl_rad, glob_rad, mapset, per, decimals, str_step, shouldBeBestAM,       \
+    isBestAM)
         {
 #pragma omp for schedule(dynamic) firstprivate(sunGeom, sunVarGeom,      \
                                                sunSlopeGeom, sunRadVar)  \
@@ -1937,8 +1940,21 @@ void calculate(double singleSlope, double singleAspect, double singleAlbedo,
 
                     q1 = gridGeom.sinlat * cos_u * sin_v +
                          gridGeom.coslat * sin_u;
-                    tan_lam_l = -cos_u * cos_v / q1;
-                    sunSlopeGeom.longit_l = atan(tan_lam_l);
+
+                    if (q1 != 0.0) {
+                        tan_lam_l = -cos_u * cos_v / q1;
+                        sunSlopeGeom.longit_l = atan(tan_lam_l);
+                        isBestAM = (tan_lam_l > 0);
+                    }
+                    else {
+                        sunSlopeGeom.longit_l = pihalf;
+                        isBestAM = true;
+                    }
+
+                    shouldBeBestAM = (0.0 < sunSlopeGeom.aspect &&
+                                      sunSlopeGeom.aspect <= M_PI);
+                    sunSlopeGeom.shift12hrs = (shouldBeBestAM != isBestAM);
+
                     sunSlopeGeom.lum_C31_l = cos(latid_l) * sunGeom.cosdecl;
                     sunSlopeGeom.lum_C33_l = sin_phi_l * sunGeom.sindecl;
 

raster/r.sun/rsunlib.c

@@ -326,6 +326,11 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
     /* func = cube; */
     sunVarGeom->isShadow = 0;
 
+    double timeOffset = 0;
+
+    if (sunSlopeGeom->shift12hrs)
+        timeOffset = M_PI;
+
     if (useShadow()) {
         length = 0;
 
@@ -365,7 +370,8 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
                    }
                  */
                 s = sunSlopeGeom->lum_C31_l *
-                        cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) +
+                        cos(-sungeom->timeAngle - sunSlopeGeom->longit_l +
+                            timeOffset) +
                     sunSlopeGeom->lum_C33_l; /* Jenco */
             }
 
@@ -390,7 +396,8 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
                    }
                  */
                 s = sunSlopeGeom->lum_C31_l *
-                        cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) +
+                        cos(-sungeom->timeAngle - sunSlopeGeom->longit_l +
+                            timeOffset) +
                     sunSlopeGeom->lum_C33_l; /* Jenco */
             }
         }
@@ -405,7 +412,7 @@ double lumcline2(struct SunGeometryConstDay *sungeom,
            }
          */
         s = sunSlopeGeom->lum_C31_l *
-                cos(-sungeom->timeAngle - sunSlopeGeom->longit_l) +
+                cos(-sungeom->timeAngle - sunSlopeGeom->longit_l + timeOffset) +
             sunSlopeGeom->lum_C33_l; /* Jenco */
     }
 

raster/r.sun/sunradstruct.h

@@ -65,6 +65,7 @@ struct SunGeometryVarSlope {
     double lum_C33_l;
     double slope;
     double aspect;
+    bool shift12hrs;
 };
 
 struct SolarRadVar {

raster/r.sun/testsuite/test_rsun.py

@@ -293,5 +293,77 @@ def test_run_outputs(self):
         )
 
 
+class TestRSunHighNorthernSlope(TestCase):
+    elevation = "elevation_high_slope"
+    slope = "slope_high_slope"
+    aspect = "aspect_high_slope"
+    beam_radiation = "beam_high_slope"
+
+    @classmethod
+    def setUpClass(cls):
+        """Use temporary region settings"""
+        cls.use_temp_region()
+        cls.runModule("g.region", raster="elevation", res=50, flags="a")
+        cls.runModule(
+            "r.mapcalc",
+            expression=(
+                f"{cls.elevation} = "
+                "-0.0002 * (x() - 637500) * (x() - 637500) +"
+                "0.0002 * (y() - 221750) * (y() - 221750)"
+            ),
+        )
+        cls.runModule(
+            "r.slope.aspect",
+            elevation=cls.elevation,
+            slope=cls.slope,
+            aspect=cls.aspect,
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region"""
+        cls.del_temp_region()
+        cls.runModule(
+            "g.remove",
+            type="raster",
+            flags="f",
+            name=[cls.elevation, cls.aspect, cls.slope],
+        )
+
+    def tearDown(self):
+        self.runModule("g.remove", type="raster", flags="f", name=self.beam_radiation)
+
+    def test_beam_radiation_no_terrain(self):
+        self.assertModule(
+            "r.sun",
+            flags="p",
+            elevation=self.elevation,
+            aspect=self.aspect,
+            slope=self.slope,
+            beam_rad=self.beam_radiation,
+            day=90,
+        )
+        self.assertRasterExists(name=self.beam_radiation)
+        values = "min=50.77406\nmax=6845.00146\nmean=3452.85688\nstddev=1893.56956"
+        self.assertRasterFitsUnivar(
+            raster=self.beam_radiation, reference=values, precision=1e-5
+        )
+
+    def test_beam_radiation_with_terrain(self):
+        self.assertModule(
+            "r.sun",
+            elevation=self.elevation,
+            aspect=self.aspect,
+            slope=self.slope,
+            beam_rad=self.beam_radiation,
+            day=90,
+        )
+        self.assertRasterExists(name=self.beam_radiation)
+        values = "min=44.05826\nmax=6845.00146\nmean=3425.53184\nstddev=1891.98052"
+        self.assertRasterFitsUnivar(
+            raster=self.beam_radiation, reference=values, precision=1e-5
+        )
+
+
 if __name__ == "__main__":
     test()