/** @file  psTime.c
 *
 *  @brief Definitions for time, time utilities, and conversion functions for use with psLib astronomy
 *  functions.
 *
 *  A collection of functions are required by psLib to manipulate time data. These functions primarily consist
 *  of conversions between specific time formats.  They use the UNIX timeval time system as the
 *  base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated.
 *
 *  @author Ross Harman, MHPCC
 *
 *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-12-17 20:47:08 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

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

#include "psTime.h"
#include "psError.h"
#include "psMemory.h"
#include "psAbort.h"
#include "psImage.h"
#include "psCoord.h"
#include "psString.h"
#include "psMetadata.h"
#include "psMetadataIO.h"
#include "psLookupTable.h"
#include "psConstants.h"
#include "psAstronomyErrors.h"

#ifndef TIME_CONFIG_FILE
#define TIME_CONFIG_FILE "../../config/psTime.config"
#pragma warning TIME_CONFIG_FILE was not defined in the makefile.
#endif

#define MAX_STRING_LENGTH 256

/** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */
#define S2R (7.272205216643039903848711535369e-5)

/** Two times pi with double precision accuracy */
#define TWOPI (2.0*M_PI)

/** Conversion from radians to degrees */
#define R2DEG = (180.0/M_PI)

                /** Maximum length of time string */
                #define MAX_TIME_STRING_LENGTH 256

                /** Seconds per minute */
                #define  SEC_PER_MINUTE 60.0

                /** Seconds per hour */
                #define  SEC_PER_HOUR (60.0*SEC_PER_MINUTE)

                /** Seconds per day */
                #define  SEC_PER_DAY (24.0*SEC_PER_HOUR)

                /** Seconds per year */
                #define  SEC_PER_YEAR (365.0*SEC_PER_DAY)

                /** Microseconds per day */
                #define USEC_PER_DAY 86400000000.0

                /** Time metadata read from config file */
                static psMetadata *timeMetadata = NULL;

/** Static function prototypes */
static char *cleanString(char *inString, int sLen);
static char* getToken(char **inString, char *delimiter, psParseErrorType *status);
psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables);

/** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null
 *  terminated copy of the original input string. */
static char *cleanString(char *inString, int sLen)
{
    char *ptrB = NULL;
    char *ptrE = NULL;
    char *cleaned = NULL;


    ptrB = inString;

    /* Skip over leading # or whitespace */
    while (isspace(*ptrB) || *ptrB=='#') {
        ptrB++;
    }

    /* Skip over trailing whitespace, null terminators, and # characters */
    ptrE = inString + sLen;
    while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') {
        ptrE--;
    }

    // Length, sLen, does not include '\0'
    sLen = ptrE - ptrB + 1;

    // Adds '\0' to end of string and +1 to sLen
    cleaned = psStringNCopy(ptrB, sLen);

    return cleaned;
}

/** Returns cleaned token based on delimiter, but not including delimiter. Also changes the pointer location
 * the beginning of the string. Tokens are newly allocated null terminated strings. */
static char* getToken(char **inString, char *delimiter, psParseErrorType *status)
{
    char *cleanToken = NULL;
    int sLen = 0;


    // Skip over leading whitespace
    while(isspace(**inString)) {
        (*inString)++;
    }

    // Length of token, not including delimiter
    sLen = strcspn(*inString, delimiter);
    if(sLen) {

        // Create new, cleaned, and null terminated token
        cleanToken = cleanString(*inString, sLen);

        // Move to end of token
        (*inString) += sLen;
    } else if(**inString!='\0' && sLen==0) {
        *status = PS_PARSE_ERROR_GENERAL;
    }

    return cleanToken;
}

/** Searches time tables in priority order and performs interpolation if input index value is within a table.
 * If the index value is out of range, the status is set accordingly. */
psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables)
{
    psU32 i = 0;
    char *tableName = NULL;
    psF64 result = 0.0;
    psLookupTable *table = NULL;
    psMetadataItem *tableMetadataItem = NULL;



    // Check time metadata. Function call reports errors.
    if(timeMetadata == NULL) {
        if(!p_psTimeInit(TIME_CONFIG_FILE))
            return 0.0;
    }

    // Search each table in priority order: ser7, eopc,finals
    for(i=0; i<nTables; i++) {

        // Get table from metadata
        tableName = metadataTableNames[i];
        tableMetadataItem = psMetadataLookup(timeMetadata, tableName);
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                    tableName);
            return 0.0;
        }
        table = (psLookupTable*)tableMetadataItem->data.V;
        PS_PTR_CHECK_NULL(table,0.0);

        // Attempt to interpolate table
        *status = PS_LOOKUP_SUCCESS;
        result = psLookupTableInterpolate(table, index, 2, status);
        if(*status == PS_LOOKUP_SUCCESS) {
            return result;
        } else if(*status == PS_LOOKUP_ERROR) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED_NAME,
                    tableName);
            return 0.0;
        }
    }

    return result;
}

bool p_psTimeInit(const char *fileName)
{
    bool foundTable = false;
    char *tableDir = NULL;
    char *tableNames = NULL;
    char *namesPtr = NULL;
    char *metadataNamesPtr = NULL;
    char *tableName = NULL;
    char *fullTableName = NULL;
    psS32 i = 0;
    psS32 j = 0;
    psS32 numTables = 0;
    psU32 nFail = 0;
    psVector *tablesFrom = NULL;
    psVector *tablesTo = NULL;
    psMetadataItem *metadataItem = NULL;
    psLookupTable *table = NULL;
    psParseErrorType status = PS_PARSE_SUCCESS;
    char metadataTableNames[4][MAX_STRING_LENGTH] = {"ser7", "eopc",  "finals", "tai"};


    // Read time config file
    timeMetadata = psMetadataParseConfig(timeMetadata, &nFail, fileName, true);
    if(timeMetadata == NULL) {
        return false;
    } else if(nFail != 0) {
        return false;
    }

    // Get number of tables
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.n");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.n");
        return false;
    }
    numTables = (psS32)metadataItem->data.S32;

    // Get lower range of tables
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.from");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.from");
        return false;
    }
    tablesFrom = psVectorCopy(tablesFrom, metadataItem->data.V, PS_TYPE_F64);
    if(tablesFrom->n != numTables) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesFrom->n, numTables);
        psFree(tablesFrom);
        return false;
    }

    // Get upper range of tables
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.to");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.to");
        return false;
    }
    tablesTo = psVectorCopy(tablesTo, metadataItem->data.V, PS_TYPE_F64);
    if(tablesTo->n != numTables) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesTo->n, numTables);
        psFree(tablesFrom);
        psFree(tablesTo);
        return false;
    }

    // Get path to time data files
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.dir");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.dir");
        psFree(tablesFrom);
        psFree(tablesTo);
        return false;
    }
    tableDir = psStringCopy(metadataItem->data.V);

    // Table file names
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.files");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.files");

        psFree(tablesFrom);
        psFree(tablesTo);
        psFree(tableDir);
        return false;
    }
    tableNames = psStringCopy(metadataItem->data.V);

    // Read time tables
    namesPtr = tableNames;
    while((tableName=getToken(&namesPtr, " ", &status)) != NULL) {

        // Form path with table name, adding one to length for last '/' that may not occur in string in cong file
        fullTableName = (char*)psAlloc(strlen(tableDir)+strlen(tableName)+1+1);

        // Old strings may come back from psAlloc(), so set initial position to EOL
        fullTableName[0]='\0';
        strcat(fullTableName, tableDir);
        strcat(fullTableName, "/");
        strcat(fullTableName, tableName);

        // Create and read table
        if(i < numTables) {
            table = psLookupTableAlloc(fullTableName, tablesFrom->data.F64[i], tablesTo->data.F64[i]);
            table = psLookupTableRead(table);
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables);
        }

        // Place tables into metadata slightly altered names as keys to create consistent naming conventions
        foundTable = false;
        for(j=0; j<numTables; j++) {
            metadataNamesPtr = strstr(tableName, metadataTableNames[j]);
            if(metadataNamesPtr != NULL) {
                psMetadataAdd(timeMetadata, PS_LIST_TAIL, strcat(metadataTableNames[j], "Table"),
                              PS_TYPE_PTR, PS_META_LOOKUPTABLE, NULL, table);
                foundTable = true;
            } else if(foundTable==false && j==numTables-1) {
                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, j, numTables);
            }
        }

        psFree(fullTableName);
        psFree(tableName);
        i++;
    }

    if(numTables != i) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i, numTables);
    }

    psFree(tableDir);
    psFree(tableNames);
    psFree(tablesFrom);
    psFree(tablesTo);

    return true;
}

bool p_psTimeFinalize(void)
{
    if(timeMetadata != NULL) {
        psFree(timeMetadata);
        timeMetadata = NULL;
    }

    return true;
}

psTime* psTimeAlloc(psTimeType type)
{
    psTime *outTime = NULL;


    // Error checks
    if(type!=PS_TIME_TAI && type!=PS_TIME_UTC) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psTime_TYPE_UNKNOWN,
                type);
        return NULL;
    }

    outTime = (psTime*)psAlloc(sizeof(psTime));

    outTime->sec = 0;
    outTime->usec = 0;
    outTime->type = type;

    return outTime;
}

psTime* psTimeGetTime(psTimeType type)
{
    struct timeval now;
    psTime *time = NULL;


    // Allocate psTime struct
    time = psTimeAlloc(type);

    if (gettimeofday(&now, (struct timezone *)0) == -1) {
        psError(PS_ERR_OS_CALL_FAILED, true,
                PS_ERRORTEXT_psTime_GET_TOD_FAILED);
        return NULL;
    }

    // Convert timeval time to psTime
    time->sec = now.tv_sec;
    time->usec = now.tv_usec;

    // Add most leapseconds to UTC time to get TAI time if necessary
    if(type == PS_TIME_TAI) {
        time->sec += psTimeGetTAIDelta(time);
    }

    return time;
}

psTime* psTimeConvert(psTime *time, psTimeType type)
{
    double delta = 0.0;


    // Error checks
    PS_PTR_CHECK_NULL(time,NULL);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);

    if (time->type == type) { // time already right type.  That was easy!
        return time;
    }

    delta = psTimeGetTAIDelta(time);

    if (type == PS_TIME_UTC) {
        time->sec = time->sec - delta;    // User wants UTC time. Subtract leapseconds.
    } else if(type == PS_TIME_TAI) {
        time->sec = time->sec + delta;    // User wants TAI time. Add leapseconds.
    } else {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type);
    }

    return time;
}

double psTimeToLST(psTime *time, double longitude)
{
    psF64  jdTdtDays    =  0.0;
    psF64  jdUt1Days    =  0.0;
    psF64  mjdUt1Days   =  0.0;
    psF64  lmstRad      =  0.0;
    psF64  fracDays     =  0.0;
    psF64  gmstRad      =  0.0;
    psF64  t            =  0.0;
    psF64  tu           =  0.0;
    psF64  ut1UtcDbl    =  0.0;
    psF64  const1       =  24110.5493771;
    psF64  const2       =  8639877.3173760;
    psF64  const3       =  307.4771600;
    psF64  const4       =  0.0931118;
    psF64  const5       = -0.0000062;
    psF64  const6       =  0.0000013;
    psTime *ut1UtcDelta = NULL;
    psTime *tdtTime     = NULL;
    psTime *taiTime     = NULL;
    psTime *utcTime     = NULL;
    psTime *ut1Time     = NULL;


    // Error checks
    PS_PTR_CHECK_NULL(time,NAN);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);

    // Calculate TAI or UTC time based on type of time user passes
    if(time->type == PS_TIME_TAI) {
        taiTime = psMemIncrRefCounter(time);
        utcTime = psTimeAlloc(PS_TIME_UTC);
        utcTime->sec = taiTime->sec - psTimeGetTAIDelta(time);
        utcTime->usec = taiTime->usec;
    } else if(time->type == PS_TIME_UTC) {
        utcTime = psMemIncrRefCounter(time);
        taiTime = psTimeAlloc(PS_TIME_TAI);
        taiTime->sec = utcTime->sec + psTimeGetTAIDelta(time);
        taiTime->usec = utcTime->usec;
    }

    // Convert Universal Time (UTC) to  UT1
    ut1UtcDbl = psTimeGetUT1Delta(taiTime);
    utcTime->type = PS_TIME_TAI; // Don't allow addition of leapseconds to ut1Time
    ut1Time = psTimeMath(utcTime, ut1UtcDbl);
    utcTime->type = PS_TIME_UTC;


    // Calculate UT1 as Julian Centuries since J2000.0
    jdUt1Days = psTimeToJD(ut1Time);
    mjdUt1Days = psTimeToMJD(ut1Time);
    t = (jdUt1Days - 2451545.0)/36525.0;

    // Calculate Terrestial Dynamical Time (TDT)
    tdtTime = psTimeMath(taiTime, 32.184);

    // Calculate TDT as Julian centuries since J2000.0
    jdTdtDays = psTimeToJD(tdtTime);
    tu = (jdTdtDays - 2451545.0)/36525.0;

    // Calculate fractional part of MJD
    fracDays = fmod(mjdUt1Days, 1.0);

    // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications.
    gmstRad = fracDays*TWOPI+(const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R;

    // Place GMST between 0 and 2*pi
    gmstRad = fmod(gmstRad, TWOPI);

    // Calculate Local Mean Sidereal Time (LMST) in radians
    lmstRad = gmstRad + longitude;

    // Free temporary structs
    psFree(ut1UtcDelta);
    psFree(ut1Time);
    psFree(tdtTime);
    psFree(utcTime);
    psFree(taiTime);

    return lmstRad;
}

double psTimeGetUT1Delta(const psTime *time)
{
    psU32 nTables = 3;
    psF64 mjd = 0.0;
    psF64 result = 0.0;
    psF64 dut2ut1 = 0.0;
    psF64 t = 0.0;
    psVector *dut = NULL;
    psMetadataItem *tableMetadataItem = NULL;
    psLookupStatusType status = PS_LOOKUP_SUCCESS;
    char *metadataTableNames[3] = {"ser7Table", "eopcTable",  "finalsTable"};


    // Error checks
    PS_PTR_CHECK_NULL(time,NAN);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);

    if(time->type != PS_TIME_TAI) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type);
        return 0.0;
    }

    // Attempt to find value through table lookup and interpolation
    mjd = psTimeToMJD(time);
    result = searchTables(mjd, 0, &status, metadataTableNames, nTables);

    // Value could not be found through table lookup and interpolation
    if(status == PS_LOOKUP_PAST_TOP) {

        // Date too early for tables. Get default time delta value from metadata, and issue warning.
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "UT1-UTC");

        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.dut");
            return 0.0;
        }
        result = tableMetadataItem->data.F64;

    } else if(status == PS_LOOKUP_PAST_BOTTOM) {
        /* Date too late for tables. Issue warning and use following formulae for predicting
           ahead of the most recent available table entry.
             ut1-utc = [0] + [1]*(MJD - [2]) - (ut2-ut1)
             [0, 1, 2] = @psLib.time.predict.dut
             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)
             t = 2000.0 + (MJD - 51544.03)/365.2422
        */
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC");

        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.dut");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.dut");
            return 0.0;
        }
        dut = (psVector*)tableMetadataItem->data.V;
        PS_PTR_CHECK_NULL(dut,0.0);

        t = 2000.0 + (mjd - 51544.03)/365.2422;
        dut2ut1 = 0.022*sin(TWOPI*t) - 0.012*cos(TWOPI*t) - 0.006*sin(4.0*M_PI*t) + 0.007*cos(4.0*M_PI*t);
        result = dut->data.F64[0] + dut->data.F64[1]*(mjd - dut->data.F64[2]) - dut2ut1;

    } else if(status != PS_LOOKUP_SUCCESS) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
        return 0.0;
    }

    return result;
}

struct psSphere* psTimeGetPoleCoords(const psTime* time)
{
    psU32 nTables = 3;
    psF64 x = 0.0;
    psF64 y = 0.0;
    psF64 mjd = 0.0;
    psF64 a = 0.0;
    psF64 c = 0.0;
    psF64 mjdPred = 0.0;
    struct psSphere* output = NULL;
    psLookupStatusType xStatus = PS_LOOKUP_SUCCESS;
    psLookupStatusType yStatus = PS_LOOKUP_SUCCESS;
    psMetadataItem *tableMetadataItem = NULL;
    char *metadataTableNames[3] = {"ser7Table", "eopcTable",  "finalsTable"};
    psVector *xp = NULL;
    psVector *yp = NULL;

    // Error checks
    PS_PTR_CHECK_NULL(time,NULL);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);

    if(time->type != PS_TIME_TAI)
    {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type);
        return NULL;
    }

    // Attempt to find value through table lookup and interpolation
    mjd = psTimeToMJD(time);
    x = searchTables(mjd, 0, &xStatus, metadataTableNames, nTables);
    y = searchTables(mjd, 0, &yStatus, metadataTableNames, nTables);

    // Value could not be found through table lookup and interpolation
    if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP)
    {

        // Date too earlier for tables. Get default polar coodinate values from metadata, and issue warning.
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "polar motion");

        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.xp");
            return NULL;
        }
        x = tableMetadataItem->data.F64;

        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.yp");
            return NULL;
        }
        y = tableMetadataItem->data.F64;

    } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM)
    {

        /* Date too late for tables. Issue warning and use following formulae for predicting
           ahead of the most recent available table entry.
              x = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c
              [0], [1], [2], [3] = @psLib.time.predict.xp
              y = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c
              [0], [1], [2], [3] = @psLib.time.predict.yp
              a = 2*pi*(mjd - pslib.time.predict.mjd)/365.25
              c = 2*pi*(mjd - pslib.time.predict.mjd)/435.0
        */
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "polar motion");


        // Get predicted MJD
        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                    "psLib.time.predict.mjd");
            return NULL;
        }
        mjdPred = tableMetadataItem->data.F64;

        // Get xp
        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp");
            return NULL;
        }
        xp = (psVector*)tableMetadataItem->data.V;
        PS_PTR_CHECK_NULL(xp,NULL);

        // Get yp
        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp");
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp");
            return NULL;
        }
        yp = (psVector*)tableMetadataItem->data.V;
        PS_PTR_CHECK_NULL(yp,NULL);

        // Calculate "a" and "c" constants
        a = TWOPI*(mjd - mjdPred)/365.25;
        c = TWOPI*(mjd - mjdPred)/435.0;

        // Calculate x and y polar coordinates
        x = xp->data.F64[0] +
            xp->data.F64[1]*cos(a) +
            xp->data.F64[2]*sin(a) +
            xp->data.F64[3]*cos(c) +
            xp->data.F64[4]*sin(c);

        y = yp->data.F64[0] +
            yp->data.F64[1]*cos(a) +
            yp->data.F64[2]*sin(a) +
            yp->data.F64[3]*cos(c) +
            yp->data.F64[4]*sin(c);

    } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS)
    {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
        return NULL;
    }

    // Create output sphere and convert arcsec to radians (i.e. x/60/60*M_PI/180)
    output = psAlloc(sizeof(psSphere));
    output->r = x * M_PI / 648000.0;
    output->d = y * M_PI / 648000.0;

    return output;
}

double psTimeGetTAIDelta(const psTime *time)
{
    psU64 i = 0;
    psF64 jd = 0.0;
    psF64 mjd = 0.0;
    psF64 out = 0.0;
    psF64 const1 = 0.0;
    psF64 const2 = 0.0;
    psF64 const3 = 0.0;
    psLookupTable* table = NULL;
    psMetadataItem *tableMetadataItem = NULL;
    psVector *stats = NULL;
    psVector *results = NULL;

    // Error checks
    PS_PTR_CHECK_NULL(time,NAN);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);

    // Check time metadata
    if(timeMetadata == NULL) {
        if(!p_psTimeInit(TIME_CONFIG_FILE)) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND, "psTime.config");
            return 0.0;
        }
    }

    // Get table from metadata
    tableMetadataItem = psMetadataLookup(timeMetadata, "taiTable");
    if(tableMetadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "taiTable");
        return 0.0;
    }
    table = (psLookupTable*)tableMetadataItem->data.V;
    PS_PTR_CHECK_NULL(table,0);

    // Determine Julian and modified Julian dates used in table lookup and time delta calculation
    jd = psTimeToJD(time);
    mjd = psTimeToMJD(time);

    stats = psVectorAlloc(table->numCols, PS_TYPE_U32);
    results = psLookupTableInterpolateAll(table, jd, stats);

    // Check for successful interpolation
    for(i=0; i<stats->n; i++) {
        if(stats->data.U32[i] == PS_LOOKUP_ERROR) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
        }
    }

    const1 = results->data.F64[0];
    const2 = results->data.F64[1];
    const3 = results->data.F64[2];
    out = const1 + (mjd - const2) * const3;

    psFree(stats);
    psFree(results);

    return out;
}


psS64 psTimeLeapseconds(const psTime *time1, const psTime *time2)
{
    psS64 diff = 0;


    // Error checks
    PS_PTR_CHECK_NULL(time1,0);
    PS_PTR_CHECK_NULL(time2,0);
    PS_INT_CHECK_RANGE(time1->usec,0,1e6-1,0);
    PS_INT_CHECK_RANGE(time2->usec,0,1e6-1,0);
    diff = abs((psS64)psTimeGetTAIDelta((psTime*)time1)-(psS64)psTimeGetTAIDelta((psTime*)time2));

    return diff;
}

double psTimeToJD(const psTime *time)
{
    double jd = 0.0;

    // Error checks
    PS_PTR_CHECK_NULL(time,NAN);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);

    // Julian date conversion
    if(time->sec < 0) {
        jd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 2440587.5; // psTime earlier than epoch
    } else {
        jd = time->sec / SEC_PER_DAY + time->usec / USEC_PER_DAY + 2440587.5; // psTime greater than epoch
    }

    return jd;
}

double psTimeToMJD(const psTime *time)
{
    double mjd = 0.0;


    // Error checks
    PS_PTR_CHECK_NULL(time,NAN);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);

    // Modified Julian date conversion
    if(time->sec < 0) {
        mjd = time->sec / SEC_PER_DAY - time->usec / USEC_PER_DAY + 40587.0; // psTime earlier than epoch
    } else {
        mjd = time->sec / SEC_PER_DAY + time->usec / USEC_PER_DAY + 40587.0; // psTime greater than epoch
    }

    return mjd;
}

char* psTimeToISOTime(const psTime *time)
{
    psS32 ms = 0;
    char *timeString = NULL;
    char *tempString = NULL;
    struct tm *tmTime = NULL;
    time_t sec;


    // Error checks
    PS_PTR_CHECK_NULL(time,NULL);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);

    tempString = psAlloc(MAX_TIME_STRING_LENGTH);
    timeString = psAlloc(MAX_TIME_STRING_LENGTH);

    ms = time->usec / 1000;
    sec = time->sec;

    // tmTime variable is statically allocated, no need to free
    tmTime = gmtime(&sec);

    // Converts psTime to YYYY-MM-DDThh:mm:ss.sss in string form
    if (!strftime(tempString, MAX_TIME_STRING_LENGTH, "%Y-%m-%dT%H:%M:%S", tmTime)) {
        psError(PS_ERR_OS_CALL_FAILED, true, PS_ERRORTEXT_psTime_CONVERT_TIME_TO_STRING_FAILED);
    }

    if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s.%3.3dZ", tempString, ms) < 0) {
        psError(PS_ERR_OS_CALL_FAILED, true, PS_ERRORTEXT_psTime_APPEND_MSEC_FAILED);
    }
    psFree(tempString);

    return timeString;
}

struct timeval psTimeToTimeval(const psTime *time)
{
    struct timeval timevalTime;

    timevalTime.tv_sec = 0;
    timevalTime.tv_usec = 0;

    // Error checks
    PS_PTR_CHECK_NULL(time,timevalTime);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,timevalTime);

    timevalTime.tv_sec = time->sec;
    timevalTime.tv_usec = time->usec;

    return timevalTime;
}

struct tm* psTimeToTM(const psTime *time)
{
    psS64 cent = 0;
    psS64 year = 0;
    psS64 month = 0;
    psS64 day = 0;
    psS64 hour = 0;
    psS64 minute = 0;
    psS64 seconds = 0;
    psS64 temp = 0;
    struct tm* tmTime = NULL;


    // Error checks
    PS_PTR_CHECK_NULL(time,NULL);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);

    seconds = time->sec%60;
    minute = time->sec/60%60;
    hour = time->sec/3600%24;
    day = (time->sec+62135596800)/86400;

    // Add 306 days to make relative to Mar 1, 0; also adjust day to be within a range (1..2**28-1) where our
    // calculations will work with 32bit ints
    if(day > (pow(2, 28)-307))
    {
        temp = (day - 146097+306)/146097+1;         // Avoid overflow if day close to maxint
        day -= temp * 146097-306;
    } else if((day += 306) <= 0)
    {
        temp = -( -day / 146097 + 1);               // Avoid ambiguity in C division of negatives
        day -= temp * 146097;
    }

    cent = (day*4-1)/146097;                        // Calc number of centuries day is after 29 Feb of yr 0
    day -= cent*146097/4;                           // 4 centuries = 146097 days
    year = (day*4-1)/1461;                          // Calc number of years into the century
    day -= year*1461/4;                             // Again March-based (4 yrs =\u02dc 146[01] days)
    month = (day*12+1093)/367;                      // Get the month (3..14 represent March through
    day -= (month*367-1094)/12;                     // February of following year)
    year += cent*100+temp*400;                      // Get the real year, which is off by

    // One if month is January or February
    if(month > 12)
    {
        year++;
        month -= 12;
    }

    // Allocate output
    tmTime = (struct tm*)psAlloc(sizeof(struct tm));

    tmTime->tm_year = year - 1900;
    tmTime->tm_mon = month - 1;
    tmTime->tm_mday = day + 1;
    tmTime->tm_hour = hour;
    tmTime->tm_min = minute;
    tmTime->tm_sec = seconds;
    tmTime->tm_isdst = -1;

    return tmTime;
}

psTime* psTimeFromJD(double time)
{
    double days = 0.0;
    double seconds = 0.0;
    psTime *outTime = NULL;


    // Allocate psTime struct
    outTime = psTimeAlloc(PS_TIME_TAI);

    // Julian date conversion courtesy of Eugene Magnier
    days = time - 2440587.5;
    seconds = days * SEC_PER_DAY;
    if(seconds < 0.0) {
        outTime->usec = (seconds - (psS64)seconds) * -1000000.0;  // psTime earlier than epoch
    } else {
        outTime->usec = (seconds - (psS64)seconds) * 1000000.0;   // psTime greater than epoch
    }
    outTime->sec = seconds;

    // Error check
    PS_INT_CHECK_RANGE(outTime->usec,0,1e6-1,outTime);

    return outTime;
}

psTime* psTimeFromMJD(double time)
{
    double days = 0.0;
    double seconds = 0.0;
    psTime *outTime = NULL;


    // Allocate psTime struct
    outTime = psTimeAlloc(PS_TIME_TAI);

    // Modified Julian date conversion courtesy of Eugene Magnier
    days = time - 40587.0;
    seconds = days * SEC_PER_DAY;
    if(seconds < 0.0) {
        outTime->usec = (seconds - (psS64)seconds) * -1000000.0;  // psTime earlier than epoch
    } else {
        outTime->usec = (seconds - (psS64)seconds) * 1000000.0;   // psTime greater than epoch
    }
    outTime->sec = seconds;

    // Error check
    PS_INT_CHECK_RANGE(outTime->usec,0,1e6-1,outTime);

    return outTime;
}

psTime* psTimeFromISOTime(const char *time)
{
    psS32 millisecond;
    struct tm tmTime;
    psTime *outTime = NULL;

    // Convert YYYY-MM-DDThh:mm:ss.sss in string form to tm time
    if (sscanf(time, "%d-%d-%dT%d:%d:%d.%d", &tmTime.tm_year, &tmTime.tm_mon, &tmTime.tm_mday,
               &tmTime.tm_hour, &tmTime.tm_min, &tmTime.tm_sec,&millisecond) < 7) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_ISOTIME_MALFORMED, time);
        return NULL;
    }

    PS_INT_CHECK_NON_NEGATIVE(tmTime.tm_year, outTime);
    PS_INT_CHECK_RANGE(tmTime.tm_mon,1,12,outTime);
    PS_INT_CHECK_RANGE(tmTime.tm_mday,1,31,outTime);
    PS_INT_CHECK_RANGE(tmTime.tm_hour,0,23,outTime);
    PS_INT_CHECK_RANGE(tmTime.tm_min,0,59,outTime);
    PS_INT_CHECK_RANGE(tmTime.tm_sec,0,59,outTime);
    PS_INT_CHECK_RANGE(millisecond,0,999,outTime);

    tmTime.tm_year -= 1900;
    tmTime.tm_mon--;
    tmTime.tm_isdst = -1;

    // Convert tm time to psTime
    outTime = psTimeFromTM(&tmTime);
    outTime->usec = millisecond * 1000;

    return outTime;
}

psTime* psTimeFromTimeval(const struct timeval *time)
{
    psTime *outTime = NULL;


    // Error check
    PS_PTR_CHECK_NULL(time,NULL);

    // Allocate psTime struct
    outTime = psTimeAlloc(PS_TIME_TAI);

    // Convert to psTime
    outTime->sec = time->tv_sec;
    outTime->usec = time->tv_usec;

    // Error check
    PS_INT_CHECK_RANGE(outTime->usec,0,1e6-1,outTime);

    return outTime;
}

psTime* psTimeFromTM(const struct tm* time)
{
    psS64 year;
    psS64 month;
    psS64 day;
    psS64 hour;
    psS64 minute;
    psS64 seconds;
    psS64 temp;
    psTime *outTime = NULL;


    // Error check
    PS_PTR_CHECK_NULL(time,NULL);

    // Allocate psTime struct
    outTime = psTimeAlloc(PS_TIME_TAI);

    // Extract data from TM struct
    year = time->tm_year + 1900;
    month = time->tm_mon + 1;
    day = time->tm_mday;
    hour = time->tm_hour;
    minute = time->tm_min;
    seconds = time->tm_sec;

    // Make month in range 3..14 (treat Jan & Feb as months 13..14 of prev year)
    if( month <= 2 )
    {
        year -= (temp = (14 - month) / 12);
        month += 12 * temp;
    } else if(month > 14)
    {
        year += (temp = (month - 3) / 12);
        month -= 12 * temp;
    }

    // Make year positive
    if (year < 0 )
    {
        day -= 146097 * (temp = (399 - year) / 400);
        year += 400 * temp;
    }

    // Add day of month, days of previous 0-11 month period that began w/March, days of previous 0-399 year
    // period that began w/March of a 400-multiple year), days of any 400-year periods before that, and 306
    // days to adjust from Mar 1, year 0-relative to Jan 1, year 1-relative. Add hours, minutes, and seconds.
    day += (month * 367 - 1094) / 12 + year % 100 * 1461 / 4 + (year/100 * 36524 + year/400) - 306;
    outTime->sec = (((day - 1) * SEC_PER_DAY) - 62135596800) + hour*SEC_PER_HOUR + minute*SEC_PER_MINUTE + seconds;

    // C's TM does not define a microsecond field. Microseconds must be manipulated by calling function.
    outTime->usec = 0;

    // Error check
    PS_INT_CHECK_RANGE(outTime->usec,0,1e6-1,outTime);

    return outTime;
}

psTime* psTimeMath(const psTime *time, psF64 delta)
{
    psF64 sec = 0.0;
    psTime *outTime = NULL;
    psTime *tempTime = NULL;


    // Error checks
    PS_PTR_CHECK_NULL(time,NULL);
    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);

    // Convert time to TAI if necessary, but without changing input arguments
    if(time->type == PS_TIME_UTC) {
        tempTime = psTimeAlloc(PS_TIME_UTC);
        tempTime->sec = time->sec;
        tempTime->usec = time->usec;
        tempTime = psTimeConvert(tempTime, PS_TIME_TAI);
    } else {
        tempTime = psMemIncrRefCounter((psTime*)time);
    }

    // Create output time
    outTime = psTimeAlloc(PS_TIME_TAI);
    sec = delta + (psF64)tempTime->sec + (psF64)tempTime->usec/1e6;
    outTime->sec = sec;
    outTime->usec = (sec - outTime->sec)*1e6;

    // Error check
    PS_INT_CHECK_RANGE(outTime->usec,0,1e6-1,outTime);

    // Convert result to same time type as input
    if(time->type == PS_TIME_UTC) {
        outTime = psTimeConvert(outTime, PS_TIME_UTC);
    }

    psFree(tempTime);

    return outTime;
}

psF64 psTimeDelta(const psTime *time1, const psTime *time2)
{
    psF64 out = 0.0;
    psF64 uSec1 = 0.0;
    psF64 uSec2 = 0.0;
    psTime *tempTime1 = NULL;
    psTime *tempTime2 = NULL;


    // Error checks
    PS_PTR_CHECK_NULL(time1,0.0);
    PS_INT_CHECK_RANGE(time1->usec,0,1e6-1,0.0);
    PS_PTR_CHECK_NULL(time2,0.0);
    PS_INT_CHECK_RANGE(time2->usec,0,1e6-1,0.0);

    // Convert time to TAI if necessary, but without changing input arguments
    if(time1->type == PS_TIME_UTC) {
        tempTime1 = psTimeAlloc(PS_TIME_UTC);
        tempTime1->sec = time1->sec;
        tempTime1->usec = time1->usec;
        tempTime1 = psTimeConvert(tempTime1, PS_TIME_TAI);
    } else {
        tempTime1 = psMemIncrRefCounter((psTime*)time1);
    }
    if(time2->type == PS_TIME_UTC) {
        tempTime2 = psTimeAlloc(PS_TIME_UTC);
        tempTime2->sec = time2->sec;
        tempTime2->usec = time2->usec;
        tempTime2 = psTimeConvert(tempTime2, PS_TIME_TAI);
    } else {
        tempTime2 = psMemIncrRefCounter((psTime*)time2);
    }

    uSec1 = tempTime1->sec >= 0 ? 1.0 : -1.0;
    uSec1 = uSec1*tempTime1->usec/1e6;
    uSec2 = tempTime2->sec >= 0 ? 1.0 : -1.0;
    uSec2 = uSec2*tempTime2->usec/1e6;
    out = (tempTime1->sec-tempTime2->sec) + (uSec1-uSec2);

    psFree(tempTime1);
    psFree(tempTime2);

    return out;
}

