| 10 | | |
| 11 | | --- main.orig 2011-08-14 06:32:49.000000000 -0600 |
| 12 | | +++ main.c 2012-02-12 19:20:15.000000000 -0700 |
| 13 | | @@ -144,12 +144,15 @@ |
| 14 | | /* use G_distance() instead ??!?! */ |
| 15 | | double distance(double x1, double x2, double y1, double y2) |
| 16 | | { |
| 17 | | + double diffX = x1 - x2; |
| 18 | | + double diffY = y1 - y2; |
| 19 | | + |
| 20 | | if (ll_correction) { |
| 21 | | - return DEGREEINMETERS * sqrt(coslatsq * (x1 - x2) * (x1 - x2) |
| 22 | | - + (y1 - y2) * (y1 - y2)); |
| 23 | | + return DEGREEINMETERS * sqrt(coslatsq * diffX * diffX |
| 24 | | + + diffY * diffY); |
| 25 | | } |
| 26 | | else { |
| 27 | | - return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); |
| 28 | | + return sqrt(diffX * diffX + diffY * diffY); |
| 29 | | } |
| 30 | | } |
| 31 | | |
| 32 | | @@ -841,36 +844,28 @@ |
| 33 | | |
| 34 | | int new_point() |
| 35 | | { |
| 36 | | - int iold, jold; |
| 37 | | - int succes = 1, succes2 = 1; |
| 38 | | - double sx, sy; |
| 39 | | - double dx, dy; |
| 40 | | - |
| 41 | | - iold = ip; |
| 42 | | - jold = jp; |
| 43 | | + int iold = ip; |
| 44 | | + int jold = jp; |
| 45 | | |
| 46 | | - while (succes) { |
| 47 | | + while (1) { |
| 48 | | yy0 += stepsinangle; |
| 49 | | xx0 += stepcosangle; |
| 50 | | |
| 51 | | |
| 52 | | /* offset 0.5 cell size to get the right cell i, j */ |
| 53 | | - sx = xx0 * invstepx + offsetx; |
| 54 | | - sy = yy0 * invstepy + offsety; |
| 55 | | - ip = (int)sx; |
| 56 | | - jp = (int)sy; |
| 57 | | + ip = (int) (xx0 * invstepx + offsetx); |
| 58 | | + jp = (int) (yy0 * invstepy + offsety); |
| 59 | | |
| 60 | | /* test outside of raster */ |
| 61 | | if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m)) |
| 62 | | return (3); |
| 63 | | |
| 64 | | if ((ip != iold) || (jp != jold)) { |
| 65 | | - dx = (double)ip *stepx; |
| 66 | | - dy = (double)jp *stepy; |
| 67 | | + double dx = (double)ip *stepx; |
| 68 | | + double dy = (double)jp *stepy; |
| 69 | | |
| 70 | | length = distance(xg0, dx, yg0, dy); /* dist from orig. grid point |
| 71 | | to the current grid point */ |
| 72 | | - succes2 = test_low_res(); |
| 73 | | - if (succes2 == 1) { |
| 74 | | + if (test_low_res() == 1) { |
| 75 | | zp = z[jp][ip]; |
| 76 | | return (1); |
| 77 | | } |
| 78 | | @@ -882,54 +877,44 @@ |
| 79 | | |
| 80 | | int test_low_res() |
| 81 | | { |
| 82 | | - int iold100, jold100; |
| 83 | | - double sx, sy; |
| 84 | | - int delx, dely, mindel; |
| 85 | | - double zp100, z2, curvature_diff; |
| 86 | | - |
| 87 | | - iold100 = ip100; |
| 88 | | - jold100 = jp100; |
| 89 | | - ip100 = floor(ip / 100.); |
| 90 | | - jp100 = floor(jp / 100.); |
| 91 | | + int iold100 = ip100; |
| 92 | | + int jold100 = jp100; |
| 93 | | + ip100 = ip * .01; |
| 94 | | + jp100 = jp * .01; |
| 95 | | /*test the new position with low resolution */ |
| 96 | | if ((ip100 != iold100) || (jp100 != jold100)) { |
| 97 | | /* G_debug(3,"%d %d %d %d\n",ip,jp, iold100,jold100); */ |
| 98 | | /* replace with approximate version |
| 99 | | curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS)); |
| 100 | | - */ |
| 101 | | + |
| 102 | | + Code folded into the 'if' statement: |
| 103 | | curvature_diff = 0.5 * length * length * invEarth; |
| 104 | | z2 = z_orig + curvature_diff + length * tanh0; |
| 105 | | zp100 = z100[jp100][ip100]; |
| 106 | | + if (zp100 <= z2) {... |
| 107 | | + */ |
| 108 | | /*G_debug(3,"%d %d %lf %lf \n",ip,jp,z2,zp100); */ |
| 109 | | |
| 110 | | - if (zp100 <= z2) |
| 111 | | + if (z100[jp100][ip100] <= z_orig + 0.5 * length * length * invEarth + |
| 112 | | length * tanh0) |
| 113 | | /*skip to the next lowres cell */ |
| 114 | | { |
| 115 | | - delx = 32000; |
| 116 | | - dely = 32000; |
| 117 | | + int mindel; |
| 118 | | + int delx = 32000; |
| 119 | | + int dely = 32000; |
| 120 | | + double sx = (xx0 * invstepx + offsetx) * .01; |
| 121 | | + double sy = (yy0 * invstepy + offsety) * .01; |
| 122 | | + |
| 123 | | if (cosangle > 0.) { |
| 124 | | - sx = xx0 * invstepx + offsetx; |
| 125 | | - delx = |
| 126 | | - floor(fabs |
| 127 | | - ((ceil(sx / 100.) - (sx / 100.)) * distcosangle)); |
| 128 | | + delx = fabs((ceil(sx) - sx) * distcosangle); |
| 129 | | } |
| 130 | | - if (cosangle < 0.) { |
| 131 | | - sx = xx0 * invstepx + offsetx; |
| 132 | | - delx = |
| 133 | | - floor(fabs |
| 134 | | - ((floor(sx / 100.) - (sx / 100.)) * distcosangle)); |
| 135 | | + else if (cosangle < 0.) { |
| 136 | | + delx = fabs((floor(sx) - sx) * distcosangle); |
| 137 | | } |
| 138 | | if (sinangle > 0.) { |
| 139 | | - sy = yy0 * invstepy + offsety; |
| 140 | | - dely = |
| 141 | | - floor(fabs |
| 142 | | - ((ceil(sy / 100.) - (sy / 100.)) * distsinangle)); |
| 143 | | + dely = fabs((ceil(sy) - sy) * distsinangle); |
| 144 | | } |
| 145 | | else if (sinangle < 0.) { |
| 146 | | - sy = yy0 * invstepy + offsety; |
| 147 | | - dely = |
| 148 | | - floor(fabs |
| 149 | | - ((floor(jp / 100.) - (sy / 100.)) * distsinangle)); |
| 150 | | + dely = fabs((floor(sy) - sy) * distsinangle); |
| 151 | | } |
| 152 | | |
| 153 | | mindel = min(delx, dely); |
| 154 | | @@ -953,17 +938,14 @@ |
| 155 | | |
| 156 | | double searching() |
| 157 | | { |
| 158 | | - double z2; |
| 159 | | - double curvature_diff; |
| 160 | | - int succes = 1; |
| 161 | | - |
| 162 | | if (zp == UNDEFZ) |
| 163 | | return 0; |
| 164 | | |
| 165 | | while (1) { |
| 166 | | - succes = new_point(); |
| 167 | | + double z2; |
| 168 | | + double curvature_diff; |
| 169 | | |
| 170 | | - if (succes != 1) { |
| 171 | | + if (new_point() != 1) { |
| 172 | | break; |
| 173 | | } |
