source: git/src/datain.c @ 85c0078

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

doc/TODO.htm,lib/,src/,tests/: In cavern, use the currently set units
when outputting measurements in warnings, errors, and the stats at
the end of the run.

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