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