#! /usr/bin/gawk -f
# Last edited on 2017-05-02 03:36:42 by jstolfi

# Convert Y coordinates of Conde eng vapor pressure plot to Kelvin temps.

BEGIN {
  
  if (inField == "") { printf "** inField?\n" > "/dev/stderr"; exit(1); }
  inField += 0;

  if (outField == "") { printf "** outField?\n" > "/dev/stderr"; exit(1); }
  outField += 0;

  split("", Ytb);
  split("", Ttb);
  ntb = 44;
  
  Ytb[ 0] = 686.0; Ttb[ 0] = 285.0;
  Ytb[ 1] = 667.0; Ttb[ 1] = 290.0;
  Ytb[ 2] = 648.0; Ttb[ 2] = 295.0;
  Ytb[ 3] = 630.0; Ttb[ 3] = 299.8;
  Ytb[ 4] = 611.0; Ttb[ 4] = 305.1;
  Ytb[ 5] = 594.0; Ttb[ 5] = 310.0;
  Ytb[ 6] = 577.0; Ttb[ 6] = 315.0;
  Ytb[ 7] = 562.0; Ttb[ 7] = 319.5;
  Ytb[ 8] = 544.0; Ttb[ 8] = 325.1;
  Ytb[ 9] = 528.0; Ttb[ 9] = 330.3;
  Ytb[10] = 514.0; Ttb[10] = 335.0;
  Ytb[11] = 500.0; Ttb[11] = 339.9;
  Ytb[12] = 486.0; Ttb[12] = 344.9;
  Ytb[13] = 472.0; Ttb[13] = 350.1;
  Ytb[14] = 459.0; Ttb[14] = 355.0;
  Ytb[15] = 446.0; Ttb[15] = 360.0;
  Ytb[16] = 433.0; Ttb[16] = 365.3;
  Ytb[17] = 422.0; Ttb[17] = 369.9;
  Ytb[18] = 410.0; Ttb[18] = 374.9;
  Ytb[19] = 398.0; Ttb[19] = 379.9;
  Ytb[20] = 386.0; Ttb[20] = 385.1;
  Ytb[21] = 375.0; Ttb[21] = 390.0;
  Ytb[22] = 365.0; Ttb[22] = 394.7;
  Ytb[23] = 354.0; Ttb[23] = 400.0;
  Ytb[24] = 344.0; Ttb[24] = 405.0;
  Ytb[25] = 334.0; Ttb[25] = 410.0;
  Ytb[26] = 324.0; Ttb[26] = 415.2;
  Ytb[27] = 315.0; Ttb[27] = 420.0;
  Ytb[28] = 306.0; Ttb[28] = 425.0;
  Ytb[29] = 298.0; Ttb[29] = 429.5;
  Ytb[30] = 288.0; Ttb[30] = 435.2;
  Ytb[31] = 280.0; Ttb[31] = 439.9;
  Ytb[32] = 272.0; Ttb[32] = 444.7;
  Ytb[33] = 262.0; Ttb[33] = 450.8;
  Ytb[34] = 247.0; Ttb[34] = 460.4;
  Ytb[35] = 231.0; Ttb[35] = 470.8;
  Ytb[36] = 217.0; Ttb[36] = 480.3;
  Ytb[37] = 204.0; Ttb[37] = 489.5;
  Ytb[38] = 190.0; Ttb[38] = 499.9;
  Ytb[39] = 176.0; Ttb[39] = 510.6;
  Ytb[40] = 164.0; Ttb[40] = 520.1;
  Ytb[41] = 152.0; Ttb[41] = 529.9;
  Ytb[42] = 140.0; Ttb[42] = 540.2;
  Ytb[43] = 129.0; Ttb[43] = 550.0;
  
  if (dump)
    { 
      # Checking:
      for (k = 1; k <= ntb-2; k++)
        { Tver = convert_conde_temp_interpolate(Ytb[k],Ytb[k-1],Ttb[k-1],Ytb[k+1],Ttb[k+1]);
          printf "%5.1f %5.1f %5.1f %6.2f\n", Ytb[k], Ttb[k], Tver, Ttb[k]-Tver > "/dev/stderr"; 
        }
    }
}

/^ *[0-9]/ { 
  Y = $(inField)+0;
  T = convert_conde_temp(Y); 
  $(outField) = sprintf(" %5.1f ", T - 273.15);
  print;
  next;
}

// { print; next; }

function convert_conde_temp(Y, k0,k1,T)
  { 
    k0 = convert_conde_temp_lookup(Y);
    k1 = k0 + 1;
    T = convert_conde_temp_interpolate(Y,Ytb[k0],Ttb[k0],Ytb[k1],Ttb[k1]);
    return T;
  }

function convert_conde_temp_lookup(Y,  k)
  { 
    # Find {k0,k1} so that {Ytb[k0],Ytb[k1]} are the values to interpolate:
    if (Y >= Ytb[0])
      { return 0; }
    else if (Y <= Ytb[ntb-1])
      { return ntb-2; }
    else
      { for (k = 1; k < ntb; k++) 
          { if (Y >= Ytb[k]) { return k-1; } }
        printf "** bug lookup\n" > "/dev/stderr"; exit(1);
      }
  } 
             
function convert_conde_temp_interpolate(Y,Y0,T0,Y1,T1, s)
  { 
    if (Y0 <= Y1) { printf "** bug Y0 = %5.1f  Y1 = %5.1f\n", Y0, Y1 > "/dev/stderr"; exit(1); }
    if (T0 > T1) { printf "** bug T0 = %5.1f T1 = %5.1f\n", T0, T1 > "/dev/stderr"; exit(1); }
    if ((Y > Y0) || (Y < Y1))
      { printf "!! extrapolating Y = %5.1f Y0 = %5.1f  Y1 = %5.1f\n", Y,Y0,Y1 > "/dev/stderr"; }
    s = (Y - Y0)/(Y1 - Y0);
    return (1-s)*T0 + s*T1;
  }