I made some tweaks to r.horizon/main.c to almost half its runtime and
remove a bug. By my tests, it's output is identical to before, but a 219
minute run has been reduced to 110 minutes. Mostly due to removal of
unneeded floor() calls. I could use some people to test and make sure I'm
not introducing any bugs, though.
Thanks,
Seth
--- main.orig 2011-08-14 06:32:49.000000000 -0600
+++ main.c 2012-02-12 19:20:15.000000000 -0700
@@ -144,12 +144,15 @@
/* use G_distance() instead ??!?! */
double distance(double x1, double x2, double y1, double y2)
{
+ double diffX = x1 - x2;
+ double diffY = y1 - y2;
+
if (ll_correction) {
- return DEGREEINMETERS * sqrt(coslatsq * (x1 - x2) * (x1 - x2)
- + (y1 - y2) * (y1 - y2));
+ return DEGREEINMETERS * sqrt(coslatsq * diffX * diffX
+ + diffY * diffY);
}
else {
- return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+ return sqrt(diffX * diffX + diffY * diffY);
}
}
@@ -841,36 +844,28 @@
int new_point()
{
- int iold, jold;
- int succes = 1, succes2 = 1;
- double sx, sy;
- double dx, dy;
-
- iold = ip;
- jold = jp;
+ int iold = ip;
+ int jold = jp;
- while (succes) {
+ while (1) {
yy0 += stepsinangle;
xx0 += stepcosangle;
/* offset 0.5 cell size to get the right cell i, j */
- sx = xx0 * invstepx + offsetx;
- sy = yy0 * invstepy + offsety;
- ip = (int)sx;
- jp = (int)sy;
+ ip = (int) (xx0 * invstepx + offsetx);
+ jp = (int) (yy0 * invstepy + offsety);
/* test outside of raster */
if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m))
return (3);
if ((ip != iold) || (jp != jold)) {
- dx = (double)ip *stepx;
- dy = (double)jp *stepy;
+ double dx = (double)ip *stepx;
+ double dy = (double)jp *stepy;
length = distance(xg0, dx, yg0, dy); /* dist from orig. grid point
to the current grid point */
- succes2 = test_low_res();
- if (succes2 == 1) {
+ if (test_low_res() == 1) {
zp = z[jp][ip];
return (1);
}
@@ -882,54 +877,44 @@
int test_low_res()
{
- int iold100, jold100;
- double sx, sy;
- int delx, dely, mindel;
- double zp100, z2, curvature_diff;
-
- iold100 = ip100;
- jold100 = jp100;
- ip100 = floor(ip / 100.);
- jp100 = floor(jp / 100.);
+ int iold100 = ip100;
+ int jold100 = jp100;
+ ip100 = ip * .01;
+ jp100 = jp * .01;
/*test the new position with low resolution */
if ((ip100 != iold100) || (jp100 != jold100)) {
/* G_debug(3,"%d %d %d %d\n",ip,jp, iold100,jold100); */
/* replace with approximate version
curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
- */
+
+ Code folded into the 'if' statement:
curvature_diff = 0.5 * length * length * invEarth;
z2 = z_orig + curvature_diff + length * tanh0;
zp100 = z100[jp100][ip100];
+ if (zp100 <= z2) {...
+ */
/*G_debug(3,"%d %d %lf %lf \n",ip,jp,z2,zp100); */
- if (zp100 <= z2)
+ if (z100[jp100][ip100] <= z_orig + 0.5 * length * length * invEarth +
length * tanh0)
/*skip to the next lowres cell */
{
- delx = 32000;
- dely = 32000;
+ int mindel;
+ int delx = 32000;
+ int dely = 32000;
+ double sx = (xx0 * invstepx + offsetx) * .01;
+ double sy = (yy0 * invstepy + offsety) * .01;
+
if (cosangle > 0.) {
- sx = xx0 * invstepx + offsetx;
- delx =
- floor(fabs
- ((ceil(sx / 100.) - (sx / 100.)) * distcosangle));
+ delx = fabs((ceil(sx) - sx) * distcosangle);
}
- if (cosangle < 0.) {
- sx = xx0 * invstepx + offsetx;
- delx =
- floor(fabs
- ((floor(sx / 100.) - (sx / 100.)) * distcosangle));
+ else if (cosangle < 0.) {
+ delx = fabs((floor(sx) - sx) * distcosangle);
}
if (sinangle > 0.) {
- sy = yy0 * invstepy + offsety;
- dely =
- floor(fabs
- ((ceil(sy / 100.) - (sy / 100.)) * distsinangle));
+ dely = fabs((ceil(sy) - sy) * distsinangle);
}
else if (sinangle < 0.) {
- sy = yy0 * invstepy + offsety;
- dely =
- floor(fabs
- ((floor(jp / 100.) - (sy / 100.)) * distsinangle));
+ dely = fabs((floor(sy) - sy) * distsinangle);
}
mindel = min(delx, dely);
@@ -953,17 +938,14 @@
double searching()
{
- double z2;
- double curvature_diff;
- int succes = 1;
-
if (zp == UNDEFZ)
return 0;
while (1) {
- succes = new_point();
+ double z2;
+ double curvature_diff;
- if (succes != 1) {
+ if (new_point() != 1) {
break;
}