src/netskel.c

   1 /* netskel.c
   2  * Survex network reduction - remove trailing traverses and concatenate
   3  * traverses between junctions
   4  * Copyright (C) 1991-2004,2005,2006,2010,2011,2012,2013,2014,2015 Olly Betts
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 2 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
  19  */
  20
  21 /* #define BLUNDER_DETECTION */
  22
  23 #if 0
  24 #define DEBUG_INVALID 1
  25 #define VALIDATE 1
  26 #define DUMP_NETWORK 1
  27 #endif
  28
  29 #ifdef HAVE_CONFIG_H
  30 # include <config.h>
  31 #endif
  32
  33 #include "validate.h"
  34 #include "debug.h"
  35 #include "cavern.h"
  36 #include "commands.h"
  37 #include "filename.h"
  38 #include "message.h"
  39 #include "filelist.h"
  40 #include "img_hosted.h"
  41 #include "netartic.h"
  42 #include "netbits.h"
  43 #include "netskel.h"
  44 #include "network.h"
  45 #include "out.h"
  46
  47 #define sqrdd(X) (sqrd((X)[0]) + sqrd((X)[1]) + sqrd((X)[2]))
  48
  49 typedef struct Stack {
  50    struct Link *join1, *join2;
  51    struct Stack *next;
  52 } stack;
  53
  54 typedef struct StackTr {
  55    struct Link *join1;
  56    struct StackTr *next;
  57 } stackTrail;
  58
  59 /* string used between things in traverse printouts eg \1 - \2 - \3 -...*/
  60 static const char *szLink = " - ";
  61 static const char *szLinkEq = " = "; /* use this one for equates */
  62
  63 #if 0
  64 #define fprint_prefix(FH, NAME) BLK((fprint_prefix)((FH), (NAME));\
  65                                     fprintf((FH), " [%p]", (void*)(NAME)); )
  66 #endif
  67
  68 static stack *ptr; /* Ptr to TRaverse linked list for in-between travs */
  69 static stackTrail *ptrTrail; /* Ptr to TRaverse linked list for trail travs*/
  70
  71 static void remove_trailing_travs(void);
  72 static void remove_travs(void);
  73 static void replace_travs(void);
  74 static void replace_trailing_travs(void);
  75 static void write_passage_models(void);
  76
  77 static void concatenate_trav(node *stn, int i);
  78
  79 static void err_stat(int cLegsTrav, double lenTrav,
  80                      double eTot, double eTotTheo,
  81                      double hTot, double hTotTheo,
  82                      double vTot, double vTotTheo);
  83
  84 extern void
  85 solve_network(void /*node *stnlist*/)
  86 {
  87    static int first_solve = 1;
  88    node *stn;
  89
  90    /* We can't average across solving to fix positions. */
  91    clear_last_leg();
  92
  93    if (stnlist == NULL) {
  94       if (first_solve) fatalerror(/*No survey data*/43);
  95       /* We've had a *solve followed by another *solve (or the implicit
  96        * *solve at the end of the data.  Don't moan about that. */
  97       return;
  98    }
  99    ptr = NULL;
 100    ptrTrail = NULL;
 101    dump_network();
 102
 103    if (first_solve && !pcs->proj_str && !proj_str_out) {
 104       /* If we haven't already solved to find some station positions, and
 105        * there's no specified coordinate system, then check if there are any
 106        * fixed points, and if there aren't, invent one at (0,0,0).
 107        *
 108        * We do this first so the solving part is just like the standard case -
 109        * this avoid problems, such as sub-nodes of the invented fix having been
 110        * removed.  It also means we can fix the "first" station, which makes
 111        * more sense to the user. */
 112       FOR_EACH_STN(stn, stnlist)
 113          if (fixed(stn)) break;
 114
 115       if (!stn) {
 116         /* If we've had a *solve and all the new survey since then is hanging,
 117          * we don't want to invent a fixed point.  We want to complain but
 118          * the easiest way to is just to continue processing and let
 119          * articulate() catch this condition as it will any hanging survey
 120          * data. */
 121          node *stnFirst = NULL;
 122
 123          /* New stations are pushed onto the head of the list, so the
 124           * first station added is the last in the list. */
 125          FOR_EACH_STN(stn, stnlist) {
 126              /* Prefer a station with legs attached when choosing one to fix
 127               * so that if there's a hanging station on a nosurvey leg we pick
 128               * the main clump of survey data. */
 129              if (stnFirst && !stnFirst->leg[0]) continue;
 130              stnFirst = stn;
 131          }
 132
 133          /* If we've got nosurvey legs, then the station we find to fix could
 134           * have no real legs attached. */
 135          SVX_ASSERT2(nosurveyhead || stnFirst->leg[0],
 136                      "no fixed stns, but we've got a zero node!");
 137          SVX_ASSERT2(stnFirst, "no stations left in net!");
 138          stn = stnFirst;
 139          printf(msg(/*Survey has no fixed points. Therefore I’ve fixed %s at (0,0,0)*/72),
 140                 sprint_prefix(stn->name));
 141          putnl();
 142          POS(stn,0) = (real)0.0;
 143          POS(stn,1) = (real)0.0;
 144          POS(stn,2) = (real)0.0;
 145          fix(stn);
 146       }
 147    }
 148
 149    first_solve = 0;
 150
 151    remove_trailing_travs();
 152    validate(); dump_network();
 153    remove_travs();
 154    validate(); dump_network();
 155    remove_subnets();
 156    validate(); dump_network();
 157    articulate();
 158    validate(); dump_network();
 159    replace_subnets();
 160    validate(); dump_network();
 161    replace_travs();
 162    validate(); dump_network();
 163    replace_trailing_travs();
 164    validate(); dump_network();
 165
 166    /* Now write out any passage models. */
 167    write_passage_models();
 168 }
 169
 170 static void
 171 remove_trailing_travs(void)
 172 {
 173    node *stn;
 174    /* TRANSLATORS: In French, Eric chose to use the terminology used by
 175     * toporobot: "sequence" for the English "traverse", which makes sense
 176     * (although toporobot actually uses this term to mean something more
 177     * specific).  Feel free to follow this lead if you can't think of a better
 178     * term - these messages mostly indicate how processing is progressing.
 179     *
 180     * A trailing traverse is a dead end back to a junction. */
 181    out_current_action(msg(/*Removing trailing traverses*/125));
 182    FOR_EACH_STN(stn, stnlist) {
 183       if (!fixed(stn) && one_node(stn)) {
 184          int i = 0;
 185          int j;
 186          node *stn2 = stn;
 187          stackTrail *trav;
 188
 189 #if PRINT_NETBITS
 190          printf("Removed trailing trav ");
 191 #endif
 192          do {
 193             struct Link *leg;
 194 #if PRINT_NETBITS
 195             print_prefix(stn2->name); printf("<%p>",stn2); fputs(szLink, stdout);
 196 #endif
 197             remove_stn_from_list(&stnlist, stn2);
 198             leg = stn2->leg[i];
 199             j = reverse_leg_dirn(leg);
 200             stn2 = leg->l.to;
 201             i = j ^ 1; /* flip direction for other leg of 2 node */
 202             /* stop if fixed or 3 or 1 node */
 203          } while (two_node(stn2) && !fixed(stn2));
 204
 205          /* put traverse on stack */
 206          trav = osnew(stackTrail);
 207          trav->join1 = stn2->leg[j];
 208          trav->next = ptrTrail;
 209          ptrTrail = trav;
 210
 211          /* We want to keep all 2-nodes using legs 0 and 1 and all one nodes
 212           * using leg 0 so we may need to swap leg j with leg 2 (for a 3 node)
 213           * or leg 1 (for a fixed 2 node) */
 214          if ((j == 0 && !one_node(stn2)) || (j == 1 && three_node(stn2))) {
 215             /* i is the direction to swap with */
 216             i = (three_node(stn2)) ? 2 : 1;
 217             /* change the other direction of leg i to use leg j */
 218             reverse_leg(stn2->leg[i])->l.reverse += j - i;
 219             stn2->leg[j] = stn2->leg[i];
 220             j = i;
 221          }
 222          stn2->leg[j] = NULL;
 223
 224 #if PRINT_NETBITS
 225          print_prefix(stn2->name); printf("<%p>",stn2); putnl();
 226 #endif
 227       }
 228    }
 229 }
 230
 231 static void
 232 remove_travs(void)
 233 {
 234    node *stn;
 235    /* TRANSLATORS: In French, Eric chose to use the terminology used by
 236     * toporobot: "sequence" for the English "traverse", which makes sense
 237     * (although toporobot actually uses this term to mean something more
 238     * specific).  Feel free to follow this lead if you can't think of a better
 239     * term - these messages mostly indicate how processing is progressing. */
 240    out_current_action(msg(/*Concatenating traverses*/126));
 241    FOR_EACH_STN(stn, stnlist) {
 242       if (fixed(stn) || three_node(stn)) {
 243          int d;
 244          for (d = 0; d <= 2; d++) {
 245             linkfor *leg = stn->leg[d];
 246             if (leg && !(leg->l.reverse & FLAG_REPLACEMENTLEG))
 247                concatenate_trav(stn, d);
 248          }
 249       }
 250    }
 251 }
 252
 253 static void
 254 concatenate_trav(node *stn, int i)
 255 {
 256    int j;
 257    stack *trav;
 258    node *stn2;
 259    linkfor *newleg, *newleg2;
 260
 261    stn2 = stn->leg[i]->l.to;
 262    /* Reject single legs as they may be already concatenated traverses */
 263    if (fixed(stn2) || !two_node(stn2)) return;
 264
 265    trav = osnew(stack);
 266    newleg2 = (linkfor*)osnew(linkrev);
 267
 268 #if PRINT_NETBITS
 269    printf("Concatenating trav "); print_prefix(stn->name); printf("<%p>",stn);
 270 #endif
 271
 272    newleg2->l.to = stn;
 273    newleg2->l.reverse = i | FLAG_REPLACEMENTLEG;
 274    trav->join1 = stn->leg[i];
 275
 276    j = reverse_leg_dirn(stn->leg[i]);
 277    SVX_ASSERT(j == 0 || j == 1);
 278
 279    newleg = copy_link(stn->leg[i]);
 280
 281    while (1) {
 282       stn = stn2;
 283
 284 #if PRINT_NETBITS
 285       fputs(szLink, stdout); print_prefix(stn->name); printf("<%p>",stn);
 286 #endif
 287
 288       /* stop if fixed or 3 or 1 node */
 289       if (fixed(stn) || !two_node(stn)) break;
 290
 291       remove_stn_from_list(&stnlist, stn);
 292
 293       i = j ^ 1; /* flip direction for other leg of 2 node */
 294
 295       stn2 = stn->leg[i]->l.to;
 296       j = reverse_leg_dirn(stn->leg[i]);
 297
 298       addto_link(newleg, stn->leg[i]);
 299    }
 300
 301    trav->join2 = stn->leg[j];
 302    trav->next = ptr;
 303    ptr = trav;
 304
 305    newleg->l.to = stn;
 306    newleg->l.reverse = j | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
 307
 308    newleg2->l.to->leg[reverse_leg_dirn(newleg2)] = newleg;
 309    /* i.e. stn->leg[i] = newleg; with original stn and i */
 310
 311    stn->leg[j] = newleg2;
 312
 313 #if PRINT_NETBITS
 314    putchar(' ');
 315    print_var(&(newleg->v));
 316    printf("\nStacked ");
 317    print_prefix(newleg2->l.to->name);
 318    printf(",%d-", reverse_leg_dirn(newleg2));
 319    print_prefix(stn->name);
 320    printf(",%d\n", j);
 321 #endif
 322 }
 323
 324 #ifdef BLUNDER_DETECTION
 325 /* expected_error is actually squared... */
 326 /* only called if fhErrStat != NULL */
 327 static void
 328 do_gross(delta e, delta v, node *stn1, node *stn2, double expected_error)
 329 {
 330    double hsqrd, rsqrd, s, cx, cy, cz;
 331    double tot;
 332    int i;
 333    int output = 0;
 334    prefix *name1 = stn1->name, *name2 = stn2->name;
 335
 336 #if 0
 337 printf( "e = ( %.2f, %.2f, %.2f )", e[0], e[1], e[2] );
 338 printf( " v = ( %.2f, %.2f, %.2f )\n", v[0], v[1], v[2] );
 339 #endif
 340    hsqrd = sqrd(v[0]) + sqrd(v[1]);
 341    rsqrd = hsqrd + sqrd(v[2]);
 342    if (rsqrd == 0.0) return;
 343
 344    cx = v[0] + e[0];
 345    cy = v[1] + e[1];
 346    cz = v[2] + e[2];
 347
 348    s = (e[0] * v[0] + e[1] * v[1] + e[2] * v[2]) / rsqrd;
 349    tot = 0;
 350    for (i = 2; i >= 0; i--) tot += sqrd(e[i] - v[i] * s);
 351
 352    if (tot <= expected_error) {
 353       if (!output) {
 354          fprint_prefix(fhErrStat, name1);
 355          fputs("->", fhErrStat);
 356          fprint_prefix(fhErrStat, name2);
 357       }
 358       fprintf(fhErrStat, " L: %.2f", sqrt(tot));
 359       /* checked - works */
 360       fprintf(fhErrStat, " (%.2fm -> %.2fm)", sqrt(sqrdd(v)), sqrt(sqrdd(v)) * (1 - s));
 361       output = 1;
 362    }
 363
 364    s = sqrd(cx) + sqrd(cy);
 365    if (s > 0.0) {
 366       s = hsqrd / s;
 367       SVX_ASSERT(s >= 0.0);
 368       s = sqrt(s);
 369       s = 1 - s;
 370       tot = sqrd(cx * s) + sqrd(cy * s) + sqrd(e[2]);
 371       if (tot <= expected_error) {
 372          double newval, oldval;
 373          if (!output) {
 374             fprint_prefix(fhErrStat, name1);
 375             fputs("->", fhErrStat);
 376             fprint_prefix(fhErrStat, name2);
 377          }
 378          fprintf(fhErrStat, " B: %.2f", sqrt(tot));
 379          /* checked - works */
 380          newval = deg(atan2(cx, cy));
 381          if (newval < 0) newval += 360;
 382          oldval = deg(atan2(v[0], v[1]));
 383          if (oldval < 0) oldval += 360;
 384          fprintf(fhErrStat, " (%.2fdeg -> %.2fdeg)", oldval, newval);
 385          output = 1;
 386       }
 387    }
 388
 389    if (hsqrd > 0.0) {
 390       double nx, ny;
 391       s = (e[0] * v[1] - e[1] * v[0]) / hsqrd;
 392       nx = cx - s * v[1];
 393       ny = cy + s * v[0];
 394       s = sqrd(nx) + sqrd(ny) + sqrd(cz);
 395       if (s > 0.0) {
 396          s = rsqrd / s;
 397          SVX_ASSERT(s >= 0);
 398          s = sqrt(s);
 399          tot = sqrd(cx - s * nx) + sqrd(cy - s * ny) + sqrd(cz - s * cz);
 400          if (tot <= expected_error) {
 401             if (!output) {
 402                fprint_prefix(fhErrStat, name1);
 403                fputs("->", fhErrStat);
 404                fprint_prefix(fhErrStat, name2);
 405             }
 406             fprintf(fhErrStat, " G: %.2f", sqrt(tot));
 407             /* checked - works */
 408             fprintf(fhErrStat, " (%.2fdeg -> %.2fdeg)",
 409                     deg(atan2(v[2], sqrt(v[0] * v[0] + v[1] * v[1]))),
 410                     deg(atan2(cz, sqrt(nx * nx + ny * ny))));
 411             output = 1;
 412          }
 413       }
 414    }
 415    if (output) fputnl(fhErrStat);
 416 }
 417 #endif
 418
 419 static void
 420 replace_travs(void)
 421 {
 422    stack *ptrOld;
 423    node *stn1, *stn2, *stn3;
 424    int i, j, k;
 425    double eTot = 0, lenTrav = 0, lenTot;
 426    double eTotTheo = 0;
 427    double vTot = 0, vTotTheo = 0, hTot = 0, hTotTheo = 0;
 428    delta e, sc;
 429    bool fEquate; /* used to indicate equates in output */
 430    int cLegsTrav = 0;
 431    bool fArtic;
 432
 433     /* TRANSLATORS: In French, Eric chose to use the terminology used by
 434      * toporobot: "sequence" for the English "traverse", which makes sense
 435      * (although toporobot actually uses this term to mean something more
 436      * specific).  Feel free to follow this lead if you can't think of a better
 437      * term - these messages mostly indicate how processing is progressing. */
 438    out_current_action(msg(/*Calculating traverses*/127));
 439
 440    if (!fhErrStat && !fSuppress)
 441       fhErrStat = safe_fopen_with_ext(fnm_output_base, EXT_SVX_ERRS, "w");
 442
 443    if (!pimg) {
 444       char *fnm = add_ext(fnm_output_base, EXT_SVX_3D);
 445       filename_register_output(fnm);
 446       pimg = img_open_write_cs(fnm, s_str(&survey_title), proj_str_out,
 447                                img_FFLAG_SEPARATOR(output_separator));
 448       if (!pimg) fatalerror(img_error(), fnm);
 449       osfree(fnm);
 450    }
 451
 452    /* First do all the one leg traverses */
 453    FOR_EACH_STN(stn1, stnlist) {
 454 #if PRINT_NETBITS
 455       printf("One leg traverses from ");
 456       print_prefix(stn1->name);
 457       printf(" [%p]\n", stn1);
 458 #endif
 459       for (i = 0; i <= 2; i++) {
 460          linkfor *leg = stn1->leg[i];
 461          if (leg && data_here(leg) &&
 462              !(leg->l.reverse & (FLAG_REPLACEMENTLEG | FLAG_FAKE))) {
 463             SVX_ASSERT(fixed(stn1));
 464             SVX_ASSERT(!fZeros(&leg->v));
 465
 466             stn2 = leg->l.to;
 467             if (TSTBIT(leg->l.flags, FLAGS_SURFACE)) {
 468                stn1->name->sflags |= BIT(SFLAGS_SURFACE);
 469                stn2->name->sflags |= BIT(SFLAGS_SURFACE);
 470             } else {
 471                stn1->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 472                stn2->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 473             }
 474             img_write_item(pimg, img_MOVE, 0, NULL,
 475                            POS(stn1, 0), POS(stn1, 1), POS(stn1, 2));
 476             if (leg->meta) {
 477                 pimg->days1 = leg->meta->days1;
 478                 pimg->days2 = leg->meta->days2;
 479             } else {
 480                 pimg->days1 = pimg->days2 = -1;
 481             }
 482             pimg->style = (leg->l.flags >> FLAGS_STYLE_BIT0) & 0x07;
 483             img_write_item(pimg, img_LINE, leg->l.flags & FLAGS_MASK,
 484                            sprint_prefix(stn1->name->up),
 485                            POS(stn2, 0), POS(stn2, 1), POS(stn2, 2));
 486             if (!(leg->l.reverse & FLAG_ARTICULATION)) {
 487 #ifdef BLUNDER_DETECTION
 488                delta err;
 489                int do_blunder;
 490 #else
 491                if (fhErrStat) {
 492                   fprint_prefix(fhErrStat, stn1->name);
 493                   fputs(szLink, fhErrStat);
 494                   fprint_prefix(fhErrStat, stn2->name);
 495                }
 496 #endif
 497                subdd(&e, &POSD(stn2), &POSD(stn1));
 498                subdd(&e, &e, &leg->d);
 499                if (fhErrStat) {
 500                   eTot = sqrdd(e);
 501                   hTot = sqrd(e[0]) + sqrd(e[1]);
 502                   vTot = sqrd(e[2]);
 503 #ifndef NO_COVARIANCES
 504                   /* FIXME: what about covariances? */
 505                   hTotTheo = leg->v[0] + leg->v[1];
 506                   vTotTheo = leg->v[2];
 507                   eTotTheo = hTotTheo + vTotTheo;
 508 #else
 509                   hTotTheo = leg->v[0] + leg->v[1];
 510                   vTotTheo = leg->v[2];
 511                   eTotTheo = hTotTheo + vTotTheo;
 512 #endif
 513 #ifdef BLUNDER_DETECTION
 514                   memcpy(&err, &e, sizeof(delta));
 515                   do_blunder = (eTot > eTotTheo);
 516                   fputs("\ntraverse ", fhErrStat);
 517                   fprint_prefix(fhErrStat, stn1->name);
 518                   fputs("->", fhErrStat);
 519                   fprint_prefix(fhErrStat, stn2->name);
 520                   fprintf(fhErrStat, " e=(%.2f, %.2f, %.2f) mag=%.2f %s\n",
 521                           e[0], e[1], e[2], sqrt(eTot),
 522                           (do_blunder ? "suspect:" : "OK"));
 523                   if (do_blunder)
 524                      do_gross(err, leg->d, stn1, stn2, eTotTheo);
 525 #endif
 526                   err_stat(1, sqrt(sqrdd(leg->d)), eTot, eTotTheo,
 527                            hTot, hTotTheo, vTot, vTotTheo);
 528                }
 529             }
 530          }
 531       }
 532    }
 533
 534    while (ptr != NULL) {
 535       /* work out where traverse should be reconnected */
 536       linkfor *leg = ptr->join1;
 537       leg = reverse_leg(leg);
 538       stn1 = leg->l.to;
 539       i = reverse_leg_dirn(leg);
 540
 541       leg = ptr->join2;
 542       leg = reverse_leg(leg);
 543       stn2 = leg->l.to;
 544       j = reverse_leg_dirn(leg);
 545
 546 #if PRINT_NETBITS
 547       printf(" Trav ");
 548       print_prefix(stn1->name);
 549       printf("<%p>[%d]%s...%s", stn1, i, szLink, szLink);
 550       print_prefix(stn2->name);
 551       printf("<%p>[%d]\n", stn2, j);
 552 #endif
 553
 554       SVX_ASSERT(fixed(stn1));
 555       SVX_ASSERT(fixed(stn2));
 556
 557       /* calculate scaling factors for error distribution */
 558       eTot = 0.0;
 559       hTot = vTot = 0.0;
 560       SVX_ASSERT(data_here(stn1->leg[i]));
 561       if (fZeros(&stn1->leg[i]->v)) {
 562          sc[0] = sc[1] = sc[2] = 0.0;
 563       } else {
 564          subdd(&e, &POSD(stn2), &POSD(stn1));
 565          subdd(&e, &e, &stn1->leg[i]->d);
 566          eTot = sqrdd(e);
 567          hTot = sqrd(e[0]) + sqrd(e[1]);
 568          vTot = sqrd(e[2]);
 569          divds(&sc, &e, &stn1->leg[i]->v);
 570       }
 571 #ifndef NO_COVARIANCES
 572       /* FIXME: what about covariances? */
 573       hTotTheo = stn1->leg[i]->v[0] + stn1->leg[i]->v[1];
 574       vTotTheo = stn1->leg[i]->v[2];
 575 #else
 576       hTotTheo = stn1->leg[i]->v[0] + stn1->leg[i]->v[1];
 577       vTotTheo = stn1->leg[i]->v[2];
 578 #endif
 579       eTotTheo = hTotTheo + vTotTheo;
 580       cLegsTrav = 0;
 581       lenTrav = 0.0;
 582       img_write_item(pimg, img_MOVE, 0, NULL,
 583                      POS(stn1, 0), POS(stn1, 1), POS(stn1, 2));
 584
 585       fArtic = stn1->leg[i]->l.reverse & FLAG_ARTICULATION;
 586       osfree(stn1->leg[i]);
 587       stn1->leg[i] = ptr->join1; /* put old link back in */
 588
 589       osfree(stn2->leg[j]);
 590       stn2->leg[j] = ptr->join2; /* and the other end */
 591
 592 #ifdef BLUNDER_DETECTION
 593       delta err;
 594       int do_blunder;
 595       memcpy(&err, &e, sizeof(delta));
 596       do_blunder = (eTot > eTotTheo);
 597       if (fhErrStat && !fArtic) {
 598          fputs("\ntraverse ", fhErrStat);
 599          fprint_prefix(fhErrStat, stn1->name);
 600          fputs("->", fhErrStat);
 601          fprint_prefix(fhErrStat, stn2->name);
 602          fprintf(fhErrStat, " e=(%.2f, %.2f, %.2f) mag=%.2f %s\n",
 603                  e[0], e[1], e[2], sqrt(eTot),
 604                  (do_blunder ? "suspect:" : "OK"));
 605       }
 606 #endif
 607       while (true) {
 608          int reached_end;
 609          prefix *leg_pfx;
 610
 611          fEquate = true;
 612          /* get next node in traverse
 613           * should have stn3->leg[k]->l.to == stn1 */
 614          stn3 = stn1->leg[i]->l.to;
 615          k = reverse_leg_dirn(stn1->leg[i]);
 616          SVX_ASSERT2(stn3->leg[k]->l.to == stn1,
 617                  "reverse leg doesn't reciprocate");
 618
 619          reached_end = (stn3 == stn2 && k == j);
 620
 621          if (data_here(stn1->leg[i])) {
 622             leg_pfx = stn1->name->up;
 623             leg = stn1->leg[i];
 624 #ifdef BLUNDER_DETECTION
 625             if (do_blunder && fhErrStat)
 626                do_gross(err, leg->d, stn1, stn3, eTotTheo);
 627 #endif
 628             if (!reached_end)
 629                adddd(&POSD(stn3), &POSD(stn1), &leg->d);
 630          } else {
 631             leg_pfx = stn3->name->up;
 632             leg = stn3->leg[k];
 633 #ifdef BLUNDER_DETECTION
 634             if (do_blunder && fhErrStat)
 635                do_gross(err, leg->d, stn1, stn3, eTotTheo);
 636 #endif
 637             if (!reached_end)
 638                subdd(&POSD(stn3), &POSD(stn1), &leg->d);
 639          }
 640
 641          lenTot = sqrdd(leg->d);
 642
 643          if (!fZeros(&leg->v)) fEquate = false;
 644          if (!reached_end) {
 645             add_stn_to_list(&stnlist, stn3);
 646             if (!fEquate) {
 647                mulsd(&e, &leg->v, &sc);
 648                adddd(&POSD(stn3), &POSD(stn3), &e);
 649             }
 650             fix(stn3);
 651          }
 652
 653          if (!(leg->l.reverse & (FLAG_REPLACEMENTLEG | FLAG_FAKE))) {
 654              if (TSTBIT(leg->l.flags, FLAGS_SURFACE)) {
 655                 stn1->name->sflags |= BIT(SFLAGS_SURFACE);
 656                 stn3->name->sflags |= BIT(SFLAGS_SURFACE);
 657              } else {
 658                 stn1->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 659                 stn3->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 660              }
 661
 662             SVX_ASSERT(!fEquate);
 663             SVX_ASSERT(!fZeros(&leg->v));
 664             if (leg->meta) {
 665                 pimg->days1 = leg->meta->days1;
 666                 pimg->days2 = leg->meta->days2;
 667             } else {
 668                 pimg->days1 = pimg->days2 = -1;
 669             }
 670             pimg->style = (leg->l.flags >> FLAGS_STYLE_BIT0) & 0x07;
 671             img_write_item(pimg, img_LINE, leg->l.flags & FLAGS_MASK,
 672                            sprint_prefix(leg_pfx),
 673                            POS(stn3, 0), POS(stn3, 1), POS(stn3, 2));
 674          }
 675
 676          /* FIXME: equate at the start of a traverse treated specially
 677           * - what about equates at end? */
 678          if (stn1->name != stn3->name && !(fEquate && cLegsTrav == 0)) {
 679             /* (node not part of same stn) &&
 680              * (not equate at start of traverse) */
 681 #ifndef BLUNDER_DETECTION
 682             if (fhErrStat && !fArtic) {
 683                if (!stn1->name->ident) {
 684                   /* FIXME: not ideal */
 685                   fputs("<fixed point>", fhErrStat);
 686                } else {
 687                   fprint_prefix(fhErrStat, stn1->name);
 688                }
 689                fputs(fEquate ? szLinkEq : szLink, fhErrStat);
 690                if (reached_end) {
 691                   if (!stn3->name->ident) {
 692                      /* FIXME: not ideal */
 693                      fputs("<fixed point>", fhErrStat);
 694                   } else {
 695                      fprint_prefix(fhErrStat, stn3->name);
 696                   }
 697                }
 698             }
 699 #endif
 700             if (!fEquate) {
 701                cLegsTrav++;
 702                lenTrav += sqrt(lenTot);
 703             }
 704          } else {
 705 #if SHOW_INTERNAL_LEGS
 706             if (fhErrStat && !fArtic) fprintf(fhErrStat, "+");
 707 #endif
 708             if (lenTot > 0.0) {
 709 #if DEBUG_INVALID
 710                fprintf(stderr, "lenTot = %8.4f ", lenTot);
 711                fprint_prefix(stderr, stn1->name);
 712                fprintf(stderr, " -> ");
 713                fprint_prefix(stderr, stn3->name);
 714 #endif
 715                BUG("during calculation of closure errors");
 716             }
 717          }
 718          if (reached_end) break;
 719
 720          i = k ^ 1; /* flip direction for other leg of 2 node */
 721
 722          stn1 = stn3;
 723       } /* endwhile */
 724
 725       if (cLegsTrav && !fArtic && fhErrStat)
 726          err_stat(cLegsTrav, lenTrav, eTot, eTotTheo,
 727                   hTot, hTotTheo, vTot, vTotTheo);
 728
 729       ptrOld = ptr;
 730       ptr = ptr->next;
 731       osfree(ptrOld);
 732    }
 733
 734    /* Leave fhErrStat open in case we're asked to close loops again... */
 735 }
 736
 737 static void
 738 err_stat(int cLegsTrav, double lenTrav,
 739          double eTot, double eTotTheo,
 740          double hTot, double hTotTheo,
 741          double vTot, double vTotTheo)
 742 {
 743    double E = sqrt(eTot / eTotTheo);
 744    double H = sqrt(hTot / hTotTheo);
 745    double V = sqrt(vTot / vTotTheo);
 746    if (!fSuppress) {
 747       double sqrt_eTot = sqrt(eTot);
 748       fputnl(fhErrStat);
 749       fprintf(fhErrStat, msg(/*Original length %6.2fm (%3d legs), moved %6.2fm (%5.2fm/leg). */145),
 750               lenTrav, cLegsTrav, sqrt_eTot, sqrt_eTot / cLegsTrav);
 751       if (lenTrav > 0.0) {
 752          fprintf(fhErrStat, msg(/*Error %6.2f%%*/146), 100 * sqrt_eTot / lenTrav);
 753       } else {
 754          /* TRANSLATORS: Here N/A means "Not Applicable" -- it means the
 755           * traverse has zero length, so error per metre is meaningless.
 756           *
 757           * There should be 4 spaces between "Error" and "N/A" so that it lines
 758           * up with the numbers in the message above. */
 759          fputs(msg(/*Error    N/A*/147), fhErrStat);
 760       }
 761       fputnl(fhErrStat);
 762       fprintf(fhErrStat, "%f\n", E);
 763       fprintf(fhErrStat, "H: %f V: %f\n", H, V);
 764       fputnl(fhErrStat);
 765    }
 766    img_write_errors(pimg, cLegsTrav, lenTrav, E, H, V);
 767 }
 768
 769 static void
 770 replace_trailing_travs(void)
 771 {
 772    stackTrail *ptrOld;
 773    node *stn1, *stn2;
 774    linkfor *leg;
 775    int i;
 776
 777    /* TRANSLATORS: In French, Eric chose to use the terminology used by
 778     * toporobot: "sequence" for the English "traverse", which makes sense
 779     * (although toporobot actually uses this term to mean something more
 780     * specific).  Feel free to follow this lead if you can't think of a better
 781     * term - these messages mostly indicate how processing is progressing.
 782     *
 783     * A trailing traverse is a dead end back to a junction. */
 784    out_current_action(msg(/*Calculating trailing traverses*/128));
 785
 786    while (ptrTrail != NULL) {
 787       leg = ptrTrail->join1;
 788       leg = reverse_leg(leg);
 789       stn1 = leg->l.to;
 790       i = reverse_leg_dirn(leg);
 791 #if PRINT_NETBITS
 792       printf(" Trailing trav ");
 793       print_prefix(stn1->name);
 794       printf("<%p>", stn1);
 795       printf("%s...\n", szLink);
 796       printf("    attachment stn is at (%f, %f, %f)\n",
 797              POS(stn1, 0), POS(stn1, 1), POS(stn1, 2));
 798 #endif
 799       /* We may have swapped the links round when we removed the leg.  If
 800        * we did then stn1->leg[i] will be in use.  The link we swapped
 801        * with is the first free leg */
 802       if (stn1->leg[i]) {
 803          /* j is the direction to swap with */
 804          int j = (stn1->leg[1]) ? 2 : 1;
 805          /* change the other direction of leg i to use leg j */
 806          reverse_leg(stn1->leg[i])->l.reverse += j - i;
 807          stn1->leg[j] = stn1->leg[i];
 808       }
 809       stn1->leg[i] = ptrTrail->join1;
 810       SVX_ASSERT(fixed(stn1));
 811       img_write_item(pimg, img_MOVE, 0, NULL,
 812                      POS(stn1, 0), POS(stn1, 1), POS(stn1, 2));
 813
 814       while (1) {
 815          prefix *leg_pfx;
 816          int j;
 817
 818          leg = stn1->leg[i];
 819          stn2 = leg->l.to;
 820          j = reverse_leg_dirn(leg);
 821          if (data_here(leg)) {
 822             leg_pfx = stn1->name->up;
 823             adddd(&POSD(stn2), &POSD(stn1), &leg->d);
 824 #if 0
 825             printf("Adding leg (%f, %f, %f)\n", leg->d[0], leg->d[1], leg->d[2]);
 826 #endif
 827          } else {
 828             leg_pfx = stn2->name->up;
 829             leg = stn2->leg[j];
 830             subdd(&POSD(stn2), &POSD(stn1), &leg->d);
 831 #if 0
 832             printf("Subtracting reverse leg (%f, %f, %f)\n", leg->d[0], leg->d[1], leg->d[2]);
 833 #endif
 834          }
 835
 836          fix(stn2);
 837          add_stn_to_list(&stnlist, stn2);
 838          if (!(leg->l.reverse & (FLAG_REPLACEMENTLEG | FLAG_FAKE))) {
 839              if (TSTBIT(leg->l.flags, FLAGS_SURFACE)) {
 840                 stn1->name->sflags |= BIT(SFLAGS_SURFACE);
 841                 stn2->name->sflags |= BIT(SFLAGS_SURFACE);
 842              } else {
 843                 stn1->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 844                 stn2->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 845              }
 846          }
 847          if (!(leg->l.reverse & (FLAG_REPLACEMENTLEG | FLAG_FAKE))) {
 848             SVX_ASSERT(!fZeros(&leg->v));
 849             if (leg->meta) {
 850                 pimg->days1 = leg->meta->days1;
 851                 pimg->days2 = leg->meta->days2;
 852             } else {
 853                 pimg->days1 = pimg->days2 = -1;
 854             }
 855             pimg->style = (leg->l.flags >> FLAGS_STYLE_BIT0) & 0x07;
 856             img_write_item(pimg, img_LINE, leg->l.flags & FLAGS_MASK,
 857                            sprint_prefix(leg_pfx),
 858                            POS(stn2, 0), POS(stn2, 1), POS(stn2, 2));
 859          }
 860
 861          /* stop if not 2 node */
 862          if (!two_node(stn2)) break;
 863
 864          stn1 = stn2;
 865          i = j ^ 1; /* flip direction for other leg of 2 node */
 866       }
 867
 868       ptrOld = ptrTrail;
 869       ptrTrail = ptrTrail->next;
 870       osfree(ptrOld);
 871    }
 872
 873    /* write out connections with no survey data */
 874    while (nosurveyhead) {
 875       nosurveylink *p = nosurveyhead;
 876       SVX_ASSERT(fixed(p->fr));
 877       SVX_ASSERT(fixed(p->to));
 878       if (TSTBIT(p->flags, FLAGS_SURFACE)) {
 879          p->fr->name->sflags |= BIT(SFLAGS_SURFACE);
 880          p->to->name->sflags |= BIT(SFLAGS_SURFACE);
 881       } else {
 882          p->fr->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 883          p->to->name->sflags |= BIT(SFLAGS_UNDERGROUND);
 884       }
 885       img_write_item(pimg, img_MOVE, 0, NULL,
 886                      POS(p->fr, 0), POS(p->fr, 1), POS(p->fr, 2));
 887       if (p->meta) {
 888           pimg->days1 = p->meta->days1;
 889           pimg->days2 = p->meta->days2;
 890       } else {
 891           pimg->days1 = pimg->days2 = -1;
 892       }
 893       pimg->style = img_STYLE_NOSURVEY;
 894       img_write_item(pimg, img_LINE, (p->flags & FLAGS_MASK),
 895                      sprint_prefix(p->fr->name->up),
 896                      POS(p->to, 0), POS(p->to, 1), POS(p->to, 2));
 897       nosurveyhead = p->next;
 898       osfree(p);
 899    }
 900
 901    /* write stations to .3d file and free legs and stations */
 902    FOR_EACH_STN(stn1, stnlist) {
 903       int d;
 904       SVX_ASSERT(fixed(stn1));
 905       if (stn1->name->stn == stn1) {
 906          int sf = stn1->name->sflags;
 907          /* take care of unused fixed points */
 908          if (!TSTBIT(sf, SFLAGS_SOLVED)) {
 909             const char * label = NULL;
 910             if (TSTBIT(sf, SFLAGS_ANON)) {
 911                label = "";
 912             } else if (stn1->name->ident) {
 913                label = sprint_prefix(stn1->name);
 914             }
 915             if (label) {
 916                /* Set flag to stop station being rewritten after *solve. */
 917                stn1->name->sflags = sf | BIT(SFLAGS_SOLVED);
 918                sf &= SFLAGS_MASK;
 919                if (stn1->name->max_export) sf |= BIT(SFLAGS_EXPORTED);
 920                img_write_item(pimg, img_LABEL, sf, label,
 921                               POS(stn1, 0), POS(stn1, 1), POS(stn1, 2));
 922             }
 923          }
 924       }
 925       /* update coords of bounding box, ignoring the base positions
 926        * of points fixed with error estimates and only counting stations
 927        * in underground surveys.
 928        *
 929        * NB We don't set SFLAGS_UNDERGROUND for the anchor station for
 930        * a point fixed with error estimates, so this test will exclude
 931        * those too, which is what we want.
 932        */
 933       if (TSTBIT(stn1->name->sflags, SFLAGS_UNDERGROUND)) {
 934          for (d = 0; d < 3; d++) {
 935             if (POS(stn1, d) < min[d]) {
 936                min[d] = POS(stn1, d);
 937                pfxLo[d] = stn1->name;
 938             }
 939             if (POS(stn1, d) > max[d]) {
 940                max[d] = POS(stn1, d);
 941                pfxHi[d] = stn1->name;
 942             }
 943          }
 944
 945          /* Range without anonymous stations at offset 3. */
 946          if (!TSTBIT(stn1->name->sflags, SFLAGS_ANON)) {
 947             for (d = 0; d < 3; d++) {
 948                if (POS(stn1, d) < min[d + 3]) {
 949                   min[d + 3] = POS(stn1, d);
 950                   pfxLo[d + 3] = stn1->name;
 951                }
 952                if (POS(stn1, d) > max[d + 3]) {
 953                   max[d + 3] = POS(stn1, d);
 954                   pfxHi[d + 3] = stn1->name;
 955                }
 956             }
 957          }
 958       }
 959
 960       d = stn1->name->shape;
 961       if (d <= 1 && !TSTBIT(stn1->name->sflags, SFLAGS_USED)) {
 962          bool unused_fixed_point = false;
 963          if (d == 0) {
 964             /* Unused fixed point without error estimates */
 965             unused_fixed_point = true;
 966          } else if (stn1->leg[0]) {
 967             prefix *pfx = stn1->leg[0]->l.to->name;
 968             if (!pfx->ident && !TSTBIT(pfx->sflags, SFLAGS_ANON)) {
 969                /* Unused fixed point with error estimates */
 970                unused_fixed_point = true;
 971             }
 972          }
 973          if (unused_fixed_point) {
 974             /* TRANSLATORS: fixed survey station that is not part of any survey
 975              */
 976             warning_in_file(stn1->name->filename, stn1->name->line,
 977                     /*Unused fixed point “%s”*/73, sprint_prefix(stn1->name));
 978          }
 979       }
 980
 981       /* For stations fixed with error estimates, we need to ignore the leg to
 982        * the "real" fixed point in the node stats.
 983        */
 984       if (stn1->leg[0] && !stn1->leg[0]->l.to->name->ident &&
 985           !TSTBIT(stn1->leg[0]->l.to->name->sflags, SFLAGS_ANON))
 986          stn1->name->shape--;
 987
 988       for (i = 0; i <= 2; i++) {
 989          leg = stn1->leg[i];
 990          /* only want to think about forwards legs */
 991          if (leg && data_here(leg)) {
 992             linkfor *legRev;
 993             node *stnB;
 994             int iB;
 995             stnB = leg->l.to;
 996             iB = reverse_leg_dirn(leg);
 997             legRev = stnB->leg[iB];
 998             SVX_ASSERT2(legRev->l.to == stn1, "leg doesn't reciprocate");
 999             SVX_ASSERT(fixed(stn1));
1000             if (!(leg->l.flags &
1001                   (BIT(FLAGS_DUPLICATE)|BIT(FLAGS_SPLAY)|
1002                    BIT(FLAGS_SURFACE)))) {
1003                /* check not an equating leg, or one inside an sdfix point */
1004                if (!(leg->l.reverse & (FLAG_REPLACEMENTLEG | FLAG_FAKE))) {
1005                   totadj += sqrt(sqrd(POS(stnB, 0) - POS(stn1, 0)) +
1006                                  sqrd(POS(stnB, 1) - POS(stn1, 1)) +
1007                                  sqrd(POS(stnB, 2) - POS(stn1, 2)));
1008                   total += sqrt(sqrdd(leg->d));
1009                   totplan += hypot(leg->d[0], leg->d[1]);
1010                   totvert += fabs(leg->d[2]);
1011                }
1012             }
1013             osfree(leg);
1014             osfree(legRev);
1015             stn1->leg[i] = stnB->leg[iB] = NULL;
1016          }
1017       }
1018    }
1019
1020    /* The station position is attached to the name, so we leave the names and
1021     * positions in place - they can then be picked up if we have a *solve
1022     * followed by more data */
1023    for (stn1 = stnlist; stn1; stn1 = stn2) {
1024       stn2 = stn1->next;
1025       stn1->name->stn = NULL;
1026       osfree(stn1);
1027    }
1028    stnlist = NULL;
1029 }
1030
1031 static void
1032 write_passage_models(void)
1033 {
1034    lrudlist * psg = model;
1035    while (psg) {
1036       lrudlist * oldp;
1037       lrud * xsect = psg->tube;
1038       int xflags = 0;
1039       while (xsect) {
1040          lrud *oldx;
1041          prefix *pfx;
1042          const char *name;
1043          pimg->l = xsect->l;
1044          pimg->r = xsect->r;
1045          pimg->u = xsect->u;
1046          pimg->d = xsect->d;
1047          if (xsect->meta) {
1048              pimg->days1 = xsect->meta->days1;
1049              pimg->days2 = xsect->meta->days2;
1050          } else {
1051              pimg->days1 = pimg->days2 = -1;
1052          }
1053
1054          pfx = xsect->stn;
1055          name = sprint_prefix(pfx);
1056          oldx = xsect;
1057          xsect = xsect->next;
1058          osfree(oldx);
1059
1060          if (!pfx->pos) {
1061              /* TRANSLATORS: e.g. the user specifies a passage cross-section at
1062               * station "entrance.27", but there is no station "entrance.27" in
1063               * the centre-line. */
1064              error_in_file(pfx->filename, pfx->line,
1065                            /*Cross section specified at non-existent station “%s”*/83,
1066                            name);
1067          } else {
1068              if (xsect == NULL) xflags = img_XFLAG_END;
1069              img_write_item(pimg, img_XSECT, xflags, name, 0, 0, 0);
1070          }
1071       }
1072       oldp = psg;
1073       psg = psg->next;
1074       osfree(oldp);
1075    }
1076    model = NULL;
1077 }