//
// Simple tool to unwrap HP 53131A time interval phase data.
//
// - Detects and fixes phase jumps caused by cycle wrapping.
// - User-specified nominal period value (default 100 ns).
// - Input is raw time interval (phase) measurement, and
// -   also accepts HP 53131A/53132A counter TI data format.
// - Assumes well-behaved TI (phase) measurements; where
// -   gradual phase drift is assumed to be frequency offset, and
// -   full-period phase jumps are assumed to be cycle wraps, but
// -   phase jumps neither large nor small are flagged for review.
// - Recreates continuous phase record as output.
//
// 24-Oct-2006 /tvb
//

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

double period;

#define TI_MAX  (1.00 * period)  // maximum expected TI value
#define TI_MIN (-1.10 * period)  // minimum expected TI value
#define TI_WRAP  (0.8 * period)  // threshold for detecting cycle wrap
#define TI_JUMP  (0.3 * period)  // threshold for detecting phase jump

double hp53131 (char *line);

void main (int argc, char *argv[])
{
    char Line[1000], *format;
    double ti, ti0, dt, adt, cycles, units, scale;
    struct { long lines, samples, races, wraps; } n;

    period = (argc > 1) ? atof(argv[1]) : 0.100e-6;
    format = (argc > 2) ? argv[2] : "%.3le";
    scale = (argc > 3) ? atof(argv[3]) : 1.0;
    fprintf(stderr, "period %lg, format %s, scale %lg\n",
            period, format, scale);

    n.lines = n.samples = n.races = n.wraps = 0;
    cycles = 0.0;

    while (fgets(Line, sizeof(Line), stdin) != NULL) {
        n.lines += 1;
        if ((n.lines % 10000) == 0) {
            fprintf(stderr, "\r%ld\r", n.lines);
        }
        if ((units = hp53131(Line)) == 0.0) {
            continue;
        }
        sscanf(Line, "%lf", &ti);
        ti *= units;

        // Integrity check - make sure all TI samples are within one period.

        if (ti > TI_MAX) {
            fprintf(stderr, "%ld: WARNING - TI too high? %lg\n", n.lines, ti);
            // continue;
        }
        if (ti < TI_MIN)  {
            fprintf(stderr, "%ld: WARNING - TI too low? %lg\n", n.lines, ti);
            // continue;
        }
        n.samples += 1;

        // Adjust for TI counter channel race.

        if (ti < 0) {
            ti += period;
            n.races += 1;
        }
        dt = (n.lines > 1) ? (ti - ti0) : 0;
        adt = fabs(dt);

        // Integrity check - a big phase jump can look like a cycle wrap.

        if (adt > TI_JUMP && adt < TI_WRAP) {
            fprintf(stderr, "%ld: WARNING - big phase jump? %lg\n", n.lines, dt);
        }

        // Adjust phase for cycle wrap.

        if (adt > TI_WRAP) {
            if (dt < 0) {
                cycles += 1;
            } else {
                cycles -= 1;
            }
            n.wraps += 1;
        }

        printf(format, (cycles * period + ti) * scale);
        putchar('\n');
        ti0 = ti;
    }

    fprintf(stderr, "%ld lines, %ld samples, %ld races, %ld wraps, %lg net cycles\n",
            n.lines, n.samples, n.races, n.wraps, cycles);
}

// Handle HP/Agilent 53131A/53132A RS-232 output format.

double hp53131 (char *line)
{
    char *p;

    // Ignore user comments
    if (line[0] == '#') {
        return 0;
    }

    // Ignore HP 53131A statistics
    if (strstr(line, ": ") != NULL) {
        return 0;
    }

    // Remove embedded commas (US version HP 53131A)
    while ((p = strchr(line, ',')) != NULL) {
        strcpy(p, p + 1);
    }

    // Handle HP microsecond and second units
    if (strstr(line, " us") != NULL) {
        return 1e-6;

    } else if (strstr(line, " s") != NULL) {
        return 1.0;
    }

    return 1.0;
}

#if 0

unwrap1 < rapco_raw.txt > phase.dat
unwrap1 < rapco_raw.txt 1e-7 %.4lf 1e6 > rapco.dat
unwrap1 1e-7 "%.1f ns" 1e9 < rapco_raw.txt > rapco.txt
unwrap1 1e-7 "%.0f ps" 1e12 < rapco_raw.txt > rapco.txt

#endif