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