r.horizon: fix distance output (#3589)

OSGeo · Apr 18, 2024 · 69fa8c1 · 69fa8c1
1 parent bb27c05
commit 69fa8c1
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 25 deletions.
diff --git a/raster/r.horizon/main.c b/raster/r.horizon/main.c
@@ -87,6 +87,7 @@ typedef struct {
     int ip, jp;
     int ip100, jp100;
     double zp;
+    double length;
 } SearchPoint;
 
 typedef struct {
@@ -969,7 +970,7 @@ int new_point(const Geometry *geometry, const OriginPoint *origin_point,
             double dx = (double)search_point->ip * geometry->stepx;
             double dy = (double)search_point->jp * geometry->stepy;
 
-            horizon->length =
+            search_point->length =
                 distance(origin_point->xg0, dx, origin_point->yg0, dy,
                          origin_point->coslatsq); /* dist from orig. grid point
                                               to the current grid point */
@@ -1000,9 +1001,9 @@ int test_low_res(const Geometry *geometry, const OriginPoint *origin_point,
            curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
          */
         double curvature_diff =
-            0.5 * horizon->length * horizon->length * invEarth;
+            0.5 * search_point->length * search_point->length * invEarth;
         double z2 = origin_point->z_orig + curvature_diff +
-                    horizon->length * horizon->tanh0;
+                    search_point->length * horizon->tanh0;
         double zp100 = z100[search_point->jp100][search_point->ip100];
         G_debug(2, "ip:%d jp:%d z2:%lf zp100:%lf \n", search_point->ip,
                 search_point->jp, z2, zp100);
@@ -1073,6 +1074,7 @@ HorizonProperties horizon_height(const Geometry *geometry,
     search_point.zp = origin_point->z_orig;
     search_point.ip100 = floor(origin_point->xg0 * geometry->invstepx / 100.);
     search_point.jp100 = floor(origin_point->yg0 * geometry->invstepy / 100.);
+    search_point.length = 0;
 
     horizon.length = 0;
     horizon.tanh0 = 0;
@@ -1092,22 +1094,23 @@ HorizonProperties horizon_height(const Geometry *geometry,
 
         /* curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS)); */
         double curvature_diff =
-            0.5 * horizon.length * horizon.length * invEarth;
+            0.5 * search_point.length * search_point.length * invEarth;
 
         double z2 = origin_point->z_orig + curvature_diff +
-                    horizon.length * horizon.tanh0;
+                    search_point.length * horizon.tanh0;
 
         if (z2 < search_point.zp) {
             horizon.tanh0 =
                 (search_point.zp - origin_point->z_orig - curvature_diff) /
-                horizon.length;
+                search_point.length;
+            horizon.length = search_point.length;
         }
 
         if (z2 >= geometry->zmax) {
             break;
         }
 
-        if (horizon.length >= origin_point->maxlength) {
+        if (search_point.length >= origin_point->maxlength) {
             break;
         }
     }

diff --git a/raster/r.horizon/testsuite/test_r_horizon.py b/raster/r.horizon/testsuite/test_r_horizon.py
@@ -87,24 +87,24 @@
 """
 
 ref5 = """azimuth,horizon_height,horizon_distance
-0.000000,0.197017,5010.039920
-20.000000,0.196832,5017.668781
-40.000000,0.196875,5017.818251
-60.000000,0.196689,5017.220346
-80.000000,0.196847,5014.299552
-100.000000,0.196645,5019.531851
-120.000000,0.196969,5014.957627
-140.000000,0.196778,5020.328674
-160.000000,0.196863,5013.431958
-180.000000,0.197017,5010.039920
-200.000000,0.196832,5014.229751
-220.000000,0.196875,5011.387034
-240.000000,0.196689,5017.220346
-260.000000,0.196847,5014.299552
-280.000000,0.196645,5019.531851
-300.000000,0.196969,5014.957627
-320.000000,0.196778,5020.328674
-340.000000,0.196863,5013.431958
+0.000000,0.197017,5000.040000
+20.000000,0.196832,5004.837660
+40.000000,0.196875,5003.728610
+60.000000,0.196689,5008.552685
+80.000000,0.196847,5004.448022
+100.000000,0.196645,5009.690609
+120.000000,0.196969,5001.279836
+140.000000,0.196778,5006.246099
+160.000000,0.196863,5004.018385
+180.000000,0.197017,5000.040000
+200.000000,0.196832,5004.837660
+220.000000,0.196875,5003.728610
+240.000000,0.196689,5008.552685
+260.000000,0.196847,5004.448022
+280.000000,0.196645,5009.690609
+300.000000,0.196969,5001.279836
+320.000000,0.196778,5006.246099
+340.000000,0.196863,5004.018385
 """