src/datain.c

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