source: git/src/netskel.c @ ad83822

RELEASE/1.0RELEASE/1.2debug-cidebug-ci-sanitisersfaster-cavernlogstereowalls-datawalls-data-hanging-as-warning
Last change on this file since ad83822 was fbae3de, checked in by Olly Betts <olly@…>, 19 years ago

Fix assertion failure when handling a survey with a station only attached
by a nosurvey leg.

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