//
// Estimate UT1 - UTC based on Earth orientation parameters.
//
// See:
//  http://maia.usno.navy.mil/
//  ftp://maia.usno.navy.mil/ser7/ser7.dat
//
// I have no idea how valid this calculation is. It seems to be
// within a second for the past 3 years but off by a factor of
// five for the last 30 years (it back-predicts 6 leap seconds
// while the actual number was more than 30).
//
// 2002-02-14 tvb
//

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

double
calc (double MJD)
{
    double T;
    double pi = 3.14159265358979323846;
    double A, C, x, y;
    double UT2_UT1;
    double UT1_UTC;

    /* T is the date in Besselian years */
    T = 2000.000 + (MJD - 51544.03) / 365.2422;

    UT2_UT1 =   0.022 * sin(2 * pi * T)
              - 0.012 * cos(2 * pi * T)
              - 0.006 * sin(4 * pi * T)
              + 0.007 * cos(4 * pi * T);

    A = 2 * pi * (MJD - 52317) / 365.25;
    C = 2 * pi * (MJD - 52317) / 435.0;

    x = 0.0590 - 0.1242 * cos(A)
               - 0.0024 * sin(A)
               - 0.0972 * cos(C)
               + 0.1089 * sin(C);

    y = 0.3280 - 0.0025 * cos(A)
               + 0.1101 * sin(A)
               + 0.1089 * cos(C)
               + 0.0972 * sin(C);

    UT1_UTC = -0.1557 - 0.00055 * (MJD - 52322) - UT2_UT1;

    return UT1_UTC;
}

#define MJD_1970_01_01 40587

void
main (int argc, char *argv[])
{
    double mjd;
    double end;
    double inc;

    if (argc > 1) {
        sscanf(argv[1], "%lf", &mjd);
    } else {
        time_t t;
        time(&t);
        mjd = MJD_1970_01_01 + t / 86400;
    }

    if (argc > 2) {
        sscanf(argv[2], "%lf", &end);
    } else {
        end = mjd;
    }

    if (argc > 3) {
        sscanf(argv[3], "%lf", &inc);
    } else {
        inc = 1;
    }

    while (mjd <= end) {
        printf("MJD = %.0lf  UT1-UTC (estimate) = %.2lf s\n", mjd, calc(mjd));
        mjd += inc;
    }
}

#if 0

http://maia.usno.navy.mil/

Modified Julian Date (MJD):
MJD = Julian Date - 2,400,000.5 days
UT2-UT1 = 0.022 Sin(2*pi*T) - 0.012 cos(2*pi*T) - 0.006 sin(4*pi*T) + 0.007 cos(4*pi*T)
where pi = 3.14159265... and T is the date in Besselian years.
T, the Besselian date, is given by T = 2000.000 + (MJD - 51544.03) / 365.2422
TT = TAI + 32.184 seconds
since 01 January 1999, the cumulative leap second count is:
TAI-UTC(BIPM) = 32 seconds (exactly)
and GPS Time-UTC(BIPM) = 13 seconds

#endif

#if 0

ftp://maia.usno.navy.mil/ser7/ser7.dat

GENERAL INFORMATION:
         To receive this information electronically, contact ser7@maia.
            usno.navy.mil or use http://maia.usno.navy.mil/.
         MJD = Julian Date - 2 400 000.5 days
         UT2-UT1 = 0.022 sin(2*pi*T) - 0.012 cos(2*pi*T)
                                     - 0.006 sin(4*pi*T) + 0.007 cos(4*pi*T)
            where pi = 3.14159265... and T is the date in Besselian years.
         TT = TAI + 32.184 seconds
         DUT1= (UT1-UTC) transmitted with time signals
             =  -0.1 seconds beginning 04 October 2001 at 0000 UTC
         Beginning 1 January 1999:
            TAI-UTC(BIPM) = 32.000 000 seconds

PREDICTIONS:
      The following formulas will not reproduce the predictions given below,
      but may be used to extend the predictions beyond the end of this table.

      x =   .0590 -  .1242 cos A -  .0024 sin A -  .0972 cos C +  .1089 sin C
      y =   .3280 -  .0025 cos A +  .1101 sin A +  .1089 cos C +  .0972 sin C
         UT1-UTC =  -.1557 -  .00055 (MJD - 52322) - (UT2-UT1)

      where A = 2*pi*(MJD-52317)/365.25 and C = 2*pi*(MJD-52317)/435.

         TAI-UTC(MJD 52318) = 32.0
      The accuracy may be estimated from the expressions:
      S x,y = 0.0042 (MJD-52317)**0.28    S t = 0.0003 (MJD-52322)**0.75
      Estimated accuracies are:  Predictions     10 d   20 d   30 d   40 d
                                 Polar coord's  0.004  0.006  0.009  0.011
                                 UT1-UTC        0.0017 0.0028 0.0039 0.0048

#endif