/** @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.74 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-10-05 03:51:41 $
 *
 *  Copyright 2004-2005 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 "psLogMsg.h"
#include "psMemory.h"
#include "psAbort.h"
#include "psImage.h"
#include "psString.h"
#include "psMetadata.h"
#include "psMetadataConfig.h"
#include "psLookupTable.h"
#include "psConstants.h"
#include "psErrorText.h"

#include "config.h"

#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 NSEC_PER_DAY 86400000000000.0

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

// Offset to convert terrestrial time(TT) to international atomic time(TAI)
#define  TAI_TT_OFFSET_SECONDS        32
#define  TAI_TT_OFFSET_NANOSECONDS    184000000

// Offset from converting to MJD
#define  MJD_EPOCH_OFFSET             40587.0

// Offset for converting to JD
#define  JD_EPOCH_OFFSET              2440587.5

// Offset of year 0000 from epoch
#define YEAR_0000_SEC                 -62125920000.0

// Offset of year 9999 from epoch
#define YEAR_9999_SEC                 253202544000.0

/** Static function prototypes */
static char *cleanString(char *inString, int sLen);
static char* getToken(char **inString, char *delimiter, psParseErrorType *status);
static psTime* convertTimeTAIUTC(psTime* time);
static psTime* convertTimeUTCTAI(psTime* time);
static psTime* convertTimeTAITT(psTime* time);
static psTime* convertTimeTTTAI(psTime* time);
static psTime* convertTimeUTCUT1(psTime* time);

/** 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 - 1;
    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+1);

    } else if(**inString!='\0' && sLen==0) {
        *status = PS_PARSE_ERROR_GENERAL;
    }

    return cleanToken;
}

// get the pslib.config filename by checking environment variable first, then original installation area.
char* p_psGetConfigFileName()
{
    char* filename = getenv("PS_CONFIG_FILE");

    if (filename == NULL) { // environment variable not found
        filename = PS_CONFIG_FILE_DEFAULT; // this should come from configure.ac
    }

    return filename;
}


// 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 p_psTimeSearchTables(psF64 index, psU64 column, char *metadataTableNames[],
                           psU32 nTables, psLookupStatusType* status)
{
    char*            tableName          = NULL;
    psF64            result             = NAN;
    psLookupTable*   table              = NULL;
    psMetadataItem*  tableMetadataItem  = NULL;

    // Check if psTime tables are already loaded
    if(!p_psTimeInit(p_psGetConfigFileName())) {
        *status = PS_LOOKUP_ERROR;
        return NAN;
    }

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

        // Get table name from list of tables to search
        tableName = metadataTableNames[i];

        // Lookup table name in time metadata
        tableMetadataItem = psMetadataLookup(timeMetadata, tableName);

        // Check if table not a metadata item
        if(tableMetadataItem == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                    tableName);
            *status = PS_LOOKUP_ERROR;
            return NAN;
        }

        // Get table from metadata
        table = (psLookupTable*)tableMetadataItem->data.V;

        // Check that table is not NULL
        PS_ASSERT_PTR_NON_NULL(table,NAN);

        // Check if index within to/from range
        if(index >= table->validFrom ) {
            if(index <= table->validTo) {
                // Attempt to interpolate table
                result = psLookupTableInterpolate(table, index, column);
                *status = PS_LOOKUP_SUCCESS;
                if(!isnan(result)) {
                    break;
                }
            } else {
                *status = PS_LOOKUP_PAST_BOTTOM;
            }
        } else {
            *status = PS_LOOKUP_PAST_TOP;
        }
    }

    return result;
}

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

    // Check if the p_psTimeInit has already been ran
    if (timeMetadata != NULL) {
        return true;
    }

    // Read config file
    timeMetadata = psMetadataConfigParse(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");
        psFree(tablesFrom);
        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 index columns for the tables
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.index");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.index");
        psFree(tablesFrom);
        psFree(tablesTo);
        return false;
    }
    tablesIndex = psVectorCopy(tablesIndex, metadataItem->data.V, PS_TYPE_S32);
    if(tablesIndex->n != numTables) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_BAD_VECTOR,tablesIndex->n,numTables);
        psFree(tablesFrom);
        psFree(tablesTo);
        psFree(tablesIndex);
        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);
        psFree(tablesIndex);
        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(tablesIndex);
        psFree(tableDir);
        return false;
    }
    tableNames = psStringCopy(metadataItem->data.V);

    // Get table format strings
    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.format");
    if(metadataItem == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.time.tables.format");
        psFree(tablesFrom);
        psFree(tablesTo);
        psFree(tablesIndex);
        psFree(tableDir);
        psFree(tableNames);
        return false;
    }
    tableFormats = psStringCopy(metadataItem->data.V);
    formatPtr = tableFormats;

    // 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);

        // Get table format
        tableFormat = getToken(&formatPtr,",",&status);
        if(tableFormat == NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                    "psLib.time.tables.format");
        }

        // Create and read table
        if(i < numTables) {
            table = psLookupTableAlloc(fullTableName, (const char*)tableFormat, tablesIndex->data.S32[i]);
            numLines = 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_DATA_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);
        psFree(tableFormat);
        psFree(table);
        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);
    psFree(tablesIndex);
    psFree(tableFormats);

    return true;
}

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

    return true;
}

static void timeFree(psTime *outTime)
{
    // There are non dynamic allocated items
}

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

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

    // Allocate memory for structure
    outTime = (psTime*)psAlloc(sizeof(psTime));
    psMemSetDeallocator(outTime, (psFreeFunc)timeFree);
    // Initialize members
    outTime->sec = 0;
    outTime->nsec = 0;
    outTime->type = type;
    outTime->leapsecond = false;

    return outTime;
}


bool psMemCheckTime(psPtr ptr)
{
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)timeFree );
}


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

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

    // Verify time structure allocated
    if(time == NULL) {
        return NULL;
    }

    // Get the system time
    //    if (gettimeofday(&now, (struct timezone *)0) == -1) {
    if (gettimeofday(&now, 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->nsec = now.tv_usec*1000;

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

    return time;
}

static psTime* convertTimeTAIUTC(psTime* time)
{
    psF64  deltaTAI     = 0.0;
    psS64  deltaSec     = 0;
    psU32  deltaNsec    = 0;
    psF64  deltaUTC     = 0.0;

    // Determine delta to convert between UTC and TAI
    deltaTAI = p_psTimeGetTAIDelta(time);
    deltaSec = (psS64)(deltaTAI);
    deltaNsec = (psU32)((deltaTAI - (psF64)deltaSec) * 1e9);

    // Determine seconds
    time->sec -= deltaSec;

    // Check for underflow in nsec
    if(deltaNsec > time->nsec) {
        // Borrow second
        time->nsec += 1e9;
        time->sec--;
    }

    // Determine nsec
    time->nsec -= deltaNsec;

    // Check for overflow in nsec
    if(time->nsec >= 1e9) {
        time->nsec -= 1e9;
        time->sec++;
    }

    // Set new type
    time->type = PS_TIME_UTC;

    // Check if leapsecond present in delta
    deltaUTC = p_psTimeGetTAIDelta(time);
    if(fabs(deltaTAI-deltaUTC) >= 1.0) {
        time->sec++;
    }

    return time;
}

static psTime* convertTimeUTCTAI(psTime* time)
{
    psF64  delta     = 0.0;
    psS64  deltaSec  = 0;
    psU32  deltaNsec = 0;

    // Determine delta to convert between UTC and TAI
    delta = p_psTimeGetTAIDelta(time);

    deltaSec = (psS64)(delta);
    deltaNsec = (psU32)((delta - (psF64)deltaSec) * 1e9);

    // Determine seconds
    time->sec += deltaSec;

    // Determine nsec
    time->nsec += deltaNsec;

    // Check for overflow in nsec
    if(time->nsec >= 1e9) {
        time->nsec -= 1e9;
        time->sec++;
    }

    // Set new type
    time->type = PS_TIME_TAI;

    return time;
}

static psTime* convertTimeTAITT(psTime* time)
{
    // Add TT offset
    time->sec += TAI_TT_OFFSET_SECONDS;
    time->nsec += TAI_TT_OFFSET_NANOSECONDS;

    // Check for overflow in nsec
    if(time->nsec >= 1e9) {
        time->nsec -= 1e9;
        time->sec++;
    }

    // Set new type
    time->type = PS_TIME_TT;

    return time;
}

static psTime* convertTimeTTTAI(psTime* time)
{
    // Subtract TT offset
    time->sec -= TAI_TT_OFFSET_SECONDS;

    // Check for nsec underflow
    if(TAI_TT_OFFSET_NANOSECONDS > time->nsec) {
        // Borrow second
        time->sec--;
        time->nsec += 1e9;
    }
    time->nsec -= TAI_TT_OFFSET_NANOSECONDS;

    // Check for overflow in nsec
    if(time->nsec >= 1e9) {
        time->nsec -= 1e9;
        time->sec++;
    }

    // Set new type
    time->type = PS_TIME_TAI;

    return time;
}

static psTime* convertTimeUTCUT1(psTime* time)
{
    psS64   ut1utc  = 0;

    // Get UT1-UTC value
    ut1utc = (psS64)(psTimeGetUT1Delta(time,PS_IERS_A) * 1e9);

    // Since UTC is within 0.9 sec of UT1 then nsec member is the member affected
    if((ut1utc < 0) && (abs(ut1utc) > time->nsec)) {
        // Borrow from sec
        time->sec--;
        if(time->leapsecond) {
            time->leapsecond = false;
        } else {
            time->leapsecond = psTimeIsLeapSecond(time);
        }
        // Add to nsec
        time->nsec += 1e9;
    }
    time->nsec += ut1utc;

    // Check for overflow in nsec
    if(time->nsec >= 1e9) {
        time->nsec -= 1e9;
        time->sec++;
        if(time->leapsecond) {
            time->leapsecond = false;
            time->sec--;
        } else {
            time->leapsecond = psTimeIsLeapSecond(time);
        }
    }

    // Set new type
    time->type = PS_TIME_UT1;

    return time;
}

psTime* psTimeConvert(psTime *time, psTimeType type)
{
    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),time);

    // If the input type is UT1 then return time and generate error message
    if(time->type == PS_TIME_UT1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Cannot convert from UT1 time type");
        return time;
    }

    // If the time to convert to is the same as psTime the return time
    if (time->type == type) {
        return time;
    }

    // Convert from TAI to UTC, TT, UT1
    if(time->type == PS_TIME_TAI) {
        // Convert from TAI to UTC
        if(type == PS_TIME_UTC) {
            time = convertTimeTAIUTC(time);
            // Convert from TAI to TT
        } else if(type == PS_TIME_TT) {
            time = convertTimeTAITT(time);
            // Convert from TAI to UT1
        } else if(type == PS_TIME_UT1) {
            // Convert to UTC first
            time = convertTimeTAIUTC(time);
            // Convert UTC to UT1
            time = convertTimeUTCUT1(time);
            // Convert from TAI to unknown time type
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type);
        }
        // Convert from TT to TAI, UTC, UT1
    } else if(time->type == PS_TIME_TT) {
        // Convert from TT to UTC
        if(type == PS_TIME_UTC) {
            // Convert to TAI time first
            time = convertTimeTTTAI(time);
            // Convert from TAI to UTC
            time = convertTimeTAIUTC(time);
            // Convert from TT to TAI
        } else if(type == PS_TIME_TAI) {
            time = convertTimeTTTAI(time);
            // Convert from TT to UT1
        } else if(type == PS_TIME_UT1) {
            // Convert to UTC first
            // Convert to TAI time first
            time = convertTimeTTTAI(time);
            // Convert from TAI to UTC
            time = convertTimeTAIUTC(time);
            // Convert from UTC to UT1
            time = convertTimeUTCUT1(time);
            // Convert from TT to unknown time type
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type);
        }
        // Convert from UTC to TAI, TT, UT1
    } else if(time->type == PS_TIME_UTC) {
        // Convert UTC to TAI
        if(type == PS_TIME_TAI) {
            time = convertTimeUTCTAI(time);
            // Convert UTC to TT
        } else if(type == PS_TIME_TT) {
            // Convert to TAI time first
            time = convertTimeUTCTAI(time);
            // Convert TAI to TT
            time = convertTimeTAITT(time);
            // Convert UTC to UT1
        } else if(type == PS_TIME_UT1) {
            time = convertTimeUTCUT1(time);
            // Convert UTC to unknown time type
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type);
        }
        // Convert unknown time type
    } else {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, time->type);
    }

    return time;
}

double psTimeToLMST(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  const1       =  24110.5493771;
    psF64  const2       =  8639877.3173760;
    psF64  const3       =  307.4771600;
    psF64  const4       =  0.0931118;
    psF64  const5       = -0.0000062;
    psF64  const6       =  0.0000013;
    psTime *tdtTime     = NULL;
    psTime *ut1Time     = NULL;

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN);

    // Verify input time is not in UT1 seconds
    if(time->type == PS_TIME_UT1) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time->type);
        return NAN;
    }

    // Determine time reference to TT
    tdtTime = psTimeAlloc(time->type);
    tdtTime->sec = time->sec;
    tdtTime->nsec = time->nsec;
    tdtTime->leapsecond = time->leapsecond;
    tdtTime = psTimeConvert(tdtTime,PS_TIME_TT);

    // Determine time reference to UT1
    ut1Time = psTimeAlloc(time->type);
    ut1Time->sec = time->sec;
    ut1Time->nsec = time->nsec;
    ut1Time->leapsecond = time->leapsecond;
    ut1Time = psTimeConvert(ut1Time,PS_TIME_UT1);

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

    // 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(ut1Time);
    psFree(tdtTime);

    return lmstRad;
}

double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin)
{
    psU32              nTables               = 2;
    psF64              mjd                   = 0.0;
    psF64              result                = 0.0;
    psU64              tableColumn           = 0;
    psF64              dut2ut1               = 0.0;
    psF64              t                     = 0.0;
    psVector*          dut                   = NULL;
    psMetadataItem*    tableMetadataItem     = NULL;
    psLookupStatusType status                = PS_LOOKUP_SUCCESS;
    char*              metadataTableNames[2] = {"dailyTable",  "finalsTable"};

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN);

    // Check for invalid bulletin specified
    if((bulletin != PS_IERS_A) && (bulletin != PS_IERS_B)) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Invalid bulletin specified %d",bulletin);
        return NAN;
    }

    // Set lookup table column based on Bullentin
    if(bulletin == PS_IERS_A) {
        tableColumn = 3;
    } else {
        tableColumn = 6;
    }

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

    // 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.
        psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"UT1-UTC");

        // Lookup value from time metadata loaded from pslib.config
        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 NAN;
        }
        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
        */
        // Generate warning of postdate information
        psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC");

        // Lookup values to calculate prediction
        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 NAN;
        }
        dut = (psVector*)tableMetadataItem->data.V;
        PS_ASSERT_PTR_NON_NULL(dut,NAN);

        // Calculate predication of future UT1-UTC
        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 NAN;
    }

    return result;
}

psSphere* p_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;
    psSphere* output = NULL;
    psLookupStatusType xStatus = PS_LOOKUP_SUCCESS;
    psLookupStatusType yStatus = PS_LOOKUP_SUCCESS;
    psMetadataItem *tableMetadataItem = NULL;
    char *metadataTableNames[3] = {"dailyTable", "eopcTable",  "finalsTable"};
    psVector *xp = NULL;
    psVector *yp = NULL;

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-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 = p_psTimeSearchTables(mjd, 0, &xStatus, metadataTableNames, nTables);
    x = p_psTimeSearchTables(mjd, 0, metadataTableNames, nTables,&xStatus);
    //    y = p_psTimeSearchTables(mjd, 0, &yStatus, metadataTableNames, nTables);
    y = p_psTimeSearchTables(mjd, 0, metadataTableNames, nTables,&yStatus);

    // 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_ASSERT_PTR_NON_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_ASSERT_PTR_NON_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*PS_PI/180)
    output = psAlloc(sizeof(psSphere));
    output->r = x * M_PI / 648000.0;
    output->d = y * M_PI / 648000.0;

    return output;
}

psF64 p_psTimeGetTAIDelta(const psTime *time)
{
    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 *results = NULL;

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN);

    // Check if psTime tables are loaded/loadable
    if(!p_psTimeInit(p_psGetConfigFileName())) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND,
                p_psGetConfigFileName());
        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_ASSERT_PTR_NON_NULL(table,0);

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

    // Set ceiling of the julian date to the last entry in the lookup table
    if(table->validTo < jd) {
        jd = table->validTo;
    }

    // Interpolation of look up table
    results = psLookupTableInterpolateAll(table, jd);

    // Check for successful interpolation
    if(results == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
    }

    // Set constants from table
    const1 = results->data.F64[1];
    const2 = results->data.F64[2];
    const3 = results->data.F64[3];

    // If const3 not equal to zero solve for difference else floor of const1
    if(fabs(const3-0.0) > FLT_EPSILON) {
        out = const1 + (mjd - const2) * const3;
    } else {
        out = floor(const1);
    }

    psFree(results);

    return out;
}

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

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time1,0);
    PS_ASSERT_PTR_NON_NULL(time2,0);
    PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0);
    PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0);
    diff = abs((psS64)p_psTimeGetTAIDelta((psTime*)time1)-(psS64)p_psTimeGetTAIDelta((psTime*)time2));

    return diff;
}

bool psTimeIsLeapSecond(const psTime* utc)
{
    psTime*    prevUtc     = NULL;
    psBool     returnValue = false;

    // Check for valid time
    PS_ASSERT_PTR_NON_NULL(utc,false);

    // Verify time is UTC type
    if(utc->type != PS_TIME_UTC) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,utc->type);
        return false;
    }

    // Allocate time to hold utc - 1 second
    prevUtc = psTimeAlloc(PS_TIME_UTC);
    prevUtc->sec = utc->sec - 1;

    // Check the absolute difference between the two times for leapsecond
    if(psTimeLeapSecondDelta(utc,prevUtc) >= 1.0) {
        returnValue = true;
    } else {
        returnValue = false;
    }

    // Free prevUtc
    psFree(prevUtc);

    return returnValue;
}

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

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN);

    // Julian date conversion
    if(time->sec < 0) {
        // psTime earlier than epoch
        jd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET;
    } else {
        // psTime greater than epoch
        jd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET;
    }
    return jd;
}

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

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NAN);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN);

    // Modified Julian date conversion
    if(time->sec < 0) {
        // psTime earlier than epoch
        mjd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET;
    } else {
        // psTime greater than epoch
        mjd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET;
    }
    return mjd;
}

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

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL);

    // Check valid year range
    PS_ASSERT_LONG_WITHIN_RANGE(time->sec,YEAR_0000_SEC,YEAR_9999_SEC,NULL)

    // Allocate temp strings
    tempString = psAlloc(MAX_TIME_STRING_LENGTH);
    timeString = psAlloc(MAX_TIME_STRING_LENGTH);

    // Convert nanoseconds to decaseconds
    ds = time->nsec / 100000000;
    sec = time->sec;

    // If leapsecond use previous day
    if(time->leapsecond) {
        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);
    }

    // Check if time is UTC and leapsecond
    if(((time->type==PS_TIME_UTC)||(time->type==PS_TIME_UT1)) && (time->leapsecond)) {
        // Modify second to be 60
        tempString[17] = '6';
        tempString[18] = '0';
    }

    // Create string with milliseconds
    if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s,%1dZ", tempString, ds) < 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 = NULL;

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->sec,0,INT32_MAX,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL);

    // Allocate structure timeval
    timevalTime = (struct timeval*)psAlloc(sizeof(struct timeval));

    // Set structure members
    timevalTime->tv_sec = time->sec;
    timevalTime->tv_usec = time->nsec / 1000;

    return timevalTime;
}

/*
struct tm* p_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_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-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 jd)
{
    psF64 days = 0.0;
    psF64 seconds = 0.0;
    psTime *outTime = NULL;

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

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

    // Error check
    PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime);

    return outTime;
}

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

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

    // Modified Julian date conversion courtesy of Eugene Magnier
    days = mjd - 40587.0;
    seconds = days * SEC_PER_DAY;

    if(seconds < 0.0) {
        outTime->nsec = (seconds - (psS64)seconds) * -1000000000.0;  // psTime earlier than epoch
    } else {
        outTime->nsec = (seconds - (psS64)seconds) * 1000000000.0;   // psTime greater than epoch
    }
    outTime->sec = seconds;

    // Error check
    PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),NULL);

    return outTime;
}

psTime* psTimeFromISO(const char *input,
                      psTimeType type)
{
    psS32 millisecond;
    struct tm tmTime;
    psTime *outTime = NULL;

    // Check for NULL string
    PS_ASSERT_PTR_NON_NULL(input,NULL);

    // Convert YYYY-MM-DDThh:mm:ss.sss in string form to tm time
    if (sscanf(input, "%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, input);
        return NULL;
    }

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

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

    // Convert tm time to psTime
    outTime = p_psTimeFromTM(&tmTime);
    outTime->nsec = millisecond * 1000000;
    outTime->type = type;

    return outTime;
}

psTime* psTimeFromTT(psS64 sec, psU32 nsec)
{
    psTime*      outTime  = NULL;

    // Verify nsec within range
    PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL);

    // Allocate psTime data
    outTime = psTimeAlloc(PS_TIME_TT);

    // Set data members
    outTime->sec = sec;
    outTime->nsec = nsec;

    // Return data structure
    return outTime;
}

psTime* psTimeFromUTC(psS64 sec, psU32 nsec, bool leapsecond)
{
    psTime*   outTime   = NULL;

    // Verify nsec within range
    PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL);

    // Allocate psTime data
    outTime = psTimeAlloc(PS_TIME_UTC);

    // Set data members
    outTime->sec = sec;
    outTime->nsec = nsec;

    // Set leapsecond flag if necessary
    outTime->leapsecond = psTimeIsLeapSecond(outTime);

    return outTime;
}

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


    // Error check
    PS_ASSERT_PTR_NON_NULL(input,NULL);

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

    // Convert to psTime
    outTime->sec = input->tv_sec;
    outTime->nsec = input->tv_usec * 1000;

    // Error check
    PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime);

    return outTime;
}

psTime* p_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_ASSERT_PTR_NON_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 )
    {
        temp = (14 - month) / 12;
        //        year -= (temp = (14 - month) / 12);
        year -= temp;
        month += 12 * temp;
    } else if(month > 14)
    {
        temp = (month - 3) / 12;
        //        year += (temp = (month - 3) / 12);
        year += temp;
        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->nsec = 0;

    // Error check
    PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime);

    return outTime;
}

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

    // Error checks
    PS_ASSERT_PTR_NON_NULL(time,NULL);
    PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-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->nsec = time->nsec;
        tempTime = psTimeConvert(tempTime, PS_TIME_TAI);
        outTime = psTimeAlloc(PS_TIME_TAI);
    } else {
        tempTime = psMemIncrRefCounter((psTime*)time);
        outTime = psTimeAlloc(time->type);
    }

    // Create output time
    sec = delta + (psF64)tempTime->sec + (psF64)tempTime->nsec/1e9;
    PS_ASSERT_LONG_WITHIN_RANGE((psS64)sec,0,PS_MAX_S64,outTime);
    outTime->sec = (psS64)sec;
    outTime->nsec = (psU32)((sec - (psF64)outTime->sec)*1e9);

    // Error check
    PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-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;
}

double 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_ASSERT_PTR_NON_NULL(time1,0.0);
    PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0.0);
    PS_ASSERT_PTR_NON_NULL(time2,0.0);
    PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0.0);

    // Verify both times of the same type
    if(time1->type != time2->type) {
        psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time1->type);
        return out;
    }

    // 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->nsec = time1->nsec;
        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->nsec = time2->nsec;
        tempTime2 = psTimeConvert(tempTime2, PS_TIME_TAI);
    } else {
        tempTime2 = psMemIncrRefCounter((psTime*)time2);
    }

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

    psFree(tempTime1);
    psFree(tempTime2);

    return out;
}

static psHash *timers = NULL;
static char timerName[256] = "\0";

static void psTimerFree ()
{

    psFree (timers);
    p_psTimeFinalize();
    return;
}

// start/restart a named timer
bool psTimerStart (char *name)
{

    psTime *start;

    if (timers == NULL) {
        timers = psHashAlloc (16);
        //        psMemSetDeallocator(timers, (psFreeFunc)psTimerFree);
    }
    strncpy(timerName, name, 256);
    start = psTimeGetNow (PS_TIME_UTC);
    psHashAdd (timers, name, start);
    psFree (start);
    return (true);
}

psF64 psTimerClear (char *name)
{
    psF64 timerTime = 0.0;
    bool status;

    if (name == NULL)
        return false;
    timerTime = psTimerMark(name);
    status = psHashRemove (timers, name);
    strncpy(timerName, "\0", 256);
    if ( !status )
        return 0.0;
    else
        return timerTime;
}

// get current elapsed time on named timer
psF64 psTimerMark (char *name)
{

    psTime *start;
    psTime *mark;
    psF64   delta;

    if (timers == NULL)
        return (0);
    start = (psTime*)psHashLookup (timers, name);
    if (start == NULL)
        return (0);
    mark = psTimeGetNow (PS_TIME_UTC);
    delta = psTimeDelta (mark, start);
    psFree (mark);
    // psFree (start); -- XXX is psHashLookup not psMemCopying?

    return (delta);
}

psF64 psTimerStop(void)
{
    psF64 timerTime = 0.0;

    if ( strncmp(timerName, "\0", 256) )
        timerTime = psTimerMark(timerName);
    strncpy(timerName, "\0", 2);
    psTimerFree();
    return timerTime;
}

void psTimeInitialize(const char *timeConfig)
{
    p_psTimeInit(timeConfig);
}

void psTimeFinalize(void)
{
    p_psTimeFinalize();
}

