#! /usr/bin/gawk -f
# Last edited on 2014-04-05 20:28:54 by stolfilocal

# Recomputes the {HUBI.S} and {WSL} columns of "00-DATA/2014-slumber-points-HUBI-BSTP.txt"

BEGIN \
  { abort = -1;
    
    # Weights for {Vh(02:30),Vh(03:30),Vh(04:30)}:
    wh0 = 0.053;
    wh1 = 0.690;
    wh2 = 0.071;
    # Soft threshold to map the {S} value to weight {WSL}:
    refS = 0.005; 
    # Next hour:
    nxh = 355; # For 2014-01-15 19:00:00
    
    NF_exp = 28; # Data lines must have this number of fields.
  } 
 
(abort >= 0) { exit abort; }

/[!]/ \
  { print; next; }

/^[ ]*([\#]|$)/ \
  { print; next; }
  
/^2014-[01][0-9]-[0-3][0-9] / \
  { # Remove trailing comments:
    gsub(/[ ]*[#].*$/, "", $0); 
    
    # Remove suffix "+" or "-" on unknown values:
    if ($0 ~ /[-+][ ]*([|]|$)/) { data_warning(("suffix '+'or '-' removed")); }
    gsub(/[-+][ ]*[|]/, " |", $0); 
    gsub(/[-+][ ]*$/, "", $0); 
    
    # Consistency checks:
    if (NF != NF_exp) { data_error(("wrong field count")); }
    if ($2 != "19:00:00") { data_error(("expected \"19:00:00\" on column 2 found \"" $(2) "\"")); }
    for (i = 3; i <= NF-1; i += 2) 
      { if ($(i) != "|") { data_error(("expected '|' on column " i " found \"" $(i) "\"")); } }
    
    # Grab relevant data fields:
    Vh0 = $(16);
    Vh1 = $(18);
    Vh2 = $(20)
    Vd = $(22);
    h = $(24);

    # More consistency:
    if (h+0 != nxh) { data_error(("expected hour " nxh " found \"" h "\"")); }
    nxh += 24;    
    
    # Compute {S}, {WSL}:
    S = vhmed(Vh0,Vh1,Vh2)/Vd/1000; 
    WSL = fmax(0.001,gauss(S/refS));
    
    # Insert in output
    $(26) = sprintf("%.4f", S); 
    $(28) = sprintf("%.3f", WSL); 
    print;
    next; 
  }

// \
  { printf "**bug [%s]\n", $0; 
    exit 1;
  }

function vhmed(Vh0,Vh1,Vh2,  Vh) 
  { Vh = (wh0*Vh0 + wh1*Vh1 + wh2*Vh2)/((Vh0 == 0 ? 0 : wh0) + (Vh1 == 0 ? 0 : wh1) + (Vh2 == 0 ? 0 : wh2));
    return Vh;
  }
  
function gauss(z)
  { if (z < 0) { z = -z; }
    if (z > 80) { return 0.0; }
    return exp(-0.5*z*z);
  }
         
function fmin(x,y)
  { return (x <= y ? x : y); }
         
function fmax(x,y)
  { return (x >= y ? x : y); }
           
function arg_error(msg)
  { printf "** %s\n", msg > "/dev/stderr"; 
    abort = 1;
    exit(abort);
  } 
 
function data_error(msg)
  { printf "%s:%s: ** %s\n", FILENAME, FNR, msg > "/dev/stderr"; 
    printf "  «%s»\n", $0 > "/dev/stderr"; 
    abort = 1;
    exit(abort);
  } 
         
function data_warning(msg)
  { printf "%s:%s: !! warning: %s\n", FILENAME, FNR, msg > "/dev/stderr"; 
    printf "  «%s»\n", $0 > "/dev/stderr"; 
  }