source: git/src/netskel.c @ a72ed95

RELEASE/1.2debug-cidebug-ci-sanitisersstereowalls-data
Last change on this file since a72ed95 was 84d5946, checked in by Olly Betts <olly@…>, 10 years ago

src/netskel.c: If there's an output coordinate system set, write it
to the 3d file.

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