source: git/src/netskel.c @ 168df28

RELEASE/1.0RELEASE/1.2debug-cidebug-ci-sanitisersstereowalls-data
Last change on this file since 168df28 was 168df28, checked in by Olly Betts <olly@…>, 18 years ago

Backport a whole pile of fixes and minor tweaks from 1.1.7.

git-svn-id: file:///home/survex-svn/survex/trunk@3114 4b37db11-9a0c-4f06-9ece-9ab7cdaee568

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