r.horizon bugfix and speedup

[sprice]

I mentioned this on the mailing list, but didn't get a response. So here's a more official report:

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;

}