/** @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.  PSLib currently uses the UNIX timeval time system as the
 *  base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. TAI
 *  time varies over time due to the earth's rotation and the movement of the continental plates.
 *
 *  @author Ross Harman, MHPCC
 *
 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-15 19:10:50 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

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

#include "psTime.h"
#include "psError.h"
#include "psMemory.h"
#include "psAbort.h"
#include "psImage.h"
#include "psCoord.h"
#include "psString.h"

#include "psAstronomyErrors.h"

#ifndef SER7_FILE
#define SER7_FILE "../../data/ser7.dat"
#pragma warning SER7_FILE was not defined in the makefile.
#endif

#ifndef TAIUTC_FILE
#define TAIUTC_FILE "../../data/tai-utc.dat"
#pragma warning TAIUTC_FILE was not defined in the makefile.
#endif

/** 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

                /** Max size of a line read */
                #define LINESIZE 1024

                /** Preprocessor macro to check for null psTime struct */
                #define CHECK_NULL_TIME(TIME,RETURN) \
                if(TIME == NULL) { \
                    psError(__func__,"NULL value not allowed"); \
                    return RETURN; \
                }

                /** Preprocessor macro to allocate psTime struct */
                #define ALLOC_TIME(TIME) \
                TIME = (psTime*)psAlloc(sizeof(psTime)); \
if (TIME == NULL) { \
    psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__); \
} \
TIME->sec = 0; \
TIME->usec = 0; \
TIME->type = PS_TIME_TAI;

/** Read TAI/UTC data file.
 *
 *  Reads TAI/UTC data file.
 *
 *  @return  psImage*: Two-dimensional array with table data.
 */
static psImage* readTaiUtcFile(
    const char *fileName                      ///< Name of file to read.
);


/** Read SER7 data file.
 *
 *  Reads SER7 data file.
 *
 *  @return  psImage*: Two-dimensional array with table data.
 */
static psImage* readSer7File(
    const char *fileName                      ///< Name of file to read.
);


/** Read finals data file.
 *
 *  Reads finals data file.
 *
 *  @return  psImage*: Two-dimensional array with table data.
 */
/* not yet used
static psImage* readFinalsFile(
    char *fileName                      ///< Name of file to read.
);
*/

/** Read EOPC data file.
 *
 *  Reads EOPC data file.
 *
 *  @return  psImage*: Two-dimensional array with table data.
 */
/* not yet used
static psImage* readEopcFile(
    char *fileName                      ///< Name of file to read.
);
*/

/** Lookup and interpolate TAI-UTC data
 *
 *  Interpolates TAI-UTC data, given time in MJD format.
 *
 *  @return  double: Interpolated value.
 */
static double lookupTaiUtcTable(
    const psTime *time                        ///< time to lookup.
);

/** Lookup and interpolate polar or UT1-UTC data
 *
 *  Interpolates TAI-UT1 or polar coordinate data, given time in MJD format.
 *
 *  @return  double: Interpolated value.
 */
static double lookupSer7Table(
    const psTime *time,                        ///< time to lookup.
    int col                              ///< Column to lookup.
);

static psImage* readTaiUtcFile(const char *fileName)
{
    char line[LINESIZE];
    int j = 0;
    int maxLines = 100;
    psF64 *ptr = NULL;
    psImage *table = NULL;
    FILE *fd = NULL;


    fd=fopen(fileName, "r");
    if(fd==NULL) {
        psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName);
        return NULL;
    }

    // Table shouldn't be larger than 100 rows. Select columns are read based on indices.
    table  = psImageAlloc(4, maxLines, PS_TYPE_F64);

    while(fgets(line, LINESIZE, fd) != NULL) {
        if(j>=maxLines) {
            psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines);
            return NULL;
            psFree(table);
        }

        ptr = table->data.F64[j++];
        *ptr = atof(line+16);                  // MJD
        *(ptr+1) = atof(line+38);              // TAI-UTC
        *(ptr+2) = atof(line+59);              // Contant #1
        *(ptr+3) = atof(line+69);              // Contant #2
    }

    *(unsigned int *)&table->numRows = j;

    p_psMemSetPersistent(table,true);
    p_psMemSetPersistent(table->data.V,true);
    p_psMemSetPersistent(table->rawDataBuffer,true);

    return table;
}

static psImage* readSer7File(const char *fileName)
{
    bool beginRecord = false;
    char line[LINESIZE];
    int j = 0;
    int maxLines = 400;
    psImage *table = NULL;
    FILE *fd = NULL;


    fd=fopen(fileName, "r");
    if(fd==NULL) {
        psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName);
        return NULL ;
    }

    // Table shouldn't be larger than 400 rows. All columns in area of interest are read.
    table  = psImageAlloc(7, maxLines, PS_TYPE_F64);

    while(fgets(line, LINESIZE, fd) != NULL) {

        // Stop parsing after finding "These" in data file after read has began
        if(beginRecord) {
            if(strstr(line, "These")!=NULL) {
                beginRecord = false;
                break;
            }

            if(j>=maxLines) {
                psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines);
                return NULL;
                psFree(table);
            }
            psF64 *ptr = table->data.F64[j++];
            sscanf(line, "%lf %lf %lf %lf %lf %lf %lf", ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6);
        }

        // Start parsing after finding "x(arcsec)" in data file
        if(strstr(line, "x(arcsec)") != NULL) {
            beginRecord = true;
        }
    }
    *(unsigned int *)&table->numRows = j;

    p_psMemSetPersistent(table,true);
    p_psMemSetPersistent(table->data.V,true);
    p_psMemSetPersistent(table->rawDataBuffer,true);

    return table;
}

/* Not yet used...
static psImage* readEopcFile(char *fileName)
{
    char line[LINESIZE];
    int i = 0;
    int j = 0;
    int maxLines = 2500;
    psImage *table = NULL;
    FILE *fd = NULL;
 
 
    fd=fopen(fileName, "r");
    if(fd==NULL) {
        psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName);
        return NULL ;
    }
 
    // Table shouldn't be larger than 2500 rows. All columns are read.
    table  = psImageAlloc(11, maxLines, PS_TYPE_F64);
    for(i=0; i<3; i++) {
        if((fgets(line, LINESIZE, fd) != NULL)) {
            psError(__func__, " : Line %d - Error reading %s", __LINE__, fileName);
        }
    }
 
    while(fgets(line, LINESIZE, fd) != NULL) {
        if(j>=maxLines) {
            psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines);
            return NULL;
            psFree(table);
        }
        psF64 *ptr = table->data.F64[j++];
        sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
               ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6, ptr+7, ptr+8, ptr+9, ptr+10);
    }
 
    *(unsigned int *)&table->numRows = j;
 
    return table;
}
 
static psImage* readFinalsFile(char *fileName)
{
    char line[LINESIZE];
    int j = 0;
    int maxLines = 15000;
    psImage *table = NULL;
    FILE *fd = NULL;
 
 
    fd=fopen(fileName, "r");
    if(fd==NULL) {
        psError(__func__, " : Line %d - Couldn't open %s", __LINE__, fileName);
        return NULL ;
    }
 
    // Table shouldn't be larger than 15000 rows. Select columns are read.
    table  = psImageAlloc(6, maxLines, PS_TYPE_F64);
    while(fgets(line, LINESIZE, fd) != NULL) {
        if(j>=maxLines) {
            psError(__func__, " : Line %d - Too many rows in file. Max size: %d", __LINE__, maxLines);
            return NULL;
            psFree(table);
        }
 
        psF64 *ptr = table->data.F64[j++];
        ptr = table->data.F64[j++];
 
        // Sometimes there are blank entries in this table, so check for a space to determine if there is no entry
        if(*(line+7) == ' ') {                   // MJD
            *ptr = 0.0;
        } else {
            *ptr = atof(line+7);
        }
 
        if(*(line+18) == ' ') {                  // Polar motion X
            *(ptr+1) = 0.0;
        } else {
            *(ptr+1) = atof(line+18);
        }
 
        if(*(line+27) == ' ') {                  // Polar motion X error
            *(ptr+2) = 0.0;
        } else {
            *(ptr+2) = atof(line+27);
        }
 
        if(*(line+37) == ' ') {                  // Polar motion Y
            *(ptr+3) = 0.0;
        } else {
            *(ptr+3) = atof(line+37);
        }
 
        if(*(line+46) == ' ') {                  // Polar motion Y error
            *(ptr+4) = 0.0;
        } else {
            *(ptr+4) = atof(line+46);
        }
 
        if(*(line+58) == ' ') {                  // UT1-UTC
            *(ptr+5) = 0.0;
        } else {
            *(ptr+5) = atof(line+58);
        }
    }
 
    *(unsigned int *)&table->numRows = j;
 
    return table;
}
*/

static double lookupTaiUtcTable(const psTime *time)
{
    int hiIdx = 0;
    int loIdx = 0;
    int numRows = 0;
    double jd = 0.0;
    double mjd = 0.0;
    double out = 0.0;
    double denom = 0.0;
    double const1 = 0.0;
    double const2 = 0.0;
    double const3 = 0.0;
    static psImage* table = NULL;

    if (table == NULL) {
        table = readTaiUtcFile(TAIUTC_FILE);
        if (table == NULL) {
            psAbort(__func__,"Failed to read tai-utc.dat file");
        }
    }

    // Number of rows is the number of rows of data read from the data file
    numRows = table->numRows;

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

    // Variable in should be in MJD format
    if(jd < table->data.F64[0][0]) {                                        // MJD time before start of table
        return 0.0;
    } else if(jd > table->data.F64[numRows-1][0]) {                         // MJD time after end of table
        return table->data.F64[numRows-1][1];
    } else {                                                                // All other times..interpolate
        while(jd > table->data.F64[hiIdx][0]) {
            hiIdx++;
            if(hiIdx >= numRows) {
                psError(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx);
                return 0.0;
            }
        }
        loIdx = hiIdx--;
        if(loIdx < 0) {
            psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx);
            return 0.0;
        }

        denom = table->data.F64[hiIdx][0]-table->data.F64[loIdx][0];
        if(fabs(denom) < FLT_EPSILON) {
            psError(__func__, " : Line %d - Divide by zero error during interpolation", __LINE__);
            return 0.0;
        } else {

            // Interpolate three constants required to calculate TAI-UTC
            const1 = table->data.F64[loIdx][1]+(table->data.F64[hiIdx][1]-table->data.F64[loIdx][1])
                     *(jd-table->data.F64[loIdx][0])/denom;

            const2 = table->data.F64[loIdx][2]+(table->data.F64[hiIdx][2]-table->data.F64[loIdx][2])
                     *(jd-table->data.F64[loIdx][0])/denom;

            const3 = table->data.F64[loIdx][3]+(table->data.F64[hiIdx][3]-table->data.F64[loIdx][3])
                     *(jd-table->data.F64[loIdx][0])/denom;
        }
    }

    out = const1 + (mjd - const2) * const3;

    return out;
}

static double lookupSer7Table(const psTime *time, int col)
{
    int hiIdx = 0;
    int loIdx = 0;
    int numRows = 0;
    double out = 0.0;
    double denom = 0.0;
    double mjdUtc = 0.0;
    static psImage* table = NULL;

    if (table == NULL) {
        table = readSer7File(SER7_FILE);
        if (table == NULL) {
            psAbort(__func__,"Failed to read the ser7.dat file.");
        }
    }

    // Number of rows is the number of rows of data read from the data file
    numRows = table->numRows;

    // Interpolation value in should be in MJD format
    mjdUtc = psTimeToMJD(time);
    if(mjdUtc < table->data.F64[0][3]) {                                    // MJD time before start of table
        psError(__func__, " : Line %d - MJD too small. Value: %d", __LINE__, mjdUtc);

    } else if(mjdUtc > table->data.F64[numRows-1][3]) {                     // MJD time after end of table
        psError(__func__, " : Line %d - MJD too big. Value: %d", __LINE__, mjdUtc);
    } else {                                                                // All other times..interpolate
        while(mjdUtc > table->data.F64[hiIdx][3]) {
            hiIdx++;
            if(hiIdx >= numRows) {
                psError(__func__, " : Line %d - High index too big. Value: %d", __LINE__, hiIdx);
                return 0.0;
            }
        }

        loIdx = hiIdx--;
        if(loIdx < 0) {
            psError(__func__, " : Line %d - Low index too small. Value: %d", __LINE__, loIdx);
        }

        denom = table->data.F64[hiIdx][3] - table->data.F64[loIdx][3];
        if(fabs(denom) < FLT_EPSILON) {
            psError(__func__, " : Line %d - Dividie by zero error during interpolation", __LINE__);
        } else {
            out = table->data.F64[loIdx][col]+(table->data.F64[hiIdx][col]-table->data.F64[loIdx][col])
                  *(mjdUtc-table->data.F64[loIdx][3])/denom;
        }
    }

    return out;
}

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


    // Allocate psTime struct
    ALLOC_TIME(time);

    if (gettimeofday(&now, (struct timezone *)0) == -1) {
        psError(__func__, " : Line %d - Failed to get time", __LINE__);
        return NULL;
    }

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

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


    if(time->type == type) {
        psError(__func__, " : Line %d - Time is already in the correct psTime format. Time type: %d", __LINE__, type);
        return NULL;
    }

    delta = psTimeGetTAIDelta(time);

    if(type == PS_TIME_UTC) {
        time->sec = time->sec - delta;    // User wants UTC time. Remove leapseconds.
    } else if(type == PS_TIME_TAI) {
        time->sec = time->sec + delta;    // User wants TAI time. Add leapseconds.
    } else {
        psError(__func__, " : Line %d - Invalid psTimeType: %d", __LINE__, 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 *tdtDelta    = NULL;
    psTime *tdtTime     = NULL;
    psTime *taiTime     = NULL;
    psTime *utcTime     = NULL;
    psTime *ut1Time     = NULL;


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

    // Convert Universal Time (UTC) to  UT1
    ALLOC_TIME(ut1UtcDelta);
    ut1UtcDelta->type = PS_TIME_UTC;

    ut1UtcDbl = psTimeGetUT1Delta(utcTime)*1e6;
    ut1UtcDelta->usec = (psU32)fabs(ut1UtcDbl);
    if(ut1UtcDbl<0) {
        ut1Time = psTimeSubtract(utcTime, ut1UtcDelta);
    } else {
        ut1Time = psTimeAdd(utcTime, ut1UtcDelta);
    }


    // 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)
    ALLOC_TIME(tdtDelta);
    tdtDelta->sec = 32;
    tdtDelta->usec = 184000;
    tdtTime = psTimeAdd(taiTime, tdtDelta);

    // 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(tdtDelta);
    psFree(utcTime);
    psFree(taiTime);

    return lmstRad;
}

double psTimeGetUT1Delta(const psTime *time)
{
    return lookupSer7Table(time, 6);
}

struct psSphere* psTimeGetPoleCoords(const psTime* time)
{

    double x = lookupSer7Table((psTime*)time, 4);
    double y = lookupSer7Table((psTime*)time, 5);

    struct psSphere* output = psAlloc(sizeof(psSphere));
    output->r = x * M_PI / 648000.0; // x in arcsec converted to radians, i.e., x/60/60*M_PI/180
    output->d = y * M_PI / 648000.0;

    return output;
}

double psTimeGetTAIDelta(const psTime *time)
{
    return lookupTaiUtcTable(time);
}

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

    // NULL error checks
    CHECK_NULL_TIME(time1,0);
    CHECK_NULL_TIME(time2,0);

    diff = abs(psTimeGetTAIDelta((psTime*)time1)-psTimeGetTAIDelta((psTime*)time2));

    return diff;
}

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


    // NULL error check
    CHECK_NULL_TIME(time,0.0);

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


    // NULL error check
    CHECK_NULL_TIME(time,0.0);

    // 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)
{
    int ms = 0;
    char *timeString = NULL;
    char *tempString = NULL;
    struct tm *tmTime = NULL;
    time_t sec;


    // NULL error check
    CHECK_NULL_TIME(time,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(__func__, " : Line %d - Failed strftime conversion", __LINE__);
    }

    if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s.%3.3dZ", tempString, ms) < 0) {
        psError(__func__, " : Line %d - Failed snprintf conversion", __LINE__);
    }
    psFree(tempString);

    return timeString;
}

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


    // NULL error check
    CHECK_NULL_TIME(time,timevalTime);

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

    return timevalTime;
}

struct tm* psTimeToTM(const psTime *time)
{
    long cent = 0;
    long year = 0;
    long month = 0;
    long day = 0;
    long hour = 0;
    long minute = 0;
    long seconds = 0;
    long temp = 0;
    struct tm* tmTime = 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));
    if (tmTime == NULL)
    {
        psAbort(__func__, " : Line %d - Failed to allocate memory", __LINE__);
    }

    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
    ALLOC_TIME(outTime);

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

    return outTime;
}

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


    // Allocate psTime struct
    ALLOC_TIME(outTime);

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

    return outTime;
}

psTime* psTimeFromISOTime(const char *time)
{
    char tempString[MAX_TIME_STRING_LENGTH];
    int month;
    int day;
    int year;
    int hour;
    int minute;
    int second;
    int millisecond;
    struct tm tmTime;
    psTime *outTime = NULL;


    // Preserve input string by creating duplicate
    strncpy(tempString, time, MAX_TIME_STRING_LENGTH);

    // Convert YYYY-MM-DDThh:mm:ss.sss in string form to tm time
    year = atoi(strtok(tempString, "-"));
    if (year < 0) {
        psError(__func__, "Years less than 0 not allowed. Value: %d", year);
        return outTime;
    }

    month = atoi(strtok(NULL, "-"));
    if (month < 1 || month > 12) {
        psError(__func__, "Month must have a value from 1 to 12. Value: %d", month);
        return outTime;
    }

    day = atoi(strtok(NULL, "T"));
    if (day < 1 || day > 31) {
        psError(__func__, "Day must have a value from 1 to 31. Value: %d", day);
        return outTime;
    }

    hour = atoi(strtok(NULL, ":"));
    if (hour < 0 || hour > 23) {
        psError(__func__, "Hour must have a value from 0 to 23. Value: %d", hour);
        return outTime;
    }

    minute = atoi(strtok(NULL, ":"));
    if (minute < 0 || minute > 59) {
        psError(__func__, "Minute must have a value from 0 to 59. Value: %d", minute);
        return outTime;
    }

    second = atoi(strtok(NULL, "."));
    if (second < 0 || second > 59) {
        psError(__func__, "Second must have a value from 0 to 59. Value: %d", second);
        return outTime;
    }

    millisecond = atoi(strtok(NULL, "Z"));
    if (millisecond < 0 || millisecond > 1000) {
        psError(__func__, "Millisecond must have a value from 0 to 999. Value: %d", millisecond);
        return outTime;
    }

    tmTime.tm_year = year - 1900;
    tmTime.tm_mon = month - 1;
    tmTime.tm_mday = day;
    tmTime.tm_hour = hour;
    tmTime.tm_min = minute;
    tmTime.tm_sec = second;
    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;


    // NULL error check
    CHECK_NULL_TIME(time,NULL);

    // Allocate psTime struct
    ALLOC_TIME(outTime)

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

    return outTime;
}

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

    // NULL error check
    CHECK_NULL_TIME(time,NULL);

    // Allocate psTime struct
    ALLOC_TIME(outTime);

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

    return outTime;
}

psTime* psTimeAdd(const psTime *tai1, const psTime *tai2)
{
    long deltaUsec = 0;
    psTime *outTime = NULL;


    if(tai1->type != tai2->type) {
        psError(__func__, " : Line %d - Argument types inconsistent for time addition", __LINE__);
        return NULL;
    }

    ALLOC_TIME(outTime);
    outTime->type = tai1->type;
    outTime->sec = tai1->sec + tai2->sec;
    deltaUsec = tai1->usec + tai2->usec;

    // Adjust time in case of microsecond overflow after subtraction
    if(deltaUsec > 1e6) {
        outTime->sec++;
        outTime->usec = 2e6 - deltaUsec;
    } else {
        outTime->usec = deltaUsec;
    }

    return outTime;
}

psTime* psTimeSubtract(const psTime *tai1, const psTime *tai2)
{
    long deltaUsec = 0;
    psTime *outTime = NULL;


    if(tai1->type != tai2->type) {
        psError(__func__, " : Line %d - Argument types inconsistent for time subtraction", __LINE__);
        return NULL;
    }

    ALLOC_TIME(outTime);
    outTime->type = tai1->type;
    outTime->sec = tai1->sec - tai2->sec;
    deltaUsec = (long)tai1->usec - (long)tai2->usec;

    // Adjust time in case of microsecond underflow after subtraction
    if(deltaUsec < 0) {
        outTime->sec--;
        outTime->usec = 1e6+deltaUsec;
    } else {
        outTime->usec = deltaUsec;
    }

    return outTime;
}

psTime* psTimeDelta(const psTime *tai1, const psTime *tai2)
{
    long deltaUsec = 0;
    psTime *outTime = NULL;


    if(tai1->type != tai2->type) {
        psError(__func__, " : Line %d - Argument types inconsistent for time subtraction", __LINE__);
        return NULL;
    }

    ALLOC_TIME(outTime);
    outTime->type = tai1->type;
    outTime->sec = tai1->sec - tai2->sec;
    deltaUsec = tai1->usec - tai2->usec;

    // Adjust time in case of microsecond underflow after subtraction
    if(deltaUsec < 0) {
        outTime->sec--;
        outTime->usec = 1e6+deltaUsec;
    } else {
        outTime->usec = deltaUsec;
    }

    outTime->sec = abs(outTime->sec);

    return outTime;
}
