source: git/src/datain.c @ 4dc9899

line_contentsstereotravis-osx
Last change on this file since 4dc9899 was 4dc9899, checked in by Olly Betts <olly@…>, 5 years ago

src/datain.c: Warn if fore and back measurements differ by more than
3 SDs, rather than just by 2 SDs.

  • Property mode set to 100644
File size: 50.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)
842{
843   char buf[64];
844   char *p;
845   sprintf(buf, "%.2f", deg(fabs(diff)));
846   p = strchr(buf, '.');
847   if (p) {
848      char *z = p;
849      while (*++p) {
850         if (*p != '0') z = p + 1;
851      }
852      *z = '\0';
853   }
854   compile_warning(msgno, buf);
855}
856
857static bool
858handle_comp_units(void)
859{
860   bool fNoComp = fTrue;
861   if (VAL(Comp) != HUGE_REAL) {
862      fNoComp = fFalse;
863      VAL(Comp) *= pcs->units[Q_BEARING];
864      if (VAL(Comp) < (real)0.0 || VAL(Comp) - M_PI * 2.0 > EPSILON) {
865         /* TRANSLATORS: Suspicious means something like 410 degrees or -20
866          * degrees */
867         compile_warning(/*Suspicious compass reading*/59);
868         VAL(Comp) = mod2pi(VAL(Comp));
869      }
870   }
871   if (VAL(BackComp) != HUGE_REAL) {
872      fNoComp = fFalse;
873      VAL(BackComp) *= pcs->units[Q_BACKBEARING];
874      if (VAL(BackComp) < (real)0.0 || VAL(BackComp) - M_PI * 2.0 > EPSILON) {
875         /* FIXME: different message for BackComp? */
876         compile_warning(/*Suspicious compass reading*/59);
877         VAL(BackComp) = mod2pi(VAL(BackComp));
878      }
879   }
880   return fNoComp;
881}
882
883static real
884handle_compass(real *p_var)
885{
886   real compvar = VAR(Comp);
887   real comp = VAL(Comp);
888   real backcomp = VAL(BackComp);
889   if (comp != HUGE_REAL) {
890      comp = (comp - pcs->z[Q_BEARING]) * pcs->sc[Q_BEARING];
891      comp -= pcs->z[Q_DECLINATION];
892   }
893   if (backcomp != HUGE_REAL) {
894      backcomp = (backcomp - pcs->z[Q_BACKBEARING])
895              * pcs->sc[Q_BACKBEARING];
896      backcomp -= pcs->z[Q_DECLINATION];
897      backcomp -= M_PI;
898      if (comp != HUGE_REAL) {
899         real diff = comp - backcomp;
900         real adj = fabs(diff) > M_PI ? M_PI : 0;
901         diff -= floor((diff + M_PI) / (2 * M_PI)) * 2 * M_PI;
902         if (sqrd(diff / 3.0) > compvar + VAR(Q_BACKBEARING)) {
903            /* fore and back readings differ by more than 3 sds */
904            /* TRANSLATORS: Degrees are the angular measurement where there are
905             * 360 in a full circle. */
906            warn_readings_differ(/*COMPASS reading and BACKCOMPASS reading disagree by %s degrees*/98, diff);
907         }
908         comp = (comp / compvar + backcomp / VAR(BackComp));
909         compvar = (compvar + VAR(BackComp)) / 4;
910         comp *= compvar;
911         comp += adj;
912      } else {
913         comp = backcomp;
914         compvar = VAR(BackComp);
915      }
916   }
917   *p_var = compvar;
918   return comp;
919}
920
921static int
922process_normal(prefix *fr, prefix *to, bool fToFirst,
923               clino_type ctype, clino_type backctype)
924{
925   real tape = VAL(Tape);
926   real clin = VAL(Clino);
927   real backclin = VAL(BackClino);
928
929   real dx, dy, dz;
930   real vx, vy, vz;
931#ifndef NO_COVARIANCES
932   real cxy, cyz, czx;
933#endif
934
935   bool fNoComp;
936
937   /* adjusted tape is negative -- probably the calibration is wrong */
938   if (tape < (real)0.0) {
939      /* TRANSLATE different message for topofil? */
940      compile_warning(/*Negative adjusted tape reading*/79);
941   }
942
943   fNoComp = handle_comp_units();
944
945   if (ctype == CTYPE_READING) {
946      bool range_0_180;
947      real z;
948      real diff_from_abs90;
949      clin *= pcs->units[Q_GRADIENT];
950      /* percentage scale */
951      if (pcs->f_clino_percent) clin = atan(clin);
952      /* We want to warn if there's a reading which it would be impossible
953       * to have read from the instrument (e.g. on a -90 to 90 degree scale
954       * you can't read "96" (it's probably a typo for "69").  However, the
955       * gradient reading from a topofil is typically in the range 0 to 180,
956       * with 90 being horizontal.
957       *
958       * Really we should allow the valid range to be specified, but for now
959       * we infer it from the zero error - if this is within 45 degrees of
960       * 90 then we assume the range is 0 to 180.
961       */
962      z = pcs->z[Q_GRADIENT];
963      range_0_180 = (z > M_PI_4 && z < 3*M_PI_4);
964      diff_from_abs90 = fabs(clin) - M_PI_2;
965      if (diff_from_abs90 > EPSILON) {
966         if (!range_0_180)
967            /* TRANSLATORS: Degrees are the angular measurement where there are
968             * 360 in a full circle. */
969            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
970      } else if (TSTBIT(pcs->infer, INFER_PLUMBS) &&
971                 diff_from_abs90 >= -EPSILON) {
972         ctype = CTYPE_INFERPLUMB;
973      }
974      if (range_0_180 && ctype != CTYPE_INFERPLUMB) {
975         /* FIXME: Warning message not ideal... */
976         if (clin < 0.0 || clin - M_PI > EPSILON)
977            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
978      }
979   }
980
981   if (backctype == CTYPE_READING) {
982      backclin *= pcs->units[Q_BACKGRADIENT];
983      /* percentage scale */
984      if (pcs->f_backclino_percent) backclin = atan(backclin);
985      if (ctype != CTYPE_READING) {
986         real diff_from_abs90 = fabs(backclin) - M_PI_2;
987         if (diff_from_abs90 > EPSILON) {
988            /* FIXME: different message for BackClino? */
989            compile_warning(/*Clino reading over 90 degrees (absolute value)*/51);
990         } else if (TSTBIT(pcs->infer, INFER_PLUMBS) &&
991                    diff_from_abs90 >= -EPSILON) {
992            backctype = CTYPE_INFERPLUMB;
993         }
994      }
995   }
996
997   /* un-infer the plumb if the backsight was just a reading */
998   if (ctype == CTYPE_INFERPLUMB && backctype == CTYPE_READING) {
999       ctype = CTYPE_READING;
1000   }
1001
1002   if (ctype != CTYPE_OMIT && backctype != CTYPE_OMIT && ctype != backctype) {
1003      /* TRANSLATORS: In data with backsights, the user has tried to give a
1004       * plumb for the foresight and a clino reading for the backsight, or
1005       * something similar. */
1006      compile_error_skip(/*CLINO and BACKCLINO readings must be of the same type*/84);
1007      return 0;
1008   }
1009
1010   if (ctype == CTYPE_PLUMB || ctype == CTYPE_INFERPLUMB ||
1011       backctype == CTYPE_PLUMB || backctype == CTYPE_INFERPLUMB) {
1012      /* plumbed */
1013      if (!fNoComp) {
1014         if (ctype == CTYPE_PLUMB ||
1015             (ctype == CTYPE_INFERPLUMB && VAL(Comp) != 0.0) ||
1016             backctype == CTYPE_PLUMB ||
1017             (backctype == CTYPE_INFERPLUMB && VAL(BackComp) != 0.0)) {
1018            /* FIXME: Different message for BackComp? */
1019            /* TRANSLATORS: A "plumbed leg" is one measured using a plumbline
1020             * (a weight on a string).  So the problem here is that the leg is
1021             * vertical, so a compass reading has no meaning! */
1022            compile_warning(/*Compass reading given on plumbed leg*/21);
1023         }
1024      }
1025
1026      dx = dy = (real)0.0;
1027      if (ctype != CTYPE_OMIT) {
1028         if (backctype != CTYPE_OMIT && (clin > 0) == (backclin > 0)) {
1029            /* We've got two UPs or two DOWNs - FIXME: not ideal message */
1030            compile_error_skip(/*CLINO and BACKCLINO readings must be of the same type*/84);
1031            return 0;
1032         }
1033         dz = (clin > (real)0.0) ? tape : -tape;
1034      } else {
1035         dz = (backclin < (real)0.0) ? tape : -tape;
1036      }
1037      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1038      vz = var(Q_POS) / 3.0 + VAR(Tape);
1039#ifndef NO_COVARIANCES
1040      /* Correct values - no covariances in this case! */
1041      cxy = cyz = czx = (real)0.0;
1042#endif
1043   } else {
1044      /* Each of ctype and backctype are either CTYPE_READING/CTYPE_HORIZ
1045       * or CTYPE_OMIT */
1046      /* clino */
1047      real L2, cosG, LcosG, cosG2, sinB, cosB, dx2, dy2, dz2, v, V;
1048      if (fNoComp) {
1049         /* TRANSLATORS: Here "legs" are survey legs, i.e. measurements between
1050          * survey stations. */
1051         compile_error_skip(/*Compass reading may not be omitted except on plumbed legs*/14);
1052         return 0;
1053      }
1054      if (tape == (real)0.0) {
1055         dx = dy = dz = (real)0.0;
1056         vx = vy = vz = (real)(var(Q_POS) / 3.0); /* Position error only */
1057#ifndef NO_COVARIANCES
1058         cxy = cyz = czx = (real)0.0;
1059#endif
1060#if DEBUG_DATAIN_1
1061         printf("Zero length leg: vx = %f, vy = %f, vz = %f\n", vx, vy, vz);
1062#endif
1063      } else {
1064         real sinGcosG;
1065         /* take into account variance in LEVEL case */
1066         real var_clin = var(Q_LEVEL);
1067         real var_comp;
1068         real comp = handle_compass(&var_comp);
1069         /* ctype != CTYPE_READING is LEVEL case */
1070         if (ctype == CTYPE_READING) {
1071            clin = (clin - pcs->z[Q_GRADIENT]) * pcs->sc[Q_GRADIENT];
1072            var_clin = VAR(Clino);
1073         }
1074         if (backctype == CTYPE_READING) {
1075            backclin = (backclin - pcs->z[Q_BACKGRADIENT])
1076               * pcs->sc[Q_BACKGRADIENT];
1077            if (ctype == CTYPE_READING) {
1078               if (sqrd((clin + backclin) / 3.0) > var_clin + VAR(BackClino)) {
1079                  /* fore and back readings differ by more than 3 sds */
1080                  /* TRANSLATORS: Degrees are the angular measurement where
1081                   * there are 360 in a full circle. */
1082                  warn_readings_differ(/*CLINO reading and BACKCLINO reading disagree by %s degrees*/99, clin + backclin);
1083               }
1084               clin = (clin / var_clin - backclin / VAR(BackClino));
1085               var_clin = (var_clin + VAR(BackClino)) / 4;
1086               clin *= var_clin;
1087            } else {
1088               clin = -backclin;
1089               var_clin = VAR(BackClino);
1090            }
1091         }
1092
1093#if DEBUG_DATAIN
1094         printf("    %4.2f %4.2f %4.2f\n", tape, comp, clin);
1095#endif
1096         cosG = cos(clin);
1097         LcosG = tape * cosG;
1098         sinB = sin(comp);
1099         cosB = cos(comp);
1100#if DEBUG_DATAIN_1
1101         printf("sinB = %f, cosG = %f, LcosG = %f\n", sinB, cosG, LcosG);
1102#endif
1103         dx = LcosG * sinB;
1104         dy = LcosG * cosB;
1105         dz = tape * sin(clin);
1106/*      printf("%.2f\n",clin); */
1107#if DEBUG_DATAIN_1
1108         printf("dx = %f\ndy = %f\ndz = %f\n", dx, dy, dz);
1109#endif
1110         dx2 = dx * dx;
1111         L2 = tape * tape;
1112         V = VAR(Tape) / L2;
1113         dy2 = dy * dy;
1114         cosG2 = cosG * cosG;
1115         sinGcosG = sin(clin) * cosG;
1116         dz2 = dz * dz;
1117         v = dz2 * var_clin;
1118#ifdef NO_COVARIANCES
1119         vx = (var(Q_POS) / 3.0 + dx2 * V + dy2 * var_comp +
1120               (.5 + sinB * sinB * cosG2) * v);
1121         vy = (var(Q_POS) / 3.0 + dy2 * V + dx2 * var_comp +
1122               (.5 + cosB * cosB * cosG2) * v);
1123         if (ctype == CTYPE_OMIT && backctype == CTYPE_OMIT) {
1124            /* if no clino, assume sd=tape/sqrt(10) so 3sds = .95*tape */
1125            vz = var(Q_POS) / 3.0 + L2 * (real)0.1;
1126         } else {
1127            vz = var(Q_POS) / 3.0 + dz2 * V + L2 * cosG2 * var_clin;
1128         }
1129         /* for Surveyor87 errors: vx=vy=vz=var(Q_POS)/3.0; */
1130#else
1131         vx = var(Q_POS) / 3.0 + dx2 * V + dy2 * var_comp +
1132            (sinB * sinB * v);
1133         vy = var(Q_POS) / 3.0 + dy2 * V + dx2 * var_comp +
1134            (cosB * cosB * v);
1135         if (ctype == CTYPE_OMIT && backctype == CTYPE_OMIT) {
1136            /* if no clino, assume sd=tape/sqrt(10) so 3sds = .95*tape */
1137            vz = var(Q_POS) / 3.0 + L2 * (real)0.1;
1138         } else {
1139            vz = var(Q_POS) / 3.0 + dz2 * V + L2 * cosG2 * var_clin;
1140         }
1141         /* usual covariance formulae are fine in no clino case since
1142          * dz = 0 so value of var_clin is ignored */
1143         cxy = sinB * cosB * (VAR(Tape) * cosG2 + var_clin * dz2)
1144               - var_comp * dx * dy;
1145         czx = VAR(Tape) * sinB * sinGcosG - var_clin * dx * dz;
1146         cyz = VAR(Tape) * cosB * sinGcosG - var_clin * dy * dz;
1147#if 0
1148         printf("vx = %6.3f, vy = %6.3f, vz = %6.3f\n", vx, vy, vz);
1149         printf("cxy = %6.3f, cyz = %6.3f, czx = %6.3f\n", cxy, cyz, czx);
1150#endif
1151#endif
1152#if DEBUG_DATAIN_1
1153         printf("In DATAIN.C, vx = %f, vy = %f, vz = %f\n", vx, vy, vz);
1154#endif
1155      }
1156   }
1157#if DEBUG_DATAIN_1
1158   printf("Just before addleg, vx = %f\n", vx);
1159#endif
1160   /*printf("dx,dy,dz = %.2f %.2f %.2f\n\n", dx, dy, dz);*/
1161   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1162#ifndef NO_COVARIANCES
1163                , cyz, czx, cxy
1164#endif
1165                );
1166   return 1;
1167}
1168
1169static int
1170process_diving(prefix *fr, prefix *to, bool fToFirst, bool fDepthChange)
1171{
1172   real tape = VAL(Tape);
1173
1174   real dx, dy, dz;
1175   real vx, vy, vz;
1176#ifndef NO_COVARIANCES
1177   real cxy = 0, cyz = 0, czx = 0;
1178#endif
1179
1180   handle_comp_units();
1181
1182   /* depth gauge readings increase upwards with default calibration */
1183   if (fDepthChange) {
1184      SVX_ASSERT(VAL(FrDepth) == 0.0);
1185      dz = VAL(ToDepth) * pcs->units[Q_DEPTH] - pcs->z[Q_DEPTH];
1186      dz *= pcs->sc[Q_DEPTH];
1187   } else {
1188      dz = VAL(ToDepth) - VAL(FrDepth);
1189      dz *= pcs->units[Q_DEPTH] * pcs->sc[Q_DEPTH];
1190   }
1191
1192   /* adjusted tape is negative -- probably the calibration is wrong */
1193   if (tape < (real)0.0) {
1194      compile_warning(/*Negative adjusted tape reading*/79);
1195   }
1196
1197   /* check if tape is less than depth change */
1198   if (tape < fabs(dz)) {
1199      /* FIXME: allow margin of error based on variances? */
1200      /* TRANSLATORS: This means that the data fed in said this.
1201       *
1202       * It could be a gross error (e.g. the decimal point is missing from the
1203       * depth gauge reading) or it could just be due to random error on a near
1204       * vertical leg */
1205      compile_warning(/*Tape reading is less than change in depth*/62);
1206   }
1207
1208   if (tape == (real)0.0 && dz == 0.0) {
1209      dx = dy = dz = (real)0.0;
1210      vx = vy = vz = (real)(var(Q_POS) / 3.0); /* Position error only */
1211   } else if (VAL(Comp) == HUGE_REAL &&
1212              VAL(BackComp) == HUGE_REAL) {
1213      /* plumb */
1214      dx = dy = (real)0.0;
1215      if (dz < 0) tape = -tape;
1216      /* FIXME: Should use FrDepth sometimes... */
1217      dz = (dz * VAR(Tape) + tape * 2 * VAR(ToDepth))
1218         / (VAR(Tape) * 2 * VAR(ToDepth));
1219      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1220      /* FIXME: Should use FrDepth sometimes... */
1221      vz = var(Q_POS) / 3.0 + VAR(Tape) * 2 * VAR(ToDepth)
1222                              / (VAR(Tape) + VAR(ToDepth));
1223   } else {
1224      real L2, sinB, cosB, dz2, D2;
1225      real var_comp;
1226      real comp = handle_compass(&var_comp);
1227      sinB = sin(comp);
1228      cosB = cos(comp);
1229      L2 = tape * tape;
1230      dz2 = dz * dz;
1231      D2 = L2 - dz2;
1232      if (D2 <= (real)0.0) {
1233         /* FIXME: Should use FrDepth sometimes... */
1234         real vsum = VAR(Tape) + 2 * VAR(ToDepth);
1235         dx = dy = (real)0.0;
1236         vx = vy = var(Q_POS) / 3.0;
1237         /* FIXME: Should use FrDepth sometimes... */
1238         vz = var(Q_POS) / 3.0 + VAR(Tape) * 2 * VAR(ToDepth) / vsum;
1239         if (dz > 0) {
1240            /* FIXME: Should use FrDepth sometimes... */
1241            dz = (dz * VAR(Tape) + tape * 2 * VAR(ToDepth)) / vsum;
1242         } else {
1243            dz = (dz * VAR(Tape) - tape * 2 * VAR(ToDepth)) / vsum;
1244         }
1245      } else {
1246         real D = sqrt(D2);
1247         /* FIXME: Should use FrDepth sometimes... */
1248         real F = VAR(Tape) * L2 + 2 * VAR(ToDepth) * D2;
1249         dx = D * sinB;
1250         dy = D * cosB;
1251
1252         vx = var(Q_POS) / 3.0 +
1253            sinB * sinB * F / D2 + var_comp * dy * dy;
1254         vy = var(Q_POS) / 3.0 +
1255            cosB * cosB * F / D2 + var_comp * dx * dx;
1256         /* FIXME: Should use FrDepth sometimes... */
1257         vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1258
1259#ifndef NO_COVARIANCES
1260         cxy = sinB * cosB * (F / D2 + var_comp * D2);
1261         /* FIXME: Should use FrDepth sometimes... */
1262         cyz = -2 * VAR(ToDepth) * dy / D;
1263         czx = -2 * VAR(ToDepth) * dx / D;
1264#endif
1265      }
1266   }
1267   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1268#ifndef NO_COVARIANCES
1269                , cxy, cyz, czx
1270#endif
1271                );
1272   return 1;
1273}
1274
1275static int
1276process_cartesian(prefix *fr, prefix *to, bool fToFirst)
1277{
1278   real dx = (VAL(Dx) * pcs->units[Q_DX] - pcs->z[Q_DX]) * pcs->sc[Q_DX];
1279   real dy = (VAL(Dy) * pcs->units[Q_DY] - pcs->z[Q_DY]) * pcs->sc[Q_DY];
1280   real dz = (VAL(Dz) * pcs->units[Q_DZ] - pcs->z[Q_DZ]) * pcs->sc[Q_DZ];
1281
1282   addlegbyname(fr, to, fToFirst, dx, dy, dz, VAR(Dx), VAR(Dy), VAR(Dz)
1283#ifndef NO_COVARIANCES
1284                , 0, 0, 0
1285#endif
1286                );
1287   return 1;
1288}
1289
1290static void
1291data_cartesian(void)
1292{
1293   prefix *fr = NULL, *to = NULL;
1294
1295   bool fMulti = fFalse;
1296
1297   reading first_stn = End;
1298
1299   reading *ordering;
1300
1301   again:
1302
1303   for (ordering = pcs->ordering ; ; ordering++) {
1304      skipblanks();
1305      switch (*ordering) {
1306       case Fr:
1307         fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1308         if (first_stn == End) first_stn = Fr;
1309         break;
1310       case To:
1311         to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1312         if (first_stn == End) first_stn = To;
1313         break;
1314       case Station:
1315         fr = to;
1316         to = read_prefix(PFX_STATION);
1317         first_stn = To;
1318         break;
1319       case Dx: case Dy: case Dz:
1320         read_reading(*ordering, fFalse);
1321         break;
1322       case Ignore:
1323         skipword(); break;
1324       case IgnoreAllAndNewLine:
1325         skipline();
1326         /* fall through */
1327       case Newline:
1328         if (fr != NULL) {
1329            if (!process_cartesian(fr, to, first_stn == To))
1330               skipline();
1331         }
1332         fMulti = fTrue;
1333         while (1) {
1334            process_eol();
1335            process_bol();
1336            if (isData(ch)) break;
1337            if (!isComm(ch)) {
1338               push_back(ch);
1339               return;
1340            }
1341         }
1342         break;
1343       case IgnoreAll:
1344         skipline();
1345         /* fall through */
1346       case End:
1347         if (!fMulti) {
1348            process_cartesian(fr, to, first_stn == To);
1349            process_eol();
1350            return;
1351         }
1352         do {
1353            process_eol();
1354            process_bol();
1355         } while (isComm(ch));
1356         goto again;
1357       default: BUG("Unknown reading in ordering");
1358      }
1359   }
1360}
1361
1362static int
1363process_cylpolar(prefix *fr, prefix *to, bool fToFirst, bool fDepthChange)
1364{
1365   real tape = VAL(Tape);
1366
1367   real dx, dy, dz;
1368   real vx, vy, vz;
1369#ifndef NO_COVARIANCES
1370   real cxy = 0;
1371#endif
1372
1373   handle_comp_units();
1374
1375   /* depth gauge readings increase upwards with default calibration */
1376   if (fDepthChange) {
1377      SVX_ASSERT(VAL(FrDepth) == 0.0);
1378      dz = VAL(ToDepth) * pcs->units[Q_DEPTH] - pcs->z[Q_DEPTH];
1379      dz *= pcs->sc[Q_DEPTH];
1380   } else {
1381      dz = VAL(ToDepth) - VAL(FrDepth);
1382      dz *= pcs->units[Q_DEPTH] * pcs->sc[Q_DEPTH];
1383   }
1384
1385   /* adjusted tape is negative -- probably the calibration is wrong */
1386   if (tape < (real)0.0) {
1387      compile_warning(/*Negative adjusted tape reading*/79);
1388   }
1389
1390   if (VAL(Comp) == HUGE_REAL && VAL(BackComp) == HUGE_REAL) {
1391      /* plumb */
1392      dx = dy = (real)0.0;
1393      vx = vy = var(Q_POS) / 3.0 + dz * dz * var(Q_PLUMB);
1394      /* FIXME: Should use FrDepth sometimes... */
1395      vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1396   } else {
1397      real sinB, cosB;
1398      real var_comp;
1399      real comp = handle_compass(&var_comp);
1400      sinB = sin(comp);
1401      cosB = cos(comp);
1402
1403      dx = tape * sinB;
1404      dy = tape * cosB;
1405
1406      vx = var(Q_POS) / 3.0 +
1407         VAR(Tape) * sinB * sinB + var_comp * dy * dy;
1408      vy = var(Q_POS) / 3.0 +
1409         VAR(Tape) * cosB * cosB + var_comp * dx * dx;
1410      /* FIXME: Should use FrDepth sometimes... */
1411      vz = var(Q_POS) / 3.0 + 2 * VAR(ToDepth);
1412
1413#ifndef NO_COVARIANCES
1414      cxy = (VAR(Tape) - var_comp * tape * tape) * sinB * cosB;
1415#endif
1416   }
1417   addlegbyname(fr, to, fToFirst, dx, dy, dz, vx, vy, vz
1418#ifndef NO_COVARIANCES
1419                , cxy, 0, 0
1420#endif
1421                );
1422   return 1;
1423}
1424
1425/* Process tape/compass/clino, diving, and cylpolar styles of survey data
1426 * Also handles topofil (fromcount/tocount or count) in place of tape */
1427static void
1428data_normal(void)
1429{
1430   prefix *fr = NULL, *to = NULL;
1431   reading first_stn = End;
1432
1433   bool fTopofil = fFalse, fMulti = fFalse;
1434   bool fRev;
1435   clino_type ctype, backctype;
1436   bool fDepthChange;
1437   unsigned long compass_dat_flags = 0;
1438
1439   reading *ordering;
1440
1441   VAL(Tape) = 0;
1442   VAL(Comp) = VAL(BackComp) = HUGE_REAL;
1443   VAL(FrCount) = VAL(ToCount) = 0;
1444   VAL(FrDepth) = VAL(ToDepth) = 0;
1445   VAL(Left) = VAL(Right) = VAL(Up) = VAL(Down) = HUGE_REAL;
1446
1447   fRev = fFalse;
1448   ctype = backctype = CTYPE_OMIT;
1449   fDepthChange = fFalse;
1450
1451   /* ordering may omit clino reading, so set up default here */
1452   /* this is also used if clino reading is the omit character */
1453   VAL(Clino) = VAL(BackClino) = 0;
1454
1455   again:
1456
1457   /* We clear these flags in the normal course of events, but if there's an
1458    * error in a reading, we might not, so make sure it has been cleared here.
1459    */
1460   pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1461   for (ordering = pcs->ordering; ; ordering++) {
1462      skipblanks();
1463      switch (*ordering) {
1464       case Fr:
1465          fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT|PFX_ANON);
1466          if (first_stn == End) first_stn = Fr;
1467          break;
1468       case To:
1469          to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT|PFX_ANON);
1470          if (first_stn == End) first_stn = To;
1471          break;
1472       case Station:
1473          fr = to;
1474          to = read_prefix(PFX_STATION);
1475          first_stn = To;
1476          break;
1477       case Dir: {
1478          typedef enum {
1479             DIR_NULL=-1, DIR_FORE, DIR_BACK
1480          } dir_tok;
1481          static sztok dir_tab[] = {
1482             {"B",     DIR_BACK},
1483             {"F",     DIR_FORE},
1484          };
1485          dir_tok tok;
1486          get_token();
1487          tok = match_tok(dir_tab, TABSIZE(dir_tab));
1488          switch (tok) {
1489           case DIR_FORE:
1490             break;
1491           case DIR_BACK:
1492             fRev = fTrue;
1493             break;
1494           default:
1495             file.lpos += strlen(buffer);
1496             compile_error_skip(-/*Found “%s”, expecting “F” or “B”*/131,
1497                                buffer);
1498             process_eol();
1499             return;
1500          }
1501          break;
1502       }
1503       case Tape:
1504          read_reading(Tape, fFalse);
1505          if (VAL(Tape) < (real)0.0)
1506             compile_warning(-/*Negative tape reading*/60);
1507          break;
1508       case Count:
1509          VAL(FrCount) = VAL(ToCount);
1510          read_reading(ToCount, fFalse);
1511          fTopofil = fTrue;
1512          break;
1513       case FrCount:
1514          read_reading(FrCount, fFalse);
1515          break;
1516       case ToCount:
1517          read_reading(ToCount, fFalse);
1518          fTopofil = fTrue;
1519          break;
1520       case Comp: case BackComp:
1521          read_bearing_or_omit(*ordering);
1522          break;
1523       case Clino: case BackClino: {
1524          reading r = *ordering;
1525          clino_type * p_ctype = (r == Clino ? &ctype : &backctype);
1526          read_reading(r, fTrue);
1527          if (VAL(r) == HUGE_REAL) {
1528             VAL(r) = handle_plumb(p_ctype);
1529             if (VAL(r) != HUGE_REAL) break;
1530             compile_error_token(-/*Expecting numeric field, found “%s”*/9);
1531             skipline();
1532             process_eol();
1533             return;
1534          }
1535          *p_ctype = CTYPE_READING;
1536          break;
1537       }
1538       case FrDepth: case ToDepth:
1539          read_reading(*ordering, fFalse);
1540          break;
1541       case Depth:
1542          VAL(FrDepth) = VAL(ToDepth);
1543          read_reading(ToDepth, fFalse);
1544          break;
1545       case DepthChange:
1546          fDepthChange = fTrue;
1547          VAL(FrDepth) = 0;
1548          read_reading(ToDepth, fFalse);
1549          break;
1550       case CompassDATComp:
1551          read_bearing_or_omit(Comp);
1552          if (is_compass_NaN(VAL(Comp))) VAL(Comp) = HUGE_REAL;
1553          break;
1554       case CompassDATBackComp:
1555          read_bearing_or_omit(BackComp);
1556          if (is_compass_NaN(VAL(BackComp))) VAL(BackComp) = HUGE_REAL;
1557          break;
1558       case CompassDATClino: case CompassDATBackClino: {
1559          reading r;
1560          clino_type * p_ctype;
1561          if (*ordering == CompassDATClino) {
1562             r = Clino;
1563             p_ctype = &ctype;
1564          } else {
1565             r = BackClino;
1566             p_ctype = &backctype;
1567          }
1568          read_reading(r, fFalse);
1569          if (is_compass_NaN(VAL(r))) {
1570             VAL(r) = HUGE_REAL;
1571             *p_ctype = CTYPE_OMIT;
1572          } else {
1573             *p_ctype = CTYPE_READING;
1574          }
1575          break;
1576       }
1577       case CompassDATLeft: case CompassDATRight:
1578       case CompassDATUp: case CompassDATDown: {
1579          /* FIXME: need to actually make use of these entries! */
1580          reading actual = Left + (*ordering - CompassDATLeft);
1581          read_reading(actual, fFalse);
1582          if (VAL(actual) < 0) VAL(actual) = HUGE_REAL;
1583          break;
1584       }
1585       case CompassDATFlags:
1586          if (ch == '#') {
1587             nextch();
1588             if (ch == '|') {
1589                filepos fp;
1590                get_pos(&fp);
1591                nextch();
1592                while (ch >= 'A' && ch <= 'Z') {
1593                   compass_dat_flags |= BIT(ch - 'A');
1594                   /* FIXME: we currently understand:
1595                    *   L (exclude from length)
1596                    *   X (exclude data)
1597                    * but should also handle at least some of:
1598                    *   C (no adjustment)
1599                    *   P (no plot)
1600                    */
1601                   nextch();
1602                }
1603                if (ch == '#') {
1604                   nextch();
1605                } else {
1606                   compass_dat_flags = 0;
1607                   set_pos(&fp);
1608                   push_back('|');
1609                   ch = '#';
1610                }
1611             } else {
1612                push_back(ch);
1613                ch = '#';
1614             }
1615          }
1616          break;
1617       case Ignore:
1618          skipword(); break;
1619       case IgnoreAllAndNewLine:
1620          skipline();
1621          /* fall through */
1622       case Newline:
1623          if (fr != NULL) {
1624             int r;
1625             int save_flags;
1626             int implicit_splay;
1627             if (fTopofil)
1628                VAL(Tape) = VAL(ToCount) - VAL(FrCount);
1629             /* Note: frdepth == todepth test works regardless of fDepthChange
1630              * (frdepth always zero, todepth is change of depth) and also
1631              * works for STYLE_NORMAL (both remain 0) */
1632             if (TSTBIT(pcs->infer, INFER_EQUATES) &&
1633                 VAL(Tape) == (real)0.0 && VAL(FrDepth) == VAL(ToDepth)) {
1634                process_equate(fr, to);
1635                goto inferred_equate;
1636             }
1637             if (fRev) {
1638                prefix *t = fr;
1639                fr = to;
1640                to = t;
1641             }
1642             if (fTopofil) {
1643                VAL(Tape) *= pcs->units[Q_COUNT] * pcs->sc[Q_COUNT];
1644             } else {
1645                VAL(Tape) *= pcs->units[Q_LENGTH];
1646                VAL(Tape) -= pcs->z[Q_LENGTH];
1647                VAL(Tape) *= pcs->sc[Q_LENGTH];
1648             }
1649             implicit_splay = TSTBIT(pcs->flags, FLAGS_IMPLICIT_SPLAY);
1650             pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1651             save_flags = pcs->flags;
1652             if (implicit_splay) {
1653                pcs->flags |= BIT(FLAGS_SPLAY);
1654             }
1655             switch (pcs->style) {
1656              case STYLE_NORMAL:
1657                r = process_normal(fr, to, (first_stn == To) ^ fRev,
1658                                   ctype, backctype);
1659                break;
1660              case STYLE_DIVING:
1661                r = process_diving(fr, to, (first_stn == To) ^ fRev,
1662                                   fDepthChange);
1663                break;
1664              case STYLE_CYLPOLAR:
1665                r = process_cylpolar(fr, to, (first_stn == To) ^ fRev,
1666                                     fDepthChange);
1667                break;
1668              default:
1669                r = 0; /* avoid warning */
1670                BUG("bad style");
1671             }
1672             pcs->flags = save_flags;
1673             if (!r) skipline();
1674
1675             /* Swap fr and to back to how they were for next line */
1676             if (fRev) {
1677                prefix *t = fr;
1678                fr = to;
1679                to = t;
1680             }
1681          }
1682
1683          fRev = fFalse;
1684          ctype = backctype = CTYPE_OMIT;
1685          fDepthChange = fFalse;
1686
1687          /* ordering may omit clino reading, so set up default here */
1688          /* this is also used if clino reading is the omit character */
1689          VAL(Clino) = VAL(BackClino) = 0;
1690
1691          inferred_equate:
1692
1693          fMulti = fTrue;
1694          while (1) {
1695              process_eol();
1696              process_bol();
1697              if (isData(ch)) break;
1698              if (!isComm(ch)) {
1699                 push_back(ch);
1700                 return;
1701              }
1702          }
1703          break;
1704       case IgnoreAll:
1705          skipline();
1706          /* fall through */
1707       case End:
1708          if (!fMulti) {
1709             int save_flags;
1710             int implicit_splay;
1711             /* Compass ignore flag is 'X' */
1712             if ((compass_dat_flags & BIT('X' - 'A'))) {
1713                process_eol();
1714                return;
1715             }
1716             if (fRev) {
1717                prefix *t = fr;
1718                fr = to;
1719                to = t;
1720             }
1721             if (fTopofil) VAL(Tape) = VAL(ToCount) - VAL(FrCount);
1722             /* Note: frdepth == todepth test works regardless of fDepthChange
1723              * (frdepth always zero, todepth is change of depth) and also
1724              * works for STYLE_NORMAL (both remain 0) */
1725             if (TSTBIT(pcs->infer, INFER_EQUATES) &&
1726                 VAL(Tape) == (real)0.0 && VAL(FrDepth) == VAL(ToDepth)) {
1727                process_equate(fr, to);
1728                process_eol();
1729                return;
1730             }
1731             if (fTopofil) {
1732                VAL(Tape) *= pcs->units[Q_COUNT] * pcs->sc[Q_COUNT];
1733             } else {
1734                VAL(Tape) *= pcs->units[Q_LENGTH];
1735                VAL(Tape) -= pcs->z[Q_LENGTH];
1736                VAL(Tape) *= pcs->sc[Q_LENGTH];
1737             }
1738             implicit_splay = TSTBIT(pcs->flags, FLAGS_IMPLICIT_SPLAY);
1739             pcs->flags &= ~(BIT(FLAGS_ANON_ONE_END) | BIT(FLAGS_IMPLICIT_SPLAY));
1740             save_flags = pcs->flags;
1741             if (implicit_splay) {
1742                pcs->flags |= BIT(FLAGS_SPLAY);
1743             }
1744             if ((compass_dat_flags & BIT('L' - 'A'))) {
1745                /* 'L' means "exclude from length" - map this to Survex's
1746                 * FLAGS_DUPLICATE. */
1747                pcs->flags |= BIT(FLAGS_DUPLICATE);
1748             }
1749             switch (pcs->style) {
1750              case STYLE_NORMAL:
1751                process_normal(fr, to, (first_stn == To) ^ fRev,
1752                               ctype, backctype);
1753                break;
1754              case STYLE_DIVING:
1755                process_diving(fr, to, (first_stn == To) ^ fRev,
1756                               fDepthChange);
1757                break;
1758              case STYLE_CYLPOLAR:
1759                process_cylpolar(fr, to, (first_stn == To) ^ fRev,
1760                                 fDepthChange);
1761                break;
1762              default:
1763                BUG("bad style");
1764             }
1765             pcs->flags = save_flags;
1766
1767             process_eol();
1768             return;
1769          }
1770          do {
1771             process_eol();
1772             process_bol();
1773          } while (isComm(ch));
1774          goto again;
1775       default:
1776          BUG("Unknown reading in ordering");
1777      }
1778   }
1779}
1780
1781static int
1782process_lrud(prefix *stn)
1783{
1784   SVX_ASSERT(next_lrud);
1785   lrud * xsect = osnew(lrud);
1786   xsect->stn = stn;
1787   xsect->l = (VAL(Left) * pcs->units[Q_LEFT] - pcs->z[Q_LEFT]) * pcs->sc[Q_LEFT];
1788   xsect->r = (VAL(Right) * pcs->units[Q_RIGHT] - pcs->z[Q_RIGHT]) * pcs->sc[Q_RIGHT];
1789   xsect->u = (VAL(Up) * pcs->units[Q_UP] - pcs->z[Q_UP]) * pcs->sc[Q_UP];
1790   xsect->d = (VAL(Down) * pcs->units[Q_DOWN] - pcs->z[Q_DOWN]) * pcs->sc[Q_DOWN];
1791   xsect->meta = pcs->meta;
1792   if (pcs->meta) ++pcs->meta->ref_count;
1793   xsect->next = NULL;
1794   *next_lrud = xsect;
1795   next_lrud = &(xsect->next);
1796
1797   return 1;
1798}
1799
1800static void
1801data_passage(void)
1802{
1803   prefix *stn = NULL;
1804   reading *ordering;
1805
1806   for (ordering = pcs->ordering ; ; ordering++) {
1807      skipblanks();
1808      switch (*ordering) {
1809       case Station:
1810         stn = read_prefix(PFX_STATION);
1811         break;
1812       case Left: case Right: case Up: case Down:
1813         read_reading(*ordering, fTrue);
1814         if (VAL(*ordering) == HUGE_REAL) {
1815            if (!isOmit(ch)) {
1816               compile_error_token(-/*Expecting numeric field, found “%s”*/9);
1817            } else {
1818               nextch();
1819            }
1820            VAL(*ordering) = -1;
1821         }
1822         break;
1823       case Ignore:
1824         skipword(); break;
1825       case IgnoreAll:
1826         skipline();
1827         /* fall through */
1828       case End: {
1829         process_lrud(stn);
1830         process_eol();
1831         return;
1832       }
1833       default: BUG("Unknown reading in ordering");
1834      }
1835   }
1836}
1837
1838static int
1839process_nosurvey(prefix *fr, prefix *to, bool fToFirst)
1840{
1841   nosurveylink *link;
1842
1843   /* Suppress "unused fixed point" warnings for these stations */
1844   fr->sflags |= BIT(SFLAGS_USED);
1845   to->sflags |= BIT(SFLAGS_USED);
1846
1847   /* add to linked list which is dealt with after network is solved */
1848   link = osnew(nosurveylink);
1849   if (fToFirst) {
1850      link->to = StnFromPfx(to);
1851      link->fr = StnFromPfx(fr);
1852   } else {
1853      link->fr = StnFromPfx(fr);
1854      link->to = StnFromPfx(to);
1855   }
1856   link->flags = pcs->flags | (STYLE_NOSURVEY << FLAGS_STYLE_BIT0);
1857   link->meta = pcs->meta;
1858   if (pcs->meta) ++pcs->meta->ref_count;
1859   link->next = nosurveyhead;
1860   nosurveyhead = link;
1861   return 1;
1862}
1863
1864static void
1865data_nosurvey(void)
1866{
1867   prefix *fr = NULL, *to = NULL;
1868
1869   bool fMulti = fFalse;
1870
1871   reading first_stn = End;
1872
1873   reading *ordering;
1874
1875   again:
1876
1877   for (ordering = pcs->ordering ; ; ordering++) {
1878      skipblanks();
1879      switch (*ordering) {
1880       case Fr:
1881          fr = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1882          if (first_stn == End) first_stn = Fr;
1883          break;
1884       case To:
1885          to = read_prefix(PFX_STATION|PFX_ALLOW_ROOT);
1886          if (first_stn == End) first_stn = To;
1887          break;
1888       case Station:
1889          fr = to;
1890          to = read_prefix(PFX_STATION);
1891          first_stn = To;
1892          break;
1893       case Ignore:
1894         skipword(); break;
1895       case IgnoreAllAndNewLine:
1896         skipline();
1897         /* fall through */
1898       case Newline:
1899         if (fr != NULL) {
1900            if (!process_nosurvey(fr, to, first_stn == To))
1901               skipline();
1902         }
1903         if (ordering[1] == End) {
1904            do {
1905               process_eol();
1906               process_bol();
1907            } while (isComm(ch));
1908            if (!isData(ch)) {
1909               push_back(ch);
1910               return;
1911            }
1912            goto again;
1913         }
1914         fMulti = fTrue;
1915         while (1) {
1916            process_eol();
1917            process_bol();
1918            if (isData(ch)) break;
1919            if (!isComm(ch)) {
1920               push_back(ch);
1921               return;
1922            }
1923         }
1924         break;
1925       case IgnoreAll:
1926         skipline();
1927         /* fall through */
1928       case End:
1929         if (!fMulti) {
1930            (void)process_nosurvey(fr, to, first_stn == To);
1931            process_eol();
1932            return;
1933         }
1934         do {
1935            process_eol();
1936            process_bol();
1937         } while (isComm(ch));
1938         goto again;
1939       default: BUG("Unknown reading in ordering");
1940      }
1941   }
1942}
1943
1944/* totally ignore a line of survey data */
1945static void
1946data_ignore(void)
1947{
1948   skipline();
1949   process_eol();
1950}
Note: See TracBrowser for help on using the repository browser.