#!/usr/bin/env awk -f


function log10(n) {
    return log(n)/log(10.0)
}


function merge_meas(est_val, est_unc, meas_val, meas_unc) {

    # add sanity (staleness) check on previous estimate
    #z_est = (meas_val - est_val)/est_unc
    #z_meas = (meas_val - est_val)/meas_unc

    # Kalman filtering compares variances,
    # convert uncertainties to square units
    est_unc = (0.5*est_unc)^2.0
    meas_unc = (0.5*meas_unc)^2.0

    K = (est_unc)/((est_unc) + (meas_unc))
    est_val = est_val + K*(meas_val - est_val)
    est_unc = (est_unc*meas_unc)/(est_unc + meas_unc)

    # convert uncertainties back to original units
    est_unc = 2.0*sqrt(est_unc)
    meas_unc = 2.0*sqrt(meas_unc)

    return sprintf(OFMT OFS OFMT ORS, est_val, est_unc)

}


function sigfig(val, unc) {

    omag_val = log10(val)
    omag_unc = log10(unc)

    omag_high = sqrt((omag_val)^2.0 + (omag_unc)^2.0)
    omag_low = ((omag_val)^-1.0 + (omag_unc)^-1.0)^-1.0

    sfig_unc = sprintf("%.f", sqrt((omag_high - omag_low)^2.0 + 6.0^2.0))
    sfig_val = sprintf("%.f", sqrt(sfig_unc^2.0 + 9.0^2.0))

    ofmt_val = "%." sfig_val "g"
    ofmt_unc = "%." sfig_unc "g"

    return sprintf(ofmt_val OFS ofmt_unc, val, unc)

}


function unc_g(val) {

    # order of magnitude
    omag_val = log10(val)

    # trim punctuation from val
    sub("[eE].*$", "", val)
    gsub("[.+-]", "", val)
    n_digits = omag_val - (length(val) - 0.0)
    n_digits *= omag_val - (length(val) - 1.0)

    return sprintf("%.2g", 10.0^(n_digits))

}


function ck_boottime(unm) {

    if (unm ~/OpenBSD/) {
        "sysctl kern.boottime" | getline t0_k
        close("sysctl kern.boottime")
        sub("^.*=", "", t0_k)
        split(t0_k, t0_k_arr)
        t0_k_arr[2] = (index("JanFebMarAprMayJunJulAugSepOctNovDec", t0_k_arr[2]) + 2)/3.0
        gsub(":", " ", t0_k_arr[4])
        t0_k = sprintf("%04d %02d %02d %s", t0_k_arr[5], t0_k_arr[2], t0_k_arr[3], t0_k_arr[4])
        t0_k = mktime(t0_k)
    }

    if (unm ~/FreeBSD/ || unm ~/Darwin/) {
        "sysctl kern.boottime" | getline t0_k
        close("sysctl kern.boottime")
        sub("^.*{", "", t0_k)
        sub("}.*$", "", t0_k)
        split(t0_k, t0_k_arr, ",")
        for (n in t0_k_arr) {
            sub("^.*= ", "", t0_k_arr[n])
        }
        t0_k = t0_k_arr[1] "." t0_k_arr[2]

    }

    t0_k_unc = ck_time(unm)

    return sprintf("%s %s", t0_k, t0_k_unc)

}


function ck_time(unm) {

    if (unm ~ /FreeBSD/ || unm ~ /Linux/) {
        "date +%s.%N" | getline t
        close("date +%s.%N")
    }
    #if (unm ~ /Darwin/) {
    else {
        "date +%s" | getline t
        close("date +%s")
    }
    #else {
    #    t = systime()
    #}

    return t

}


function ck_uptime(unm) {

    # check uptime cmd
    "uptime" | getline t_up_cmd
    close("uptime")

    # evaluate measured uptime
    t_up_meas[1] = 0.0
    t_up_meas[2] = ck_time(unm)

    sub("^.*up ", "", t_up_cmd)
    sub(", load.*$", "", t_up_cmd)
    split(t_up_cmd, t_up_cmd_arr, ",")

    for (i in t_up_cmd_arr) {

        if (t_up_cmd_arr[i] ~ /day/) {
            split(t_up_cmd_arr[i], days)
            t_up_meas[1] += 86400.0*(days[1] + 0.0)
            (t_up_meas[2] > 86400.0) ? t_up_meas[2] = 86400.0 : t_up_meas[2] += 0.0
        }

        if (t_up_cmd_arr[i] ~ /hr/) {
            split(t_up_cmd_arr[i], hrs)
            t_up_meas[1] += 3600.0*(hrs[1] + 0.0)
            (t_up_meas[2] > 3600.0) ? t_up_meas[2] = 3600.0 : t_up_meas[2] += 0.0
        }

        if (t_up_cmd_arr[i] ~ /min/) {
            split(t_up_cmd_arr[i], mins)
            t_up_meas[1] += (60.0*(mins[1] + 0.0))
            (t_up_meas[2] > 60.0) ? t_up_meas[2] = 60.0 : t_up_meas[2] += 0.0
        }

        #if (t_up_cmd_arr[i] ~ /sec/) {
        #    split(t_up_cmd_arr[i], sec)
        #    t_up_meas[1] += sec[1]
        #    (t_up_meas[2] > 1.0) ? t_up_meas[2] = 1.0 : t_up_meas[2] += 0.0
        #}

        if (t_up_cmd_arr[i] ~ /:/) {
            split(t_up_cmd_arr[i], hrs_min, ":")
            t_up_meas[1] += 3600.0*(hrs_min[1] + 0.0)
            (t_up_meas[2] > 3600.0) ? t_up_meas[2] = 3600.0 : t_up_meas[2] += 0.0
            t_up_meas[1] += 60.0*(hrs_min[2] + 0.0)
            (t_up_meas[2] > 60.0) ? t_up_meas[2] = 60.0 : t_up_meas[2] += 0.0
        }

    }

    if (uname ~ /FreeBSD/ || uname ~ /Linux/) {
        split(ck_proc_uptime(uname), t_pu_meas)
        split(merge_meas(t_up_meas[1], t_up_meas[2], t_pu_meas[1], 0.1), t_up_meas)
    }

    return sigfig(t_up_meas[1], t_up_meas[2])

}


function ck_proc_uptime(unm) {

    proc_file = ""
    if (unm ~ /FreeBSD/) {
        proc_file = "/compat/linux/proc/uptime"
    }
    if (unm ~ /Linux/) {
        proc_file = "/proc/uptime"
    }

    getline proc_uptime < proc_file
    close(proc_file)

    return proc_uptime

}


BEGIN {

    OFMT="%f"
    pi = 4.0*atan2(1.0, 1.0)
    c0 = 299792458  # m/sec

    "uname" | getline uname
    close("uname")

    ## check ARGV for previous estimates
    if (ARGC) {
        ARGV[1] ? t_prev[1] = ARGV[1] : t_prev[1] = 0.0
        ARGV[2] ? t_prev[2] = ARGV[2] : t_prev[2] = ck_time(uname)
        ARGV[3] ? t_up_prev[1] = ARGV[3] : t_up_prev[1] = 0.0
        ARGV[4] ? t_up_prev[2] = ARGV[4] : t_up_prev[2] = ck_time(uname)
        ARGV[5] ? t0_prev[1] = ARGV[5] : split(ck_boottime(uname), t0_prev)
        ARGV[6] ? t0_prev[2] = ARGV[6] : split(ck_boottime(uname), t0_prev)
    }

    ## check systime, check uptime, check systime again
    t_meas[1] = ck_time(uname)
    split(ck_uptime(uname), t_up_meas)
    t_meas[2] = ck_time(uname)

    ## convert to [val, unc] vector
    t_meas[2] -= t_meas[1]
    t_meas[1] += 0.5*t_meas[2]
    if ( t_meas[2] == 0.0 ) {
        t_meas[2] = unc_g(t_meas[1])
        #t_meas[2] = sqrt(t_meas[2]^2.0 + 1.0^2.0)
    }
    t_meas_str = sigfig(t_meas[1], t_meas[2])

    ## estimate dt
    dt[1] = t_meas[1] - t_prev[1]
    dt[2] = t_meas[2] + t_prev[2]

    ## predict uptime based on dt, t_up(new) = t_up(old) + dt
    t_up_est[1] = t_up_prev[1] + dt[1]
    t_up_est[2] = t_up_prev[2] + dt[2]

    ## check for difference between dt and dt_up
    t_up_est[2] += sqrt((t_up_est[1] - t_up_prev[1] - dt[1])^2.0)
    t_up_prev[2] += sqrt((t_up_est[1] - t_up_prev[1] - dt[1])^2.0)

    ## alt. t_up prediction based on t0_prev
    t_up_alt[1] = t_meas[1] - t0_prev[1]
    t_up_alt[2] = t_meas[2] + t0_prev[2]

    ## dt_up_alt
    dt_up_alt[1] = t_up_alt[1] - t_up_prev[1]
    dt_up_alt[2] = t_up_alt[2] + t_up_prev[2]

    ## additional error terms
    t_up_alt[2] += sqrt((dt_up_alt[1] - dt[1])^2.0)
    t_up_alt[2] += sqrt((t_up_alt[1] - t_up_est[1])^2.0)
    t_up_est[2] += sqrt((t_up_est[1] - t_up_alt[1])^2.0)

    ## merge alt. prediction
    split(merge_meas(t_up_est[1], t_up_est[2], t_up_alt[1], t_up_alt[2]), t_up_est)

    ## predict boot time, t0
    t0_est[1] = t_meas[1] - t_up_est[1]
    t0_est[2] = t_meas[2] + t_up_est[2]

    ## merge predicted and measured uptimes
    split(merge_meas(t_up_est[1], t_up_est[2], t_up_meas[1], t_up_meas[2]), t_up_est)
    t_up_est_str = sigfig(t_up_est[1], t_up_est[2])

    ## evaluate measured boot time, t0
    t0_meas[1] = t_meas[1] - t_up_est[1]
    t0_meas[2] = t_meas[2] + t_up_est[2]

    ## additional error terms
    t0_meas[2] += sqrt((t0_meas[1] - t0_est[1])^2.0)
    t0_est[2] += sqrt((t0_meas[1] - t0_est[1])^2.0)

    ## merge previous and updated boot time, t0
    split(merge_meas(t0_meas[1], t0_meas[2], t0_est[1], t0_est[2]), t0_est)

    t0_est_str = sigfig(t0_est[1], t0_est[2])

    print(t_meas_str, t_up_est_str, t0_est_str)

}