//
// Calculate 3-cornered hat directly from TSC 5110A OCR files.
//
// 27-Mar-2003 tvb
//

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

#if 0

Sample input file:

  0.1 2.49+/-.078x10-11
  0.2 1.81+/-.057x10-11
  0.4 1.83+/-.057x10-11

#endif

#define ROWMAX 25

double RowToTau (int Row)
{
    double m, e;
    int i;

    i = Row / 3;
    e = i - 2;
    switch (Row % 3) {
        case 0 : m = 1.0; break;
        case 1 : m = 2.0; break;
        case 2 : m = 4.0; break;
    }

    return m * pow(10.0, e);
}

int TauToRow (double Tau)
{
    char buf[100];
    int Row, i, n;

    sprintf(buf, "%.0le", Tau);
    switch (buf[0]) {
    case '1' :
        i = 0;
        break;

    case '2' :
        i = 1;
        break;

    case '4' :
        i = 2;
        break;

    default :
        fprintf(stderr, "Unexpected tau %s %g\n", buf, Tau);
        exit(1);
    }
    n = atoi(&buf[2]);

    //
    // Calculate Row so that
    // tau .01 is row 0, .02 is row 1, .04 is row 2,
    // tau .1 is row 3, .2 is row 4, .4 is row 5,
    // tau 1 is row 6, etc.
    // 8 decades plus 1 = 25 (.01 .1 1 10 100 1k 10k 100k 1M)

    Row = 6 + (3 * n) + i;
    if (Row < 0 || Row >= ROWMAX) {
        fprintf(stderr, "Unexpected tau %g (Row %d)\n", Tau, Row);
        exit(1);
    }

    return Row;
}

void ReadSigma (char *Path, double *Data)
{
    FILE *File;
    char Line[100];

    File = fopen(Path, "r");
    if (File == NULL) {
        fprintf(stderr, "Cannot open: %s\n", Path);
        exit(1);
    }

    while (fgets(Line, sizeof Line, File) != NULL) {
        int n;
        int Row;
        double Tau, Sigma, Error;
        double Mantissa, Uncertainty, Exponent;

        if (strstr(Line, "+/-") != 0) {
            n = sscanf(Line, "%lf%lf+/-%lfx10%lf",
                &Tau, &Mantissa, &Uncertainty, &Exponent);
            if (n != 4) {
                fprintf(stderr, "Bad format on line: %s", Line);
                // exit(1);

            } else {
                Row = TauToRow(Tau);
                Sigma = Mantissa * pow(10.0, Exponent);
                Error = Uncertainty * pow(10.0, Exponent);
                Data[Row] = Sigma;
            }
        }
    }

    return;
}

void main (int argc, char *argv[])
{
    double A[ROWMAX], B[ROWMAX], C[ROWMAX];
    double AB[ROWMAX], BC[ROWMAX], AC[ROWMAX];
    int i;

    if (argc < 4) {
        fprintf(stderr, "Usage: AB BC AC\n");
        exit(1);
    }

    for (i = 0; i < ROWMAX; i += 1) {
        AB[i] = BC[i] = AC[i] = 0.0;
    }

    ReadSigma(argv[1], AB);
    ReadSigma(argv[2], BC);
    ReadSigma(argv[3], AC);

    printf("%6s %10s %10s %10s %10s %10s %10s\n",
           "tau", "AB", "BC", "AC", "A", "B", "C");

#define SQR(x) (x * x)
#define SQRT(x) ( x <= 0.0 ? sqrt(-x) : sqrt(x) )

    for (i = 0; i < ROWMAX; i += 1) {

        A[i] = SQRT( (0 + SQR(AB[i]) - SQR(BC[i]) + SQR(AC[i])) / 2.0 );
        B[i] = SQRT( (0 + SQR(AB[i]) + SQR(BC[i]) - SQR(AC[i])) / 2.0 );
        C[i] = SQRT( (0 - SQR(AB[i]) + SQR(BC[i]) + SQR(AC[i])) / 2.0 );

        printf("%.3le %.3le %.3le %.3le %.3le %.3le %.3le\n",
               RowToTau(i), AB[i], BC[i], AC[i], A[i], B[i], C[i]);
    }

    return;
}