r.horizon: fix bug in point mode computation (#3288)

In point mode, some variables were not set properly, leading to major errors. Whenever it entered https://github.com/OSGeo/grass/blob/main/raster/r.horizon/main.c#L940 condition, it resulted in an error due to global variables sinangle and cosangle equal to 0. Raster mode was fine, those variables are initialized to proper values.

The fix uses com_par function (which was previously unused) with minor modification of its input. It sets those variables (and some others) and can be used in both modes, so now the variables in both modes should be set consistently.

The test computes horizon for an artificial, circular area (circle with z value 0 and 1000 around), so the horizon angles should be close to identical. This test would have previously resulted in zero horizon height everywhere.

raster/r.horizon/main.c

@@ -87,7 +87,7 @@ double amax1(double, double);
 double amin1(double, double);
 int min(int, int);
 int max(int, int);
-void com_par(double angle);
+void com_par(void);
 int is_shadow(void);
 double horizon_height(void);
 void calculate_shadow(void);
@@ -718,13 +718,11 @@ int max(int arg1, int arg2)
 
 /**********************************************************/
 
-void com_par(double angle)
+void com_par()
 {
-    sinangle = sin(angle);
     if (fabs(sinangle) < 0.0000001) {
         sinangle = 0.;
     }
-    cosangle = cos(angle);
     if (fabs(cosangle) < 0.0000001) {
         cosangle = 0.;
     }
@@ -836,8 +834,9 @@ void calculate_shadow(void)
 
         delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
 
-        stepsinangle = stepxy * delt_nor / delt_dist;
-        stepcosangle = stepxy * delt_east / delt_dist;
+        sinangle = delt_nor / delt_dist;
+        cosangle = delt_east / delt_dist;
+        com_par();
 
         shadow_angle = horizon_height();
 
@@ -1190,25 +1189,8 @@ void calculate(double xcoord, double ycoord, int buffer_e, int buffer_w,
                         sqrt(delt_east * delt_east + delt_nor * delt_nor);
 
                     sinangle = delt_nor / delt_dist;
-                    if (fabs(sinangle) < 0.0000001) {
-                        sinangle = 0.;
-                    }
                     cosangle = delt_east / delt_dist;
-                    if (fabs(cosangle) < 0.0000001) {
-                        cosangle = 0.;
-                    }
-                    distsinangle = 32000;
-                    distcosangle = 32000;
-
-                    if (sinangle != 0.) {
-                        distsinangle = 100. / (distxy * sinangle);
-                    }
-                    if (cosangle != 0.) {
-                        distcosangle = 100. / (distxy * cosangle);
-                    }
-
-                    stepsinangle = stepxy * sinangle;
-                    stepcosangle = stepxy * cosangle;
+                    com_par();
 
                     z_orig = zp = z[j][i];
                     maxlength = (zmax - z_orig) / TANMINANGLE;

raster/r.horizon/testsuite/test_r_horizon.py

@@ -63,18 +63,50 @@
 160.000000,0.015356
 """
 
+ref4 = """azimuth,horizon_height
+0.000000,0.197017
+20.000000,0.196832
+40.000000,0.196875
+60.000000,0.196689
+80.000000,0.196847
+100.000000,0.196645
+120.000000,0.196969
+140.000000,0.196778
+160.000000,0.196863
+180.000000,0.197017
+200.000000,0.196832
+220.000000,0.196875
+240.000000,0.196689
+260.000000,0.196847
+280.000000,0.196645
+300.000000,0.196969
+320.000000,0.196778
+340.000000,0.196863
+"""
+
 
 class TestHorizon(TestCase):
+    circle = "circle"
     horizon = "test_horizon_from_elevation"
     horizon_output = "test_horizon_output_from_elevation"
 
     @classmethod
     def setUpClass(cls):
         cls.use_temp_region()
         cls.runModule("g.region", raster="elevation")
+        cls.runModule(
+            "r.circle",
+            flags="b",
+            output=cls.circle,
+            coordinates=(637505, 221755),
+            min=5000,
+            multiplier=1000,
+        )
+        cls.runModule("r.null", map=cls.circle, null=0)
 
     @classmethod
     def tearDownClass(cls):
+        cls.runModule("g.remove", flags="f", type="raster", name=cls.circle)
         cls.del_temp_region()
 
     def tearDown(self):
@@ -112,6 +144,20 @@ def test_point_mode_multiple_direction(self):
         stdout = module.outputs.stdout
         self.assertMultiLineEqual(first=ref2, second=stdout)
 
+    def test_point_mode_multiple_direction_artificial(self):
+        """Test mode with 1 point and multiple directions with artificial surface"""
+        module = SimpleModule(
+            "r.horizon",
+            elevation=self.circle,
+            coordinates=(637505, 221755),
+            output=self.horizon,
+            direction=0,
+            step=20,
+        )
+        self.assertModule(module)
+        stdout = module.outputs.stdout
+        self.assertMultiLineEqual(first=ref4, second=stdout)
+
     def test_raster_mode_one_direction(self):
         """Test mode with 1 point and one direction"""
         module = SimpleModule(