src/datain.c

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