//
// adev3 - Demo program to compute several Allan deviation statistics.
//
//    ADEV - Allan deviation (classic or overlapping)
//    MDEV - Modified Allan deviation
//    HDEV - Hadamard deviation (classic or overlapping)
//
// See Usage[] below for instructions.
//
// This was written primarily to explore the accuracy and speed of the
// algorithms. Calculations match the results from Stable32 (the gold
// standard of commercial clock statistics programs). Input is a file
// of phase (time interval error) measurements. Output is one line or
// a table of Allan deviations (with sample counts in parenthesis).
//
// Build (Win32) with cl /Ox (highly optimized).
//
// 28-Dec-2008 (tvb) Tom Van Baak
//

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

char Usage[] =
    "Tool to compute ADEV, MDEV, and HDEV from raw phase data\n"
    "\n"
    "Usage: adev3 [options] table-style < input\n"
    "\n"
    "Options:\n"
    "    /f:#   - input field number (default 1)\n"
    "    /s:#   - input scaling factor (default 1.0)\n"
    "    /t0:#  - tau0 value, sampling interval (default 1 second)\n"
    "    /xl    - display table in Excel CSV format\n"
    "    /dp:#  - number of decimal places\n"
    "    /ld:#  - number of digits for term counts\n"
    "    /b2b   - use traditional back-to-back non-overlapping calculations\n"
    "    /q     - quiet (no headings or timings)\n"
    "\n"
    " Style of Allan deviation table:\n"
    "     #      - calculate and display just one tau\n"
    "     i,j    - all tau from i to j\n"
    "     i,j,k  - tau from i to j in k increments\n"
    "     all    - all tau\n"
    "     many   - many tau (72 tau per decade, equal density)\n"
    "     m#     - many tau (# tau per decade, equal density)\n"
    "     two    - octave (1, 2, 4, 8, etc)\n"
    "     ten    - decade (1, 10, 100, 1000, etc)\n"
    "     124    - decade (1, 2, 4, 10, 20, 40, etc)\n"
    "     125    - decade (1, 2, 5, 10, 20, 50, etc)\n"
    "     nine   - 9 steps per decade (1,2,3...9,10,20,30...90, etc)\n"
    "\n"
    "Examples\n"
    "    adev3 < phase.dat\n"
    "    adev3 two < phase.dat\n"
    "    adev3 /xl nine < phase.dat\n"
    ;

struct {
    int dp;         // decimal places
    int field;      // data field number
    int ld;         // integer digits
    int ovlp;       // overlapping
    int quiet;      // quiet mode (no headings)
    double scale;   // input data scale factor
    char *style;    // table style
    long tau0;      // input data interval
    int xl;         // Excel csv output
} opt;

#define FPRINTF !opt.quiet && fprintf

#define MIN_SAMPLES 3  // sanity check for too large tau or too little data.

//
// Calculate Allan deviation (ADEV) from given phase array.
//  if ovlp == 0 then use classic calculation (back-to-back intervals).
//  if ovlp == 1 then use fully overlapping intervals.
//

double calc_adev (double data[], long count, long tau, long *terms, int ovlp)
{
    long i, n, stride;
    double sum, v;

    stride = ovlp ? 1 : tau;
    sum = n = 0;
    for (i = 0; (i + 2*tau) < count; i += stride) {
        v = data[i + 2*tau] - 2 * data[i + tau] + data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 2.0;

    *terms = n;
    return (n < MIN_SAMPLES) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Calculate Hadamard deviation (HDEV) from given phase array.
//  if ovlp == 0 then use classic calculation (back-to-back intervals).
//  if ovlp == 1 then use fully overlapping intervals.
//

double calc_hdev (double data[], long count, long tau, long *terms, int ovlp)
{
    long i, n, stride;
    double sum, v;

    stride = ovlp ? 1 : tau;
    sum = n = 0;
    for (i = 0; (i + 3*tau) < count; i += stride) {
        v = data[i + 3*tau] - 3 * data[i + 2*tau] + 3 * data[i + tau] - data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 6.0;

    *terms = n;
    return (n < MIN_SAMPLES) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Calculate Modified Allan deviation (MDEV) from given phase array.
// - This is a much faster (unrolled, single-loop) version.
// - Note long tau*tau can exceed 2^31 so tau is promoted to double.
//

double calc_mdev (double data[], long count, long tau, long *terms)
{
    long i, n;
    double sum, v;

    sum = v = n = 0;
    for (i = 0; (i + 2*tau) < count && i < tau; i += 1) {
        v += data[i + 2*tau] - 2 * data[i + tau] + data[i];
    }
    sum += v * v;
    n += 1;

    for (i = 0; (i + 3*tau) < count; i += 1) {
        v += data[i + 3*tau] - 3*data[i + 2*tau] + 3*data[i + tau] - data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 2.0 * (double)tau * (double)tau;

    *terms = n;
    return (n < MIN_SAMPLES) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Perform all three Allan calculations and display results.
//

int stop;           // hack to stop open-ended table calculations
double *phase;      // array of data read from input file
long count;
double *read_data (FILE *f, long *count);

void sigma_tau (long tau)
{
    struct { double adev, mdev, hdev; } a;
    struct {  long  adev, mdev, hdev; } n;
    static long clock_start;

    if (stop) return;

    // First time -- read phase data from stdin the first time.

    if (phase == NULL) {
        FPRINTF(stderr, "(tau0 %ld)\n", opt.tau0);
        FPRINTF(stderr, "(scale %lg)\n", opt.scale);
        phase = read_data(stdin, &count);
        FPRINTF(stderr, "(read %ld samples)\n", count);
        FPRINTF(stderr, "(table style = '%s')\n", opt.style);
        FPRINTF(stderr, "(output %d decimal places)\n", opt.dp);
        FPRINTF(stderr, "(term counts %d digits)\n", opt.ld);
        FPRINTF(stderr, "(overlapping %s)\n", opt.ovlp ? "yes" : "no");
        FPRINTF(stderr, "Displaying Tau, ADEV, MDEV, HDEV (each with term counts)\n");
        clock_start = clock();
    }

    a.adev = calc_adev(phase, count, tau, &n.adev, opt.ovlp) / opt.tau0;
    a.mdev = calc_mdev(phase, count, tau, &n.mdev) / opt.tau0;
    a.hdev = calc_hdev(phase, count, tau, &n.hdev, opt.ovlp) / opt.tau0;

    stop = (a.adev == 0.0) && (a.mdev == 0.0) && (a.hdev == 0.0);

    // Last time -- display calculation timing.

    if (stop) {
        long clock_end;
        double elapsed;

        clock_end = clock();
        elapsed = (double)(clock_end - clock_start) / CLOCKS_PER_SEC;
        FPRINTF(stderr, "(elapsed time %.2lf seconds)\n", elapsed);
        free(phase);
        return;
    }

    // Output in Excel format or fixed-width table format.

    if (opt.xl) {
        printf("%ld", tau);
        if (opt.tau0 != 1) printf(",%ld", tau * opt.tau0);
        printf(",%.6le,%ld", a.adev, n.adev);
        printf(",%.6le,%ld", a.mdev, n.mdev);
        printf(",%.6le,%ld", a.hdev, n.hdev);
        printf("\n");

    } else {
        printf("%*ld", opt.ld, tau);
        if (opt.tau0 != 1) {
            printf("%*ld", opt.ld + 4, tau * opt.tau0);
        }
        printf(" %.*le (%*ld)", opt.dp, a.adev, opt.ld, n.adev);
        printf(" %.*le (%*ld)", opt.dp, a.mdev, opt.ld, n.mdev);
        printf(" %.*le (%*ld)", opt.dp, a.hdev, opt.ld, n.hdev);
        printf("\n");
    }
}

int main (int argc, char *argv[])
{
    int i;
    long tau, decade;
    int get_options (int argc, char *argv[]);

    i = get_options(argc, argv);
    argc -= i;
    argv += i;
    opt.style = (argc > 1) ? argv[1] : "124";

    // Table: octaves of tau.

    if (strcmp(opt.style, "two") == 0 ||
        strncmp(opt.style, "octal", 3) == 0) {
        for (tau = 1; !stop; tau *= 2) {
            sigma_tau(tau);
        }
    } else

    // Table: decades of tau.

    if (strcmp(opt.style, "ten") == 0 ||
        strncmp(opt.style, "decimal", 3) == 0) {
        for (tau = 1; !stop; tau *= 10) {
            sigma_tau(tau);
        }
    } else

    // Table: all possible tau.

    if (strncmp(opt.style, "all", 3) == 0) {
        for (tau = 1; !stop; tau += 1) {
            sigma_tau(tau);
        }
    } else

    // Table: many tau (equal density of tau per decade); e.g., m50

    if (strncmp(opt.style, "m", 1) == 0) {
        long i, tau_prev, density;

        if (strcmp(opt.style, "many") == 0) {
            density = 72;
        } else {
            density = atoi(&opt.style[1]);
            if (density < 1) density = 1;
        }
        FPRINTF(stderr, "(many tau %d/decade)\n", density);
        tau_prev = 0;
        for (decade = 1; !stop; decade *= 10) {
            for (i = 0; i < density; i += 1) {
                tau = (long) (decade * pow(10.0, (double) i / density));
                if (tau > tau_prev) {
                    sigma_tau(tau);
                    tau_prev = tau;
                }
            }
        }
    } else

    // Table: nine tau (each digit) per decade.

    if (strncmp(opt.style, "nine", 3) == 0 ||
        strcmp(opt.style, "129") == 0) {
        for (decade = 1; !stop; decade *= 10) {
            for (tau = 1; tau < 10; tau += 1) {
                sigma_tau(tau * decade);
            }
        }
    } else

    // Table: 1-2-4 steps per decade.

    if (strcmp(opt.style, "124") == 0) {
        for (decade = 1; !stop; decade *= 10) {
            sigma_tau(1 * decade);
            sigma_tau(2 * decade);
            sigma_tau(4 * decade);
        }
    } else

    // Table: 1-2-5 steps per decade.

    if (strcmp(opt.style, "125") == 0) {
        for (decade = 1; !stop; decade *= 10) {
            sigma_tau(1 * decade);
            sigma_tau(2 * decade);
            sigma_tau(5 * decade);
        }
    } else

    // Range of tau.

    if (strstr(opt.style, ",") != NULL) {
        long i, j, k;
        i = j = k = 1;
        sscanf(opt.style, "%ld,%ld,%ld", &i, &j, &k);
        for (tau = i; tau <= j; tau += k) {
            sigma_tau(tau);
        }
    } else

    // One user-specified tau number.

    if (opt.style[0] >= '0' && opt.style[0] <= '9') {
        tau = atoi(opt.style);
        sigma_tau(tau);
    } else

    {
        fprintf(stderr, "Unknown table style or invalid tau: %s\n", opt.style);
        exit(1);
    }

    return 0;   // exit(0);
}

//
// Read file into memory.
// - One time interval reading per line, pre-scaled to seconds units.
// - Sample data (e.g., from a time interval counter on GPS receiver):
//   1.6862E-08
//   2.6083E-08
//   5.1941E-09
//   2.0985E-08
//   9.2359E-09
//   ...
//

double *read_data (FILE *f, long *count)
{
    double *data;
    long index, size;
    char line[1000], *p;
    char *get_field (char *line, int field);

    size = 1024 * sizeof data[0];
    data = malloc(size);
    for (index = 0; fgets(line, sizeof line, f) != NULL; index += 1) {
        if (index == size / (long)sizeof data[0]) {
            size *= 2;
            data = realloc(data, size);
        }
        if (data == NULL) {
            fprintf(stderr, "malloc failed (%ld bytes)\n", size);
            exit(1);
        }
        p = get_field(line, opt.field);
        if (p == NULL) {
            fprintf(stderr, "Line %ld: missing field %d: %s\n", index + 1, opt.field, line);
            exit(1);
        }
        data[index] = atof(p) * opt.scale;
    }
    *count = index;
    return data;
}

//
// Get field within line (first field is field 1).
// - Returns pointer to start of field within given input line.
// - Note: field is terminated with \0 (input line is modified).
//

char *get_field (char *line, int field)
{
    char *p, *start;
    int f;

    start = NULL;
    p = line;
    for (f = 0; *p && f < field; f += 1) {
        while (*p && isspace(*p)) {
            p++; // Skip leading white space
        }
        if (*p) {
            start = p;
        }
        while (*p && !isspace(*p)) {
            p++; // Skip over data field
        }
    }
    if (start && *p) {
        *p = '\0';
        return start;
    }
    return NULL;
}

int get_options (int argc, char *argv[])
{
    int args;

    // Set default options.

    opt.dp = 4;
    opt.ld = 6;
    opt.ovlp = 1;
    opt.field = 1;
    opt.scale = 1.0;
    opt.tau0 = 1;

    // Get user options.

    args = 0;
    while (argc > 1 && argv[1][0] == '/') {
        if (strcmp(argv[1], "/?") == 0 ||
            strcmp(argv[1], "/h") == 0 ||
            strcmp(argv[1], "/help") == 0) {
            fprintf(stderr, Usage);
            exit(0);
        } else

        if (strcmp(argv[1], "/q") == 0) {
            opt.quiet = 1;
        } else

        if (strcmp(argv[1], "/xl") == 0) {
            opt.xl = 1;
        } else

        if (strcmp(argv[1], "/b2b") == 0) {
            opt.ovlp = 0;
        } else

        if (strncmp(argv[1], "/f:", 3) == 0) {
            opt.field = atoi(strchr(argv[1], ':') + 1);
        } else

        if (strncmp(argv[1], "/s:", 3) == 0) {
            opt.scale = atof(strchr(argv[1], ':') + 1);
        } else

        if (strncmp(argv[1], "/dp:", 4) == 0) {
            opt.dp = atol(strchr(argv[1], ':') + 1);
        } else

        if (strncmp(argv[1], "/ld:", 4) == 0) {
            opt.ld = atol(strchr(argv[1], ':') + 1);
        } else

        if (strncmp(argv[1], "/t0:", 4) == 0 ||
            strncmp(argv[1], "/tau:", 5) == 0) {
            opt.tau0 = atol(strchr(argv[1], ':') + 1);
        } else

        {
            fprintf(stderr, "Error: unknown option: %s\n", argv[1]);
            exit(1);
        }
        argc -= 1;
        argv += 1;
        args += 1;
    }

    // Minimal sanity checks.

    if (opt.scale == 0.0) opt.scale = 1.0;
    if (opt.dp < 0) opt.dp = 0;
    if (opt.ld < 1) opt.ld = 1;
    if (opt.tau0 < 0) opt.tau0 = 1;

    return args;
}