Changes between Initial Version and Version 2 of Ticket #1575


Ignore:
Timestamp:
Feb 14, 2012, 3:02:22 PM (8 years ago)
Author:
hamish
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #1575

    • Property Keywords r.horizon added
    • Property Platform changed from Unspecified to All
    • Property Component changed from Default to Raster
  • Ticket #1575 – Description

    initial v2  
    88Thanks,
    99Seth
    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         }