source: git/src/datain.c @ dda0ca7

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

src/datain.c,src/netbits.c,src/netbits.h,tests/: If the same leg is
repeated consecutively, average the readings and treat as a single
leg.

  • Property mode set to 100644
File size: 49.9 KB
Line 
1/* datain.c
2 * Reads in survey files, dealing with special characters, keywords & data
3 * Copyright (C) 1991-2003,2005,2009,2010,2011,2012,2013,2014,2015 Olly Betts
4 * Copyright (C) 2004 Simeon Warner
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#ifdef HAVE_CONFIG_H
22#include <config.h>
23#endif
24
25#include <limits.h>
26#include <stdarg.h>
27
28#include "debug.h"
29#include "cavern.h"
30#include "date.h"
31#include "filename.h"
32#include "message.h"
33#include "filelist.h"
34#include "netbits.h"
35#include "netskel.h"
36#include "readval.h"
37#include "datain.h"
38#include "commands.h"
39#include "out.h"
40#include "str.h"
41
42#define EPSILON (REAL_EPSILON * 1000)
43
44#define var(I) (pcs->Var[(I)])
45
46/* true if x is not-a-number value in Compass (999.0 or -999.0)    */
47/* Compass uses 999.0 but understands Karst data which used -999.0 */
48#define is_compass_NaN(x) ( fabs(fabs(x)-999.0) <  EPSILON )
49
50int ch;
51
52typedef enum {
53    CTYPE_OMIT, CTYPE_READING, CTYPE_PLUMB, CTYPE_INFERPLUMB, CTYPE_HORIZ
54} clino_type;
55
56/* Don't explicitly initialise as we can't set the jmp_buf - this has
57 * static scope so will be initialised like this anyway */
58parse file /* = { NULL, NULL, 0, fFalse, NULL } */ ;
59
60bool f_export_ok;
61
62static real value[Fr - 1];
63#define VAL(N) value[(N)-1]
64static real variance[Fr - 1];
65#define VAR(N) variance[(N)-1]
66
67/* style functions */
68static void data_normal(void);
69static void data_cartesian(void);
70static void data_passage(void);
71static void data_nosurvey(void);
72static void data_ignore(void);
73
74void
75get_pos(filepos *fp)
76{
77   fp->ch = ch;
78   fp->offset = ftell(file.fh);
79   if (fp->offset == -1)
80      fatalerror_in_file(file.filename, 0, /*Error reading file*/18);
81}
82
83void
84set_pos(const filepos *fp)
85{
86   ch = fp->ch;
87   if (fseek(file.fh, fp->offset, SEEK_SET) == -1)
88      fatalerror_in_file(file.filename, 0, /*Error reading file*/18);
89}
90
91static void
92push_back(int c)
93{
94   if (c != EOF && ungetc(c, file.fh) == EOF)
95      fatalerror_in_file(file.filename, 0, /*Error reading file*/18);
96}
97
98static void
99report_parent(parse * p) {
100    if (p->parent)
101        report_parent(p->parent);
102    /* Force re-report of include tree for further errors in
103     * parent files */
104    p->reported_where = fFalse;
105    /* TRANSLATORS: %s is replaced by the filename of the parent file, and %u
106     * by the line number in that file. */
107    fprintf(STDERR, msg(/*In file included from %s:%u:\n*/5), p->filename, p->line);
108}
109
110static void
111error_list_parent_files(void)
112{
113   if (!file.reported_where && file.parent) {
114      report_parent(file.parent);
115      /* Suppress reporting of full include tree for further errors
116       * in this file */
117      file.reported_where = fTrue;
118   }
119}
120
121void
122compile_error(int en, ...)
123{
124   int col = 0;
125   va_list ap;
126   va_start(ap, en);
127   error_list_parent_files();
128   if (en < 0) {
129      en = -en;
130      col = ftell(file.fh) - file.lpos;
131   }
132   v_report(1, file.filename, file.line, col, en, ap);
133   va_end(ap);
134}
135
136void
137compile_error_skip(int en, ...)
138{
139   int col = 0;
140   va_list ap;
141   va_start(ap, en);
142   error_list_parent_files();
143   if (en < 0) {
144      en = -en;
145      col = ftell(file.fh) - file.lpos;
146   }
147   v_report(1, file.filename, file.line, col, en, ap);
148   va_end(ap);
149   skipline();
150}
151
152void
153compile_error_at(const char * filename, unsigned line, int en, ...)
154{
155   va_list ap;
156   va_start(ap, en);
157   v_report(1, filename, line, 0, en, ap);
158   va_end(ap);
159}
160
161void
162compile_error_pfx(const prefix * pfx, int en, ...)
163{
164   va_list ap;
165   va_start(ap, en);
166   v_report(1, pfx->filename, pfx->line, 0, en, ap);
167   va_end(ap);
168}
169
170void
171compile_error_token(int en)
172{
173   char *p = NULL;
174   int len = 0;
175   skipblanks();
176   while (!isBlank(ch) && !isEol(ch)) {
177      s_catchar(&p, &len, (char)ch);
178      nextch();
179   }
180   compile_error(en, p ? p : "");
181   osfree(p);
182}
183
184void
185compile_warning(int en, ...)
186{
187   int col = 0;
188   va_list ap;
189   va_start(ap, en);
190   error_list_parent_files();
191   if (en < 0) {
192      en = -en;
193      col = ftell(file.fh) - file.lpos;
194   }
195   v_report(0, file.filename, file.line, col, en, ap);
196   va_end(ap);
197}
198
199void
200compile_warning_at(const char * filename, unsigned line, int en, ...)
201{
202   va_list ap;
203   va_start(ap, en);
204   v_report(0, filename, line, 0, en, ap);
205   va_end(ap);
206}
207
208void
209compile_warning_pfx(const prefix * pfx, int en, ...)
210{
211   va_list ap;
212   va_start(ap, en);
213   v_report(0, pfx->filename, pfx->line, 0, en, ap);
214   va_end(ap);
215}
216
217/* This function makes a note where to put output files */
218static void
219using_data_file(const char *fnm)
220{
221   if (!fnm_output_base) {
222      /* was: fnm_output_base = base_from_fnm(fnm); */
223      fnm_output_base = baseleaf_from_fnm(fnm);
224   } else if (fnm_output_base_is_dir) {
225      /* --output pointed to directory so use the leaf basename in that dir */
226      char *lf, *p;
227      lf = baseleaf_from_fnm(fnm);
228      p = use_path(fnm_output_base, lf);
229      osfree(lf);
230      osfree(fnm_output_base);
231      fnm_output_base = p;
232      fnm_output_base_is_dir = 0;
233   }
234}
235
236static void
237skipword(void)
238{
239   while (!isBlank(ch) && !isEol(ch)) nextch();
240}
241
242extern void
243skipblanks(void)
244{
245   while (isBlank(ch)) nextch();
246}
247
248extern void
249skipline(void)
250{
251   while (!isEol(ch)) nextch();
252}
253
254#ifndef NO_PERCENTAGE
255static long int filelen;
256#endif
257
258static void
259process_bol(void)
260{
261#ifndef NO_PERCENTAGE
262   /* print %age of file done */
263   if (filelen > 0) {
264      filepos fp;
265      get_pos(&fp);
266      printf("%d%%\r", (int)(100 * fp.offset / filelen));
267   }
268#endif
269
270   nextch();
271   skipblanks();
272}
273
274static void
275process_eol(void)
276{
277   int eolchar;
278
279   skipblanks();
280
281   if (!isEol(ch)) {
282      if (!isComm(ch)) compile_error(-/*End of line not blank*/15);
283      skipline();
284   }
285
286   eolchar = ch;
287   file.line++;
288   /* skip any different eol characters so we get line counts correct on
289    * DOS text files and similar, but don't count several adjacent blank
290    * lines as one */
291   while (ch != EOF) {
292      nextch();
293      if (ch == eolchar || !isEol(ch)) {
294         push_back(ch);
295         break;
296      }
297   }
298   file.lpos = ftell(file.fh);
299}
300
301static bool
302process_non_data_line(void)
303{
304   process_bol();
305
306   if (isData(ch)) return fFalse;
307
308   if (isKeywd(ch)) {
309      nextch();
310      handle_command();
311   }
312
313   process_eol();
314
315   return fTrue;
316}
317
318static void
319read_reading(reading r, bool f_optional)
320{
321   int n_readings;
322   q_quantity q;
323   switch (r) {
324      case Tape: q = Q_LENGTH; break;
325      case Comp: q = Q_BEARING; break;
326      case BackComp: q = Q_BACKBEARING; break;
327      case Clino: q = Q_GRADIENT; break;
328      case BackClino: q = Q_BACKGRADIENT; break;
329      case FrDepth: case ToDepth: q = Q_DEPTH; break;
330      case Dx: q = Q_DX; break;
331      case Dy: q = Q_DY; break;
332      case Dz: q = Q_DZ; break;
333      case FrCount: case ToCount: q = Q_COUNT; break;
334      case Left: q = Q_LEFT; break;
335      case Right: q = Q_RIGHT; break;
336      case Up: q = Q_UP; break;
337      case Down: q = Q_DOWN; break;
338      default:
339        q = Q_NULL; /* Suppress compiler warning */;
340        BUG("Unexpected case");
341   }
342   VAL(r) = read_numeric(f_optional, &n_readings);
343   VAR(r) = var(q);
344   if (n_readings > 1) VAR(r) /= sqrt(n_readings);
345}
346
347static void
348read_bearing_or_omit(reading r)
349{
350   int n_readings;
351   q_quantity q = Q_NULL;
352   VAL(r) = read_numeric_or_omit(&n_readings);
353   switch (r) {
354      case Comp: q = Q_BEARING; break;
355      case BackComp: q = Q_BACKBEARING; break;
356      default:
357        q = Q_NULL; /* Suppress compiler warning */;
358        BUG("Unexpected case");
359   }
360   VAR(r) = var(q);
361   if (n_readings > 1) VAR(r) /= sqrt(n_readings);
362}
363
364/* For reading Compass MAK files which have a freeform syntax */
365static void
366nextch_handling_eol(void)
367{
368   nextch();
369   while (ch != EOF && isEol(ch)) {
370      process_eol();
371      nextch();
372   }
373}
374
375#define LITLEN(S) (sizeof(S"") - 1)
376#define has_ext(F,L,E) ((L) > LITLEN(E) + 1 &&\
377                        (F)[(L) - LITLEN(E) - 1] == FNM_SEP_EXT &&\
378                        strcasecmp((F) + (L) - LITLEN(E), E) == 0)
379extern void
380data_file(const char *pth, const char *fnm)
381{
382   int begin_lineno_store;
383   parse file_store;
384   volatile enum {FMT_SVX, FMT_DAT, FMT_MAK} fmt = FMT_SVX;
385
386   {
387      char *filename;
388      FILE *fh;
389      size_t len;
390
391      if (!pth) {
392         /* file specified on command line - don't do special translation */
393         fh = fopenWithPthAndExt(pth, fnm, EXT_SVX_DATA, "rb", &filename);
394      } else {
395         fh = fopen_portable(pth, fnm, EXT_SVX_DATA, "rb", &filename);
396      }
397
398      if (fh == NULL) {
399         compile_error(/*Couldn’t open file “%s”*/24, fnm);
400         return;
401      }
402
403      len = strlen(filename);
404      if (has_ext(filename, len, "dat")) {
405         fmt = FMT_DAT;
406      } else if (has_ext(filename, len, "mak")) {
407         fmt = FMT_MAK;
408      }
409
410      file_store = file;
411      if (file.fh) file.parent = &file_store;
412      file.fh = fh;
413      file.filename = filename;
414      file.line = 1;
415      file.lpos = 0;
416      file.reported_where = fFalse;
417   }
418
419   if (fPercent) printf("%s:\n", fnm);
420
421   using_data_file(file.filename);
422
423   begin_lineno_store = pcs->begin_lineno;
424   pcs->begin_lineno = 0;
425
426#ifndef NO_PERCENTAGE
427   /* Try to find how long the file is...
428    * However, under ANSI fseek( ..., SEEK_END) may not be supported */
429   filelen = 0;
430   if (fPercent) {
431      if (fseek(file.fh, 0l, SEEK_END) == 0) {
432         filepos fp;
433         get_pos(&fp);
434         filelen = fp.offset;
435      }
436      rewind(file.fh); /* reset file ptr to start & clear any error state */
437   }
438#endif
439
440   if (fmt == FMT_DAT) {
441      short *t;
442      int i;
443      settings *pcsNew;
444
445      pcsNew = osnew(settings);
446      *pcsNew = *pcs; /* copy contents */
447      pcsNew->begin_lineno = 0;
448      pcsNew->next = pcs;
449      pcs = pcsNew;
450      default_units(pcs);
451      default_calib(pcs);
452
453      pcs->style = STYLE_NORMAL;
454      pcs->units[Q_LENGTH] = METRES_PER_FOOT;
455      t = ((short*)osmalloc(ossizeof(short) * 257)) + 1;
456
457      t[EOF] = SPECIAL_EOL;
458      memset(t, 0, sizeof(short) * 33);
459      for (i = 33; i < 127; i++) t[i] = SPECIAL_NAMES;
460      t[127] = 0;
461      for (i = 128; i < 256; i++) t[i] = SPECIAL_NAMES;
462      t['\t'] |= SPECIAL_BLANK;
463      t[' '] |= SPECIAL_BLANK;
464      t['\032'] |= SPECIAL_EOL; /* Ctrl-Z, so olde DOS text files are handled ok */
465      t['\n'] |= SPECIAL_EOL;
466      t['\r'] |= SPECIAL_EOL;
467      t['.'] |= SPECIAL_DECIMAL;
468      t['-'] |= SPECIAL_MINUS;
469      t['+'] |= SPECIAL_PLUS;
470      pcs->Translate = t;
471      pcs->Case = OFF;
472      pcs->Truncate = INT_MAX;
473      pcs->infer = BIT(INFER_EQUATES)|BIT(INFER_EXPORTS)|BIT(INFER_PLUMBS);
474   } else if (fmt == FMT_MAK) {
475      short *t;
476      int i;
477      settings *pcsNew;
478
479      pcsNew = osnew(settings);
480      *pcsNew = *pcs; /* copy contents */
481      pcsNew->begin_lineno = 0;
482      pcsNew->next = pcs;
483      pcs = pcsNew;
484
485      t = ((short*)osmalloc(ossizeof(short) * 257)) + 1;
486
487      t[EOF] = SPECIAL_EOL;
488      memset(t, 0, sizeof(short) * 33);
489      for (i = 33; i < 127; i++) t[i] = SPECIAL_NAMES;
490      t[127] = 0;
491      for (i = 128; i < 256; i++) t[i] = SPECIAL_NAMES;
492      t['['] = t[','] = t[';'] = 0;
493      t['\t'] |= SPECIAL_BLANK;
494      t[' '] |= SPECIAL_BLANK;
495      t['\032'] |= SPECIAL_EOL; /* Ctrl-Z, so olde DOS text files are handled ok */
496      t['\n'] |= SPECIAL_EOL;
497      t['\r'] |= SPECIAL_EOL;
498      t['.'] |= SPECIAL_DECIMAL;
499      t['-'] |= SPECIAL_MINUS;
500      t['+'] |= SPECIAL_PLUS;
501      pcs->Translate = t;
502      pcs->Case = OFF;
503      pcs->Truncate = INT_MAX;
504   }
505
506#ifdef HAVE_SETJMP_H
507   /* errors in nested functions can longjmp here */
508   if (setjmp(file.jbSkipLine)) {
509      skipline();
510      process_eol();
511   }
512#endif
513
514   if (fmt == FMT_DAT) {
515      while (!feof(file.fh) && !ferror(file.fh)) {
516         static reading compass_order[] = {
517            Fr, To, Tape, CompassDATComp, CompassDATClino,
518            CompassDATLeft, CompassDATRight, CompassDATUp, CompassDATDown,
519            CompassDATFlags, IgnoreAll
520         };
521         static reading compass_order_backsights[] = {
522            Fr, To, Tape, CompassDATComp, CompassDATClino,
523            CompassDATLeft, CompassDATRight, CompassDATUp, CompassDATDown,
524            CompassDATBackComp, CompassDATBackClino,
525            CompassDATFlags, IgnoreAll
526         };
527         /* <Cave name> */
528         process_bol();
529         skipline();
530         process_eol();
531         /* SURVEY NAME: <Short name> */
532         get_token();
533         get_token();
534         /* if (ch != ':') ... */
535         nextch();
536         get_token();
537         skipline();
538         process_eol();
539         /* SURVEY DATE: 7 10 79  COMMENT:<Long name> */
540         get_token();
541         get_token();
542         copy_on_write_meta(pcs);
543         if (ch == ':') {
544             int year, month, day;
545
546             nextch();
547
548             /* NB order is *month* *day* year */
549             month = read_uint();
550             day = read_uint();
551             year = read_uint();
552             /* Note: Larry says a 2 digit year is always 19XX */
553             if (year < 100) year += 1900;
554
555             pcs->meta->days1 = pcs->meta->days2 = days_since_1900(year, month, day);
556         } else {
557             pcs->meta->days1 = pcs->meta->days2 = -1;
558         }
559         skipline();
560         process_eol();
561         /* SURVEY TEAM: */
562         get_token();
563         get_token();
564         skipline();
565         process_eol();
566         /* <Survey team> */
567         nextch();
568         skipline();
569         process_eol();
570         /* DECLINATION: 1.00  FORMAT: DDDDLUDRADLN  CORRECTIONS: 2.00 3.00 4.00 */
571         get_token();
572         nextch(); /* : */
573         skipblanks();
574         pcs->z[Q_DECLINATION] = -read_numeric(fFalse, NULL);
575         pcs->z[Q_DECLINATION] *= pcs->units[Q_DECLINATION];
576         get_token();
577         pcs->ordering = compass_order;
578         if (strcmp(buffer, "FORMAT") == 0) {
579            nextch(); /* : */
580            get_token();
581            if (strlen(buffer) >= 12 && buffer[11] == 'B') {
582               /* We have backsights for compass and clino */
583               pcs->ordering = compass_order_backsights;
584            }
585            get_token();
586         }
587         if (strcmp(buffer, "CORRECTIONS") == 0) {
588            nextch(); /* : */
589            pcs->z[Q_BEARING] = -rad(read_numeric(fFalse, NULL));
590            pcs->z[Q_GRADIENT] = -rad(read_numeric(fFalse, NULL));
591            pcs->z[Q_LENGTH] = -read_numeric(fFalse, NULL);
592         } else {
593            pcs->z[Q_BEARING] = 0;
594            pcs->z[Q_GRADIENT] = 0;
595            pcs->z[Q_LENGTH] = 0;
596         }
597         skipline();
598         process_eol();
599         /* BLANK LINE */
600         process_bol();
601         skipline();
602         process_eol();
603         /* heading line */
604         process_bol();
605         skipline();
606         process_eol();
607         /* BLANK LINE */
608         process_bol();
609         skipline();
610         process_eol();
611         while (!feof(file.fh)) {
612            process_bol();
613            if (ch == '\x0c') {
614               nextch();
615               process_eol();
616               break;
617            }
618            data_normal();
619         }
620         clear_last_leg();
621      }
622      {
623         settings *pcsParent = pcs->next;
624         SVX_ASSERT(pcsParent);
625         pcs->ordering = NULL;
626         free_settings(pcs);
627         pcs = pcsParent;
628      }
629   } else if (fmt == FMT_MAK) {
630      nextch_handling_eol();
631      while (!feof(file.fh) && !ferror(file.fh)) {
632         if (ch == '#') {
633            /* include a file */
634            int ch_store;
635            char *dat_pth = path_from_fnm(file.filename);
636            char *dat_fnm = NULL;
637            int dat_fnm_len;
638            nextch_handling_eol();
639            while (ch != ',' && ch != ';' && ch != EOF) {
640               while (isEol(ch)) process_eol();
641               s_catchar(&dat_fnm, &dat_fnm_len, (char)ch);
642               nextch_handling_eol();
643            }
644            while (ch != ';' && ch != EOF) {
645               prefix *name;
646               nextch_handling_eol();
647               name = read_prefix(PFX_STATION|PFX_OPT);
648               if (name) {
649                  skipblanks();
650                  if (ch == '[') {
651                     /* fixed pt */
652                     node *stn;
653                     real x, y, z;
654                     name->sflags |= BIT(SFLAGS_FIXED);
655                     nextch_handling_eol();
656                     while (!isdigit(ch) && ch != '+' && ch != '-' &&
657                            ch != '.' && ch != ']' && ch != EOF) {
658                        nextch_handling_eol();
659                     }
660                     x = read_numeric(fFalse, NULL);
661                     while (!isdigit(ch) && ch != '+' && ch != '-' &&
662                            ch != '.' && ch != ']' && ch != EOF) {
663                        nextch_handling_eol();
664                     }
665                     y = read_numeric(fFalse, NULL);
666                     while (!isdigit(ch) && ch != '+' && ch != '-' &&
667                            ch != '.' && ch != ']' && ch != EOF) {
668                        nextch_handling_eol();
669                     }
670                     z = read_numeric(fFalse, NULL);
671                     stn = StnFromPfx(name);
672                     if (!fixed(stn)) {
673                        POS(stn, 0) = x;
674                        POS(stn, 1) = y;
675                        POS(stn, 2) = z;
676                        fix(stn);
677                     } else {
678                        if (x != POS(stn, 0) || y != POS(stn, 1) ||
679                            z != POS(stn, 2)) {
680                           compile_error(/*Station already fixed or equated to a fixed point*/46);
681                        } else {
682                           compile_warning(/*Station already fixed at the same coordinates*/55);
683                        }
684                     }
685                     while (ch != ']' && ch != EOF) nextch_handling_eol();
686                     if (ch == ']') {
687                        nextch_handling_eol();
688                        skipblanks();
689                     }
690                  } else {
691                     /* FIXME: link station - ignore for now */
692                     /* FIXME: perhaps issue warning? */
693                  }
694                  while (ch != ',' && ch != ';' && ch != EOF)
695                     nextch_handling_eol();
696               }
697            }
698            if (dat_fnm) {
699               ch_store = ch;
700               data_file(dat_pth, dat_fnm);
701               ch = ch_store;
702               osfree(dat_fnm);
703            }
704         } else {
705            /* FIXME: also check for % and $ later */
706            nextch_handling_eol();
707         }
708      }
709      {
710         settings *pcsParent = pcs->next;
711         SVX_ASSERT(pcsParent);
712         free_settings(pcs);
713         pcs = pcsParent;
714      }
715   } else {
716      while (!feof(file.fh) && !ferror(file.fh)) {
717         if (!process_non_data_line()) {
718            f_export_ok = fFalse;
719            switch (pcs->style) {
720             case STYLE_NORMAL:
721             case STYLE_DIVING:
722             case STYLE_CYLPOLAR:
723               data_normal();
724               break;
725             case STYLE_CARTESIAN:
726               data_cartesian();
727               break;
728             case STYLE_PASSAGE:
729               data_passage();
730               break;
731             case STYLE_NOSURVEY:
732               data_nosurvey();
733               break;
734             case STYLE_IGNORE:
735               data_ignore();
736               break;
737             default:
738               BUG("bad style");
739            }
740         }
741      }
742      clear_last_leg();
743   }
744
745   /* don't allow *BEGIN at the end of a file, then *EXPORT in the
746    * including file */
747   f_export_ok = fFalse;
748
749#ifndef NO_PERCENTAGE
750   if (fPercent) putnl();
751#endif
752
753   if (pcs->begin_lineno) {
754      error_in_file(file.filename, pcs->begin_lineno,
755                    /*BEGIN with no matching END in this file*/23);
756      /* Implicitly close any unclosed BEGINs from this file */
757      do {
758         settings *pcsParent = pcs->next;
759         SVX_ASSERT(pcsParent);
760         free_settings(pcs);
761         pcs = pcsParent;
762      } while (pcs->begin_lineno);
763   }
764
765   pcs->begin_lineno = begin_lineno_store;
766
767   if (ferror(file.fh))
768      fatalerror_in_file(file.filename, 0, /*Error reading file*/18);
769
770   (void)fclose(file.fh);
771
772   file = file_store;
773
774   /* don't free this - it may be pointed to by prefix.file */
775   /* osfree(file.filename); */
776}
777
778static real
779mod2pi(real a)
780{
781   return a - floor(a / (2 * M_PI)) * (2 * M_PI);
782}
783
784static real
785handle_plumb(clino_type *p_ctype)
786{
787   typedef enum {
788      CLINO_NULL=-1, CLINO_UP, CLINO_DOWN, CLINO_LEVEL
789   } clino_tok;
790   static sztok clino_tab[] = {
791      {"D",     CLINO_DOWN},
792      {"DOWN",  CLINO_DOWN},
793      {"H",     CLINO_LEVEL},
794      {"LEVEL", CLINO_LEVEL},
795      {"U",     CLINO_UP},
796      {"UP",    CLINO_UP},
797      {NULL,    CLINO_NULL}
798   };
799   static real clinos[] = {(real)M_PI_2, (real)(-M_PI_2), (real)0.0};
800   clino_tok tok;
801
802   skipblanks();
803   if (isalpha(ch)) {
804      filepos fp;
805      get_pos(&fp);
806      get_token();
807      tok = match_tok(clino_tab, TABSIZE(clino_tab));
808      if (tok != CLINO_NULL) {
809         *p_ctype = (tok == CLINO_LEVEL ? CTYPE_HORIZ : CTYPE_PLUMB);
810         return clinos[tok];
811      }
812      set_pos(&fp);
813   } else if (isSign(ch)) {
814      int chOld = ch;
815      nextch();
816      if (toupper(ch) == 'V') {
817         nextch();
818         *p_ctype = CTYPE_PLUMB;
819         return (!isMinus(chOld) ? M_PI_2 : -M_PI_2);
820      }
821
822      if (isOmit(chOld)) {
823         *p_ctype = CTYPE_OMIT;
824         /* no clino reading, so assume 0 with large sd */
825         return (real)0.0;
826      }
827   } else if (isOmit(ch)) {
828      /* OMIT char may not be a SIGN char too so we need to check here as
829       * well as above... */
830      nextch();
831      *p_ctype = CTYPE_OMIT;
832      /* no clino reading, so assume 0 with large sd */
833      return (real)0.0;
834   }
835   return HUGE_REAL;
836}
837
838static void
839warn_readings_differ(int msgno, real diff)
840{
841   char buf[64];
842   char *p;
843   sprintf(buf, "%.2f", deg(fabs(diff)));
844   p = strchr(buf, '.');
845   if (p) {
846      char *z = p;
847      while (*++p) {
848         if (*p != '0') z = p + 1;
849      }
850      *z = '\0';
851   }
852   compile_warning(msgno, buf);
853}
854
855static bool
856handle_comp_units(void)
857{
858   bool fNoComp = fTrue;
859   if (VAL(Comp) != HUGE_REAL) {
860      fNoComp = fFalse;
861      VAL(Comp) *= pcs->units[Q_BEARING];
862      if (VAL(Comp) < (real)0.0 || VAL(Comp) - M_PI * 2.0 > EPSILON) {
863         /* TRANSLATORS: Suspicious means something like 410 degrees or -20
864          * degrees */
865         compile_warning(/*Suspicious compass reading*/59);
866         VAL(Comp) = mod2pi(VAL(Comp));
867      }
868   }
869   if (VAL(BackComp) != HUGE_REAL) {
870      fNoComp = fFalse;
871      VAL(BackComp) *= pcs->units[Q_BACKBEARING];
872      if (VAL(BackComp) < (real)0.0 || VAL(BackComp) - M_PI * 2.0 > EPSILON) {
873         /* FIXME: different message for BackComp? */
874         compile_warning(/*Suspicious compass reading*/59);
875         VAL(BackComp) = mod2pi(VAL(BackComp));
876      }
877   }
878   return fNoComp;
879}
880
881static real
882handle_compass(real *p_var)
883{
884   real compvar = VAR(Comp);
885   real comp = VAL(Comp);
886   real backcomp = VAL(BackComp);
887   if (comp != HUGE_REAL) {
888      comp = (comp - pcs->z[Q_BEARING]) * pcs->sc[Q_BEARING];
889      comp -= pcs->z[Q_DECLINATION];
890   }
891   if (backcomp != HUGE_REAL) {
892      backcomp = (backcomp - pcs->z[Q_BACKBEARING])
893              * pcs->sc[Q_BACKBEARING];
894      backcomp -= pcs->z[Q_DECLINATION];
895      backcomp -= M_PI;
896      if (comp != HUGE_REAL) {
897         real diff = comp - backcomp;
898         real adj = fabs(diff) > M_PI ? M_PI : 0;
899         diff -= floor((diff + M_PI) / (2 * M_PI)) * 2 * M_PI;
900         if (sqrd(diff / 2.0) > compvar + VAR(Q_BACKBEARING)) {
901            /* fore and back readings differ by more than 2 sds */
902            /* TRANSLATORS: Degrees are the angular measurement where there are
903             * 360 in a full circle. */
904            warn_readings_differ(/*COMPASS reading and BACKCOMPASS reading disagree by %s degrees*/98, diff);
905         }
906         comp = (comp / compvar + backcomp / VAR(BackComp));
907         compvar = (compvar + VAR(BackComp)) / 4;
908         comp *= compvar;
909         comp += adj;
910      } else {
911         comp = backcomp;
912         compvar = VAR(BackComp);
913      }
914   }
915   *p_var = compvar;
916   return comp;
917}
918
919static int
920process_normal(prefix *fr, prefix *to, bool fToFirst,
921               clino_type ctype, clino_type backctype)
922{
923   real tape = VAL(Tape);
924   real clin = VAL(Clino);
925   real backclin = VAL(BackClino);
926
927   real dx, dy, dz;
928   real vx, vy, vz;
929#ifndef NO_COVARIANCES
930   real cxy, cyz, czx;
931#endif
932
933   bool fNoComp;
934
935   /* adjusted tape is negative -- probably the calibration is wrong */
936   if (tape < (real)0.0) {
937      /* TRANSLATE different message for topofil? */
938      compile_warning(/*Negative adjusted tape reading*/79);
939   }
940
941   fNoComp = handle_comp_units();
942
943   if (ctype == CTYPE_READING) {
944      bool range_0_180;
945      real z;
946      real diff_from_abs90;
947      clin *= pcs->units[Q_GRADIENT];
948      /* percentage scale */
949      if (pcs->f_clino_percent) clin = atan(clin);
950      /* We want to warn if there's a reading which it would be impossible
951       * to have read from the instrument (e.g. on a -90 to 90 degree scale
952       * you can't read "96" (it's probably a typo for "69").  However, the
953       * gradient reading from a topofil is typically in the range 0 to 180,
954       * with 90 being horizontal.
955       *
956       * Really we should allow the valid range to be specified, but for now
957       * we infer it from the zero error - if this is within 45 degrees of
958       * 90 then we assume the range is 0 to 180.
959       */
960      z = pcs->z[Q_GRADIENT];
961      range_0_180 = (z > M_PI_4 && z < 3*M_PI_4);
962      diff_from_abs90 = fabs(clin) - M_PI_2;
963      if (diff_from_abs90 > EPSILON) {
964         if (!range_0_180)
965            /* TRANSLATORS: Degrees are the angular measurement where there are
966             * 360 in a full circle. */
967            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
968      } else if (TSTBIT(pcs->infer, INFER_PLUMBS) &&
969                 diff_from_abs90 >= -EPSILON) {
970         ctype = CTYPE_INFERPLUMB;
971      }
972      if (range_0_180 && ctype != CTYPE_INFERPLUMB) {
973         /* FIXME: Warning message not ideal... */
974         if (clin < 0.0 || clin - M_PI > EPSILON)
975            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
976      }
977   }
978
979   if (backctype == CTYPE_READING) {
980      backclin *= pcs->units[Q_BACKGRADIENT];
981      /* percentage scale */
982      if (pcs->f_backclino_percent) backclin = atan(backclin);
983      if (ctype != CTYPE_READING) {
984         real diff_from_abs90 = fabs(backclin) - M_PI_2;
985         if (diff_from_abs90 > EPSILON) {
986            /* FIXME: different message for BackClino? */
987            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
988         } else if (TSTBIT(pcs->infer, INFER_PLUMBS) &&
989                    diff_from_abs90 >= -EPSILON) {
990            backctype = CTYPE_INFERPLUMB;
991         }
992      }
993   }
994
995   /* un-infer the plumb if the backsight was just a reading */
996   if (ctype == CTYPE_INFERPLUMB && backctype == CTYPE_READING) {
997       ctype = CTYPE_READING;
998   }
999
1000   if (ctype != CTYPE_OMIT && backctype != CTYPE_OMIT && ctype != backctype) {
1001      /* TRANSLATORS: In data with backsights, the user has tried to give a
1002       * plumb for the foresight and a clino reading for the backsight, or
1003       * something similar. */
1004      compile_error_skip(/*CLINO and BACKCLINO readings must be of the same type*/84);
1005      return 0;
1006   }
1007
1008   if (ctype == CTYPE_PLUMB || ctype == CTYPE_INFERPLUMB ||
1009       backctype == CTYPE_PLUMB || backctype == CTYPE_INFERPLUMB) {
1010      /* plumbed */
1011      if (!fNoComp) {
1012         if (ctype == CTYPE_PLUMB ||
1013             (ctype == CTYPE_INFERPLUMB && VAL(Comp) != 0.0) ||
1014             backctype == CTYPE_PLUMB ||
1015             (backctype == CTYPE_INFERPLUMB && VAL(BackComp) != 0.0)) {
1016            /* FIXME: Different message for BackComp? */
1017            /* TRANSLATORS: A "plumbed leg" is one measured using a plumbline
1018             * (a weight on a string).  So the problem here is that the leg is
1019             * vertical, so a compass reading has no meaning! */
1020            compile_warning(/*Compass reading given on plumbed leg*/21);
1021         }
1022      }
1023
1024      dx = dy = (real)0.0;
1025      if (ctype != CTYPE_OMIT) {
1026         if (backctype != CTYPE_OMIT && (clin > 0) == (backclin > 0)) {
1027            /* We've got two UPs or two DOWNs - FIXME: not ideal message */
1028            compile_error_skip(/*CLINO and BACKCLINO readings must be of the same type*/84);
1029            return 0;
1030         }
1031         dz = (clin > (real)0.0) ? tape : -tape;
1032      } else {
1033         dz = (backclin < (real)0.0) ? tape : -tape;
1034      }
1035      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1036      vz = var(Q_POS) / 3.0 + VAR(Tape);
1037#ifndef NO_COVARIANCES
1038      /* Correct values - no covariances in this case! */
1039      cxy = cyz = czx = (real)0.0;
1040#endif
1041   } else {
1042      /* Each of ctype and backctype are either CTYPE_READING/CTYPE_HORIZ
1043       * or CTYPE_OMIT */
1044      /* clino */
1045      real L2, cosG, LcosG, cosG2, sinB, cosB, dx2, dy2, dz2, v, V;
1046      if (fNoComp) {
1047         /* TRANSLATORS: Here "legs" are survey legs, i.e. measurements between
1048          * survey stations. */
1049         compile_error_skip(/*Compass reading may not be omitted except on plumbed legs*/14);
1050         return 0;
1051      }
1052      if (tape == (real)0.0) {
1053         dx = dy = dz = (real)0.0;
1054         vx = vy = vz = (real)(var(Q_POS) / 3.0); /* Position error only */
1055#ifndef NO_COVARIANCES
1056         cxy = cyz = czx = (real)0.0;
1057#endif
1058#if DEBUG_DATAIN_1
1059         printf("Zero length leg: vx = %f, vy = %f, vz = %f\n", vx, vy, vz);
1060#endif
1061      } else {
1062         real sinGcosG;
1063         /* take into account variance in LEVEL case */
1064         real var_clin = var(Q_LEVEL);
1065         real var_comp;
1066         real comp = handle_compass(&var_comp);
1067         /* ctype != CTYPE_READING is LEVEL case */
1068         if (ctype == CTYPE_READING) {
1069            clin = (clin - pcs->z[Q_GRADIENT]) * pcs->sc[Q_GRADIENT];
1070            var_clin = VAR(Clino);
1071         }
1072         if (backctype == CTYPE_READING) {
1073            backclin = (backclin - pcs->z[Q_BACKGRADIENT])
1074               * pcs->sc[Q_BACKGRADIENT];
1075            if (ctype == CTYPE_READING) {
1076               if (sqrd((clin + backclin) / 3.0) > var_clin + VAR(BackClino)) {
1077                  /* fore and back readings differ by more than 3 sds */
1078                  /* TRANSLATORS: Degrees are the angular measurement where
1079                   * there are 360 in a full circle. */
1080                  warn_readings_differ(/*CLINO reading and BACKCLINO reading disagree by %s degrees*/99, clin + backclin);
1081               }
1082               clin = (clin / var_clin - backclin / VAR(BackClino));
1083               var_clin = (var_clin + VAR(BackClino)) / 4;
1084               clin *= var_clin;
1085            } else {
1086               clin = -backclin;
1087               var_clin = VAR(BackClino);
1088            }
1089         }
1090
1091#if DEBUG_DATAIN
1092         printf("    %4.2f %4.2f %4.2f\n", tape, comp, clin);
1093#endif
1094         cosG = cos(clin);
1095         LcosG = tape * cosG;
1096         sinB = sin(comp);
1097         cosB = cos(comp);
1098#if DEBUG_DATAIN_1
1099         printf("sinB = %f, cosG = %f, LcosG = %f\n", sinB, cosG, LcosG);
1100#endif
1101         dx = LcosG * sinB;
1102         dy = LcosG * cosB;
1103         dz = tape * sin(clin);
1104/*      printf("%.2f\n",clin); */
1105#if DEBUG_DATAIN_1
1106         printf("dx = %f\ndy = %f\ndz = %f\n", dx, dy, dz);
1107#endif
1108         dx2 = dx * dx;
1109         L2 = tape * tape;
1110         V = VAR(Tape) / L2;
1111         dy2 = dy * dy;
1112         cosG2 = cosG * cosG;
1113         sinGcosG = sin(clin) * cosG;
1114         dz2 = dz * dz;
1115         v = dz2 * var_clin;
1116#ifdef NO_COVARIANCES
1117         vx = (var(Q_POS) / 3.0 + dx2 * V + dy2 * var_comp +
1118               (.5 + sinB * sinB * cosG2) * v);
1119         vy = (var(Q_POS) / 3.0 + dy2 * V + dx2 * var_comp +
1120               (.5 + cosB * cosB * cosG2) * v);
1121         if (ctype == CTYPE_OMIT && backctype == CTYPE_OMIT) {
1122            /* if no clino, assume sd=tape/sqrt(10) so 3sds = .95*tape */
1123            vz = var(Q_POS) / 3.0 + L2 * (real)0.1;
1124         } else {
1125            vz = var(Q_POS) / 3.0 + dz2 * V + L2 * cosG2 * var_clin;
1126         }
1127         /* for Surveyor87 errors: vx=vy=vz=var(Q_POS)/3.0; */
1128#else
1129         vx = var(Q_POS) / 3.0 + dx2 * V + dy2 * var_comp +
1130            (sinB * sinB * v);
1131         vy = var(Q_POS) / 3.0 + dy2 * V + dx2 * var_comp +
1132            (cosB * cosB * v);
1133         if (ctype == CTYPE_OMIT && backctype == CTYPE_OMIT) {
1134            /* if no clino, assume sd=tape/sqrt(10) so 3sds = .95*tape */
1135            vz = var(Q_POS) / 3.0 + L2 * (real)0.1;
1136         } else {
1137            vz = var(Q_POS) / 3.0 + dz2 * V + L2 * cosG2 * var_clin;
1138         }
1139         /* usual covariance formulae are fine in no clino case since
1140          * dz = 0 so value of var_clin is ignored */
1141         cxy = sinB * cosB * (VAR(Tape) * cosG2 + var_clin * dz2)
1142               - var_comp * dx * dy;
1143         czx = VAR(Tape) * sinB * sinGcosG - var_clin * dx * dz;
1144         cyz = VAR(Tape) * cosB * sinGcosG - var_clin * dy * dz;
1145#if 0
1146         printf("vx = %6.3f, vy = %6.3f, vz = %6.3f\n", vx, vy, vz);
1147         printf("cxy = %6.3f, cyz = %6.3f, czx = %6.3f\n", cxy, cyz, czx);
1148#endif
1149#endif
1150#if DEBUG_DATAIN_1
1151         printf("In DATAIN.C, vx = %f, vy = %f, vz = %f\n", vx, vy, vz);
1152#endif
1153      }
1154   }
1155#if DEBUG_DATAIN_1
1156   printf("Just before addleg, vx = %f\n", vx);
1157#endif
1158   /*printf("dx,dy,dz = %.2f %.2f %.2f\n\n", dx, dy, dz);*/
1159   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1160#ifndef NO_COVARIANCES
1161                , cyz, czx, cxy
1162#endif
1163                );
1164   return 1;
1165}
1166
1167static int
1168process_diving(prefix *fr, prefix *to, bool fToFirst, bool fDepthChange)
1169{
1170   real tape = VAL(Tape);
1171
1172   real dx, dy, dz;
1173   real vx, vy, vz;
1174#ifndef NO_COVARIANCES
1175   real cxy = 0, cyz = 0, czx = 0;
1176#endif
1177
1178   handle_comp_units();
1179
1180   /* depth gauge readings increase upwards with default calibration */
1181   if (fDepthChange) {
1182      SVX_ASSERT(VAL(FrDepth) == 0.0);
1183      dz = VAL(ToDepth) * pcs->units[Q_DEPTH] - pcs->z[Q_DEPTH];
1184      dz *= pcs->sc[Q_DEPTH];
1185   } else {
1186      dz = VAL(ToDepth) - VAL(FrDepth);
1187      dz *= pcs->units[Q_DEPTH] * pcs->sc[Q_DEPTH];
1188   }
1189
1190   /* adjusted tape is negative -- probably the calibration is wrong */
1191   if (tape < (real)0.0) {
1192      compile_warning(/*Negative adjusted tape reading*/79);
1193   }
1194
1195   /* check if tape is less than depth change */
1196   if (tape < fabs(dz)) {
1197      /* FIXME: allow margin of error based on variances? */
1198      /* TRANSLATORS: This means that the data fed in said this.
1199       *
1200       * It could be a gross error (e.g. the decimal point is missing from the
1201       * depth gauge reading) or it could just be due to random error on a near
1202       * vertical leg */
1203      compile_warning(/*Tape reading is less than change in depth*/62);
1204   }
1205
1206   if (tape == (real)0.0 && dz == 0.0) {
1207      dx = dy = dz = (real)0.0;
1208      vx = vy = vz = (real)(var(Q_POS) / 3.0); /* Position error only */
1209   } else if (VAL(Comp) == HUGE_REAL &&
1210              VAL(BackComp) == HUGE_REAL) {
1211      /* plumb */
1212      dx = dy = (real)0.0;
1213      if (dz < 0) tape = -tape;
1214      /* FIXME: Should use FrDepth sometimes... */
1215      dz = (dz * VAR(Tape) + tape * 2 * VAR(ToDepth))
1216         / (VAR(Tape) * 2 * VAR(ToDepth));
1217      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1218      /* FIXME: Should use FrDepth sometimes... */
1219      vz = var(Q_POS) / 3.0 + VAR(Tape) * 2 * VAR(ToDepth)
1220                              / (VAR(Tape) + VAR(ToDepth));
1221   } else {
1222      real L2, sinB, cosB, dz2, D2;
1223      real var_comp;
1224      real comp = handle_compass(&var_comp);
1225      sinB = sin(comp);
1226      cosB = cos(comp);
1227      L2 = tape * tape;
1228      dz2 = dz * dz;
1229      D2 = L2 - dz2;
1230      if (D2 <= (real)0.0) {
1231         /* FIXME: Should use FrDepth sometimes... */
1232         real vsum = VAR(Tape) + 2 * VAR(ToDepth);
1233         dx = dy = (real)0.0;
1234         vx = vy = var(Q_POS) / 3.0;
1235         /* FIXME: Should use FrDepth sometimes... */
1236         vz = var(Q_POS) / 3.0 + VAR(Tape) * 2 * VAR(ToDepth) / vsum;
1237         if (dz > 0) {
1238            /* FIXME: Should use FrDepth sometimes... */
1239            dz = (dz * VAR(Tape) + tape * 2 * VAR(ToDepth)) / vsum;
1240         } else {
1241            dz = (dz * VAR(Tape) - tape * 2 * VAR(ToDepth)) / vsum;
1242         }
1243      } else {
1244         real D = sqrt(D2);
1245         /* FIXME: Should use FrDepth sometimes... */
1246         real F = VAR(Tape) * L2 + 2 * VAR(ToDepth) * D2;
1247         dx = D * sinB;
1248         dy = D * cosB;
1249
1250         vx = var(Q_POS) / 3.0 +
1251            sinB * sinB * F / D2 + var_comp * dy * dy;
1252         vy = var(Q_POS) / 3.0 +
1253            cosB * cosB * F / D2 + var_comp * dx * dx;
1254         /* FIXME: Should use FrDepth sometimes... */
1255         vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1256
1257#ifndef NO_COVARIANCES
1258         cxy = sinB * cosB * (F / D2 + var_comp * D2);
1259         /* FIXME: Should use FrDepth sometimes... */
1260         cyz = -2 * VAR(ToDepth) * dy / D;
1261         czx = -2 * VAR(ToDepth) * dx / D;
1262#endif
1263      }
1264   }
1265   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1266#ifndef NO_COVARIANCES
1267                , cxy, cyz, czx
1268#endif
1269                );
1270   return 1;
1271}
1272
1273static int
1274process_cartesian(prefix *fr, prefix *to, bool fToFirst)
1275{
1276   real dx = (VAL(Dx) * pcs->units[Q_DX] - pcs->z[Q_DX]) * pcs->sc[Q_DX];
1277   real dy = (VAL(Dy) * pcs->units[Q_DY] - pcs->z[Q_DY]) * pcs->sc[Q_DY];
1278   real dz = (VAL(Dz) * pcs->units[Q_DZ] - pcs->z[Q_DZ]) * pcs->sc[Q_DZ];
1279
1280   addlegbyname(fr, to, fToFirst, dx, dy, dz, VAR(Dx), VAR(Dy), VAR(Dz)
1281#ifndef NO_COVARIANCES
1282                , 0, 0, 0
1283#endif
1284                );
1285   return 1;
1286}
1287
1288static void
1289data_cartesian(void)
1290{
1291   prefix *fr = NULL, *to = NULL;
1292
1293   bool fMulti = fFalse;
1294
1295   reading first_stn = End;
1296
1297   reading *ordering;
1298
1299   again:
1300
1301   for (ordering = pcs->ordering ; ; ordering++) {
1302      skipblanks();
1303      switch (*ordering) {
1304       case Fr:
1305         fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1306         if (first_stn == End) first_stn = Fr;
1307         break;
1308       case To:
1309         to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1310         if (first_stn == End) first_stn = To;
1311         break;
1312       case Station:
1313         fr = to;
1314         to = read_prefix(PFX_STATION);
1315         first_stn = To;
1316         break;
1317       case Dx: case Dy: case Dz:
1318         read_reading(*ordering, fFalse);
1319         break;
1320       case Ignore:
1321         skipword(); break;
1322       case IgnoreAllAndNewLine:
1323         skipline();
1324         /* fall through */
1325       case Newline:
1326         if (fr != NULL) {
1327            if (!process_cartesian(fr, to, first_stn == To))
1328               skipline();
1329         }
1330         fMulti = fTrue;
1331         while (1) {
1332            process_eol();
1333            process_bol();
1334            if (isData(ch)) break;
1335            if (!isComm(ch)) {
1336               push_back(ch);
1337               return;
1338            }
1339         }
1340         break;
1341       case IgnoreAll:
1342         skipline();
1343         /* fall through */
1344       case End:
1345         if (!fMulti) {
1346            process_cartesian(fr, to, first_stn == To);
1347            process_eol();
1348            return;
1349         }
1350         do {
1351            process_eol();
1352            process_bol();
1353         } while (isComm(ch));
1354         goto again;
1355       default: BUG("Unknown reading in ordering");
1356      }
1357   }
1358}
1359
1360static int
1361process_cylpolar(prefix *fr, prefix *to, bool fToFirst, bool fDepthChange)
1362{
1363   real tape = VAL(Tape);
1364
1365   real dx, dy, dz;
1366   real vx, vy, vz;
1367#ifndef NO_COVARIANCES
1368   real cxy = 0;
1369#endif
1370
1371   handle_comp_units();
1372
1373   /* depth gauge readings increase upwards with default calibration */
1374   if (fDepthChange) {
1375      SVX_ASSERT(VAL(FrDepth) == 0.0);
1376      dz = VAL(ToDepth) * pcs->units[Q_DEPTH] - pcs->z[Q_DEPTH];
1377      dz *= pcs->sc[Q_DEPTH];
1378   } else {
1379      dz = VAL(ToDepth) - VAL(FrDepth);
1380      dz *= pcs->units[Q_DEPTH] * pcs->sc[Q_DEPTH];
1381   }
1382
1383   /* adjusted tape is negative -- probably the calibration is wrong */
1384   if (tape < (real)0.0) {
1385      compile_warning(/*Negative adjusted tape reading*/79);
1386   }
1387
1388   if (VAL(Comp) == HUGE_REAL && VAL(BackComp) == HUGE_REAL) {
1389      /* plumb */
1390      dx = dy = (real)0.0;
1391      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1392      /* FIXME: Should use FrDepth sometimes... */
1393      vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1394   } else {
1395      real sinB, cosB;
1396      real var_comp;
1397      real comp = handle_compass(&var_comp);
1398      sinB = sin(comp);
1399      cosB = cos(comp);
1400
1401      dx = tape * sinB;
1402      dy = tape * cosB;
1403
1404      vx = var(Q_POS) / 3.0 +
1405         VAR(Tape) * sinB * sinB + var_comp * dy * dy;
1406      vy = var(Q_POS) / 3.0 +
1407         VAR(Tape) * cosB * cosB + var_comp * dx * dx;
1408      /* FIXME: Should use FrDepth sometimes... */
1409      vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1410
1411#ifndef NO_COVARIANCES
1412      cxy = (VAR(Tape) - var_comp * tape * tape) * sinB * cosB;
1413#endif
1414   }
1415   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1416#ifndef NO_COVARIANCES
1417                , cxy, 0, 0
1418#endif
1419                );
1420   return 1;
1421}
1422
1423/* Process tape/compass/clino, diving, and cylpolar styles of survey data
1424 * Also handles topofil (fromcount/tocount or count) in place of tape */
1425static void
1426data_normal(void)
1427{
1428   prefix *fr = NULL, *to = NULL;
1429   reading first_stn = End;
1430
1431   bool fTopofil = fFalse, fMulti = fFalse;
1432   bool fRev;
1433   clino_type ctype, backctype;
1434   bool fDepthChange;
1435   unsigned long compass_dat_flags = 0;
1436
1437   reading *ordering;
1438
1439   VAL(Tape) = 0;
1440   VAL(Comp) = VAL(BackComp) = HUGE_REAL;
1441   VAL(FrCount) = VAL(ToCount) = 0;
1442   VAL(FrDepth) = VAL(ToDepth) = 0;
1443   VAL(Left) = VAL(Right) = VAL(Up) = VAL(Down) = HUGE_REAL;
1444
1445   fRev = fFalse;
1446   ctype = backctype = CTYPE_OMIT;
1447   fDepthChange = fFalse;
1448
1449   /* ordering may omit clino reading, so set up default here */
1450   /* this is also used if clino reading is the omit character */
1451   VAL(Clino) = VAL(BackClino) = 0;
1452
1453   again:
1454
1455   /* We clear these flags in the normal course of events, but if there's an
1456    * error in a reading, we might not, so make sure it has been cleared here.
1457    */
1458   pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1459   for (ordering = pcs->ordering; ; ordering++) {
1460      skipblanks();
1461      switch (*ordering) {
1462       case Fr:
1463          fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT|PFX_ANON);
1464          if (first_stn == End) first_stn = Fr;
1465          break;
1466       case To:
1467          to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT|PFX_ANON);
1468          if (first_stn == End) first_stn = To;
1469          break;
1470       case Station:
1471          fr = to;
1472          to = read_prefix(PFX_STATION);
1473          first_stn = To;
1474          break;
1475       case Dir: {
1476          typedef enum {
1477             DIR_NULL=-1, DIR_FORE, DIR_BACK
1478          } dir_tok;
1479          static sztok dir_tab[] = {
1480             {"B",     DIR_BACK},
1481             {"F",     DIR_FORE},
1482          };
1483          dir_tok tok;
1484          get_token();
1485          tok = match_tok(dir_tab, TABSIZE(dir_tab));
1486          switch (tok) {
1487           case DIR_FORE:
1488             break;
1489           case DIR_BACK:
1490             fRev = fTrue;
1491             break;
1492           default:
1493             file.lpos += strlen(buffer);
1494             compile_error_skip(-/*Found “%s”, expecting “F” or “B”*/131,
1495                                buffer);
1496             process_eol();
1497             return;
1498          }
1499          break;
1500       }
1501       case Tape:
1502          read_reading(Tape, fFalse);
1503          if (VAL(Tape) < (real)0.0)
1504             compile_warning(-/*Negative tape reading*/60);
1505          break;
1506       case Count:
1507          VAL(FrCount) = VAL(ToCount);
1508          read_reading(ToCount, fFalse);
1509          fTopofil = fTrue;
1510          break;
1511       case FrCount:
1512          read_reading(FrCount, fFalse);
1513          break;
1514       case ToCount:
1515          read_reading(ToCount, fFalse);
1516          fTopofil = fTrue;
1517          break;
1518       case Comp: case BackComp:
1519          read_bearing_or_omit(*ordering);
1520          break;
1521       case Clino: case BackClino: {
1522          reading r = *ordering;
1523          clino_type * p_ctype = (r == Clino ? &ctype : &backctype);
1524          read_reading(r, fTrue);
1525          if (VAL(r) == HUGE_REAL) {
1526             VAL(r) = handle_plumb(p_ctype);
1527             if (VAL(r) != HUGE_REAL) break;
1528             compile_error_token(-/*Expecting numeric field, found “%s”*/9);
1529             skipline();
1530             process_eol();
1531             return;
1532          }
1533          *p_ctype = CTYPE_READING;
1534          break;
1535       }
1536       case FrDepth: case ToDepth:
1537          read_reading(*ordering, fFalse);
1538          break;
1539       case Depth:
1540          VAL(FrDepth) = VAL(ToDepth);
1541          read_reading(ToDepth, fFalse);
1542          break;
1543       case DepthChange:
1544          fDepthChange = fTrue;
1545          VAL(FrDepth) = 0;
1546          read_reading(ToDepth, fFalse);
1547          break;
1548       case CompassDATComp:
1549          read_bearing_or_omit(Comp);
1550          if (is_compass_NaN(VAL(Comp))) VAL(Comp) = HUGE_REAL;
1551          break;
1552       case CompassDATBackComp:
1553          read_bearing_or_omit(BackComp);
1554          if (is_compass_NaN(VAL(BackComp))) VAL(BackComp) = HUGE_REAL;
1555          break;
1556       case CompassDATClino: case CompassDATBackClino: {
1557          reading r;
1558          clino_type * p_ctype;
1559          if (*ordering == CompassDATClino) {
1560             r = Clino;
1561             p_ctype = &ctype;
1562          } else {
1563             r = BackClino;
1564             p_ctype = &backctype;
1565          }
1566          read_reading(r, fFalse);
1567          if (is_compass_NaN(VAL(r))) {
1568             VAL(r) = HUGE_REAL;
1569             *p_ctype = CTYPE_OMIT;
1570          } else {
1571             *p_ctype = CTYPE_READING;
1572          }
1573          break;
1574       }
1575       case CompassDATLeft: case CompassDATRight:
1576       case CompassDATUp: case CompassDATDown: {
1577          /* FIXME: need to actually make use of these entries! */
1578          reading actual = Left + (*ordering - CompassDATLeft);
1579          read_reading(actual, fFalse);
1580          if (VAL(actual) < 0) VAL(actual) = HUGE_REAL;
1581          break;
1582       }
1583       case CompassDATFlags:
1584          if (ch == '#') {
1585             nextch();
1586             if (ch == '|') {
1587                filepos fp;
1588                get_pos(&fp);
1589                nextch();
1590                while (ch >= 'A' && ch <= 'Z') {
1591                   compass_dat_flags |= BIT(ch - 'A');
1592                   /* FIXME: we currently understand:
1593                    *   L (exclude from length)
1594                    *   X (exclude data)
1595                    * but should also handle at least some of:
1596                    *   C (no adjustment)
1597                    *   P (no plot)
1598                    */
1599                   nextch();
1600                }
1601                if (ch == '#') {
1602                   nextch();
1603                } else {
1604                   compass_dat_flags = 0;
1605                   set_pos(&fp);
1606                   push_back('|');
1607                   ch = '#';
1608                }
1609             } else {
1610                push_back(ch);
1611                ch = '#';
1612             }
1613          }
1614          break;
1615       case Ignore:
1616          skipword(); break;
1617       case IgnoreAllAndNewLine:
1618          skipline();
1619          /* fall through */
1620       case Newline:
1621          if (fr != NULL) {
1622             int r;
1623             int save_flags;
1624             int implicit_splay;
1625             if (fTopofil)
1626                VAL(Tape) = VAL(ToCount) - VAL(FrCount);
1627             /* Note: frdepth == todepth test works regardless of fDepthChange
1628              * (frdepth always zero, todepth is change of depth) and also
1629              * works for STYLE_NORMAL (both remain 0) */
1630             if (TSTBIT(pcs->infer, INFER_EQUATES) &&
1631                 VAL(Tape) == (real)0.0 && VAL(FrDepth) == VAL(ToDepth)) {
1632                process_equate(fr, to);
1633                goto inferred_equate;
1634             }
1635             if (fRev) {
1636                prefix *t = fr;
1637                fr = to;
1638                to = t;
1639             }
1640             if (fTopofil) {
1641                VAL(Tape) *= pcs->units[Q_COUNT] * pcs->sc[Q_COUNT];
1642             } else {
1643                VAL(Tape) *= pcs->units[Q_LENGTH];
1644                VAL(Tape) -= pcs->z[Q_LENGTH];
1645                VAL(Tape) *= pcs->sc[Q_LENGTH];
1646             }
1647             implicit_splay = TSTBIT(pcs->flags, FLAGS_IMPLICIT_SPLAY);
1648             pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1649             save_flags = pcs->flags;
1650             if (implicit_splay) {
1651                pcs->flags |= BIT(FLAGS_SPLAY);
1652             }
1653             switch (pcs->style) {
1654              case STYLE_NORMAL:
1655                r = process_normal(fr, to, (first_stn == To) ^ fRev,
1656                                   ctype, backctype);
1657                break;
1658              case STYLE_DIVING:
1659                r = process_diving(fr, to, (first_stn == To) ^ fRev,
1660                                   fDepthChange);
1661                break;
1662              case STYLE_CYLPOLAR:
1663                r = process_cylpolar(fr, to, (first_stn == To) ^ fRev,
1664                                     fDepthChange);
1665                break;
1666              default:
1667                r = 0; /* avoid warning */
1668                BUG("bad style");
1669             }
1670             pcs->flags = save_flags;
1671             if (!r) skipline();
1672
1673             /* Swap fr and to back to how they were for next line */
1674             if (fRev) {
1675                prefix *t = fr;
1676                fr = to;
1677                to = t;
1678             }
1679          }
1680
1681          fRev = fFalse;
1682          ctype = backctype = CTYPE_OMIT;
1683          fDepthChange = fFalse;
1684
1685          /* ordering may omit clino reading, so set up default here */
1686          /* this is also used if clino reading is the omit character */
1687          VAL(Clino) = VAL(BackClino) = 0;
1688
1689          inferred_equate:
1690
1691          fMulti = fTrue;
1692          while (1) {
1693              process_eol();
1694              process_bol();
1695              if (isData(ch)) break;
1696              if (!isComm(ch)) {
1697                 push_back(ch);
1698                 return;
1699              }
1700          }
1701          break;
1702       case IgnoreAll:
1703          skipline();
1704          /* fall through */
1705       case End:
1706          if (!fMulti) {
1707             int save_flags;
1708             int implicit_splay;
1709             /* Compass ignore flag is 'X' */
1710             if ((compass_dat_flags & BIT('X' - 'A'))) {
1711                process_eol();
1712                return;
1713             }
1714             if (fRev) {
1715                prefix *t = fr;
1716                fr = to;
1717                to = t;
1718             }
1719             if (fTopofil) VAL(Tape) = VAL(ToCount) - VAL(FrCount);
1720             /* Note: frdepth == todepth test works regardless of fDepthChange
1721              * (frdepth always zero, todepth is change of depth) and also
1722              * works for STYLE_NORMAL (both remain 0) */
1723             if (TSTBIT(pcs->infer, INFER_EQUATES) &&
1724                 VAL(Tape) == (real)0.0 && VAL(FrDepth) == VAL(ToDepth)) {
1725                process_equate(fr, to);
1726                process_eol();
1727                return;
1728             }
1729             if (fTopofil) {
1730                VAL(Tape) *= pcs->units[Q_COUNT] * pcs->sc[Q_COUNT];
1731             } else {
1732                VAL(Tape) *= pcs->units[Q_LENGTH];
1733                VAL(Tape) -= pcs->z[Q_LENGTH];
1734                VAL(Tape) *= pcs->sc[Q_LENGTH];
1735             }
1736             implicit_splay = TSTBIT(pcs->flags, FLAGS_IMPLICIT_SPLAY);
1737             pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1738             save_flags = pcs->flags;
1739             if (implicit_splay) {
1740                pcs->flags |= BIT(FLAGS_SPLAY);
1741             }
1742             if ((compass_dat_flags & BIT('L' - 'A'))) {
1743                /* 'L' means "exclude from length" - map this to Survex's
1744                 * FLAGS_DUPLICATE. */
1745                pcs->flags |= BIT(FLAGS_DUPLICATE);
1746             }
1747             switch (pcs->style) {
1748              case STYLE_NORMAL:
1749                process_normal(fr, to, (first_stn == To) ^ fRev,
1750                               ctype, backctype);
1751                break;
1752              case STYLE_DIVING:
1753                process_diving(fr, to, (first_stn == To) ^ fRev,
1754                               fDepthChange);
1755                break;
1756              case STYLE_CYLPOLAR:
1757                process_cylpolar(fr, to, (first_stn == To) ^ fRev,
1758                                 fDepthChange);
1759                break;
1760              default:
1761                BUG("bad style");
1762             }
1763             pcs->flags = save_flags;
1764
1765             process_eol();
1766             return;
1767          }
1768          do {
1769             process_eol();
1770             process_bol();
1771          } while (isComm(ch));
1772          goto again;
1773       default:
1774          BUG("Unknown reading in ordering");
1775      }
1776   }
1777}
1778
1779static int
1780process_lrud(prefix *stn)
1781{
1782   SVX_ASSERT(next_lrud);
1783   lrud * xsect = osnew(lrud);
1784   xsect->stn = stn;
1785   xsect->l = (VAL(Left) * pcs->units[Q_LEFT] - pcs->z[Q_LEFT]) * pcs->sc[Q_LEFT];
1786   xsect->r = (VAL(Right) * pcs->units[Q_RIGHT] - pcs->z[Q_RIGHT]) * pcs->sc[Q_RIGHT];
1787   xsect->u = (VAL(Up) * pcs->units[Q_UP] - pcs->z[Q_UP]) * pcs->sc[Q_UP];
1788   xsect->d = (VAL(Down) * pcs->units[Q_DOWN] - pcs->z[Q_DOWN]) * pcs->sc[Q_DOWN];
1789   xsect->meta = pcs->meta;
1790   if (pcs->meta) ++pcs->meta->ref_count;
1791   xsect->next = NULL;
1792   *next_lrud = xsect;
1793   next_lrud = &(xsect->next);
1794
1795   return 1;
1796}
1797
1798static void
1799data_passage(void)
1800{
1801   prefix *stn = NULL;
1802   reading *ordering;
1803
1804   for (ordering = pcs->ordering ; ; ordering++) {
1805      skipblanks();
1806      switch (*ordering) {
1807       case Station:
1808         stn = read_prefix(PFX_STATION);
1809         break;
1810       case Left: case Right: case Up: case Down:
1811         read_reading(*ordering, fTrue);
1812         if (VAL(*ordering) == HUGE_REAL) {
1813            if (!isOmit(ch)) {
1814               compile_error_token(-/*Expecting numeric field, found “%s”*/9);
1815            } else {
1816               nextch();
1817            }
1818            VAL(*ordering) = -1;
1819         }
1820         break;
1821       case Ignore:
1822         skipword(); break;
1823       case IgnoreAll:
1824         skipline();
1825         /* fall through */
1826       case End: {
1827         process_lrud(stn);
1828         process_eol();
1829         return;
1830       }
1831       default: BUG("Unknown reading in ordering");
1832      }
1833   }
1834}
1835
1836static int
1837process_nosurvey(prefix *fr, prefix *to, bool fToFirst)
1838{
1839   nosurveylink *link;
1840
1841   /* Suppress "unused fixed point" warnings for these stations */
1842   fr->sflags |= BIT(SFLAGS_USED);
1843   to->sflags |= BIT(SFLAGS_USED);
1844
1845   /* add to linked list which is dealt with after network is solved */
1846   link = osnew(nosurveylink);
1847   if (fToFirst) {
1848      link->to = StnFromPfx(to);
1849      link->fr = StnFromPfx(fr);
1850   } else {
1851      link->fr = StnFromPfx(fr);
1852      link->to = StnFromPfx(to);
1853   }
1854   link->flags = pcs->flags | (STYLE_NOSURVEY << FLAGS_STYLE_BIT0);
1855   link->meta = pcs->meta;
1856   if (pcs->meta) ++pcs->meta->ref_count;
1857   link->next = nosurveyhead;
1858   nosurveyhead = link;
1859   return 1;
1860}
1861
1862static void
1863data_nosurvey(void)
1864{
1865   prefix *fr = NULL, *to = NULL;
1866
1867   bool fMulti = fFalse;
1868
1869   reading first_stn = End;
1870
1871   reading *ordering;
1872
1873   again:
1874
1875   for (ordering = pcs->ordering ; ; ordering++) {
1876      skipblanks();
1877      switch (*ordering) {
1878       case Fr:
1879          fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1880          if (first_stn == End) first_stn = Fr;
1881          break;
1882       case To:
1883          to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1884          if (first_stn == End) first_stn = To;
1885          break;
1886       case Station:
1887          fr = to;
1888          to = read_prefix(PFX_STATION);
1889          first_stn = To;
1890          break;
1891       case Ignore:
1892         skipword(); break;
1893       case IgnoreAllAndNewLine:
1894         skipline();
1895         /* fall through */
1896       case Newline:
1897         if (fr != NULL) {
1898            if (!process_nosurvey(fr, to, first_stn == To))
1899               skipline();
1900         }
1901         if (ordering[1] == End) {
1902            do {
1903               process_eol();
1904               process_bol();
1905            } while (isComm(ch));
1906            if (!isData(ch)) {
1907               push_back(ch);
1908               return;
1909            }
1910            goto again;
1911         }
1912         fMulti = fTrue;
1913         while (1) {
1914            process_eol();
1915            process_bol();
1916            if (isData(ch)) break;
1917            if (!isComm(ch)) {
1918               push_back(ch);
1919               return;
1920            }
1921         }
1922         break;
1923       case IgnoreAll:
1924         skipline();
1925         /* fall through */
1926       case End:
1927         if (!fMulti) {
1928            (void)process_nosurvey(fr, to, first_stn == To);
1929            process_eol();
1930            return;
1931         }
1932         do {
1933            process_eol();
1934            process_bol();
1935         } while (isComm(ch));
1936         goto again;
1937       default: BUG("Unknown reading in ordering");
1938      }
1939   }
1940}
1941
1942/* totally ignore a line of survey data */
1943static void
1944data_ignore(void)
1945{
1946   skipline();
1947   process_eol();
1948}
Note: See TracBrowser for help on using the repository browser.