//
// gh.exe
//
// Estimate acceleration of gravity for given latitude
// and elevation/altitude. Nominal value is 9.8 m/s^2.
//
// Then compute the relativistic effect on clock rate and
// display as both frequency offset (clocks run fast at
// higher altitudes) and time drift (nanoseconds gained
// per day).
//
// While we're at it calculate period of pendulum swing.
//
// 19-Sep-2005 /tvb
//

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

#define PI (4 * atan(1.0))

#define C 299792458.0
#define C2 (C * C)

#define SIN(a) ( sin(2.0 * PI * (a) / 360.0) )
#define SIN2(a) ( SIN(a) * SIN(a) )

#define G1 9.7803184
#define B1 0.0053024
#define B2 (-0.0000059)
#define H1 1.1e-6
#define H2 3.086e-6

#define METERS_TO_FEET(m) ((m) / 0.0254 / 12.0)
#define METERS_TO_INCHES(m) ((m) / 0.0254)

void main (int argc, char *argv[])
{
    double lat, ele, alt, g, h, df, T, L1, L2;

    if (argc == 1) {
        fprintf(stderr, "Usage: %s [latitude] [elevation] [altitude]\n", argv[0]);
        exit(1);
    }

    // Get degrees latitude.
    // Get meters elevation of ground above MSL.
    // Get meters altitude in air above ground.

    lat = (argc > 1) ? atof(argv[1]) : 45.0;
    ele = (argc > 2) ? atof(argv[2]) : 0;
    alt = (argc > 3) ? atof(argv[3]) : 0;

    // Compute estimate of g.

    g = G1 * (1 + B1 * SIN2(lat) + B2 * SIN2(lat * 2.0));

    // Add to g at higher elevations due to extra mass.
    // Subtract from g when above the earth due to extra distance.
    // Not sure if this is the right way to handle mountains or airplanes.

    g += (H1 * ele);
    g -= (H2 * alt);

    // Compute general relativistic clock effect.

    h = ele + alt;
    df = g * h / C2;

    // Prevent equation misuse.

    if (h > 10000) {
        fprintf(stderr, "I don't trust these formulas at space altitudes.\n");
        (argc > 4) ? 0 : exit(0);
    }

    // Calculate 1 and 2 seconds pendulum swing.

    T = 2 * PI * sqrt (1.0 / g);
    L1 = g * pow(1.0 / (2 * PI), 2.0);
    L2 = g * pow(2.0 / (2 * PI), 2.0);

    // Display results.

    printf("%12.2lf degrees latitude\n", lat);
    printf("%12lg meters high (%lg ele + %lg alt) %.0lf feet\n",
           h, ele, alt, METERS_TO_FEET(h));
    printf("%12.6lf g\n", g);
    printf("%12.6lf s period of meter pendulum\n", T);
    printf("%12.6lf m length of second swing pendulum (%.3lf inches)\n",
           L1, METERS_TO_INCHES(L1));
    printf("%12.6lf m length of second period pendulum (%.3lf inches)\n",
           L2, METERS_TO_INCHES(L2));
    printf("%12.4le GR clock effect per meter\n", g / C2);
    printf("%12.4le net frequency offset\n", df);
    printf("%12.3lf ps/hour time dilation\n", df * 3600.0 * 1e12);
    printf("%12.3lf ns/day time dilation\n", df * 86400.0 * 1e9);
}