IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 8, 2004, 2:36:13 PM (22 years ago)
Author:
harman
Message:

Time config files updates

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astronomy/psTime.c

    r2273 r2301  
    55 *
    66 *  A collection of functions are required by psLib to manipulate time data. These functions primarily consist
    7  *  of conversions between specific time formats.  PSLib currently uses the UNIX timeval time system as the
    8  *  base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. TAI
    9  *  time varies over time due to the earth's rotation and the movement of the continental plates.
     7 *  of conversions between specific time formats.  They use the UNIX timeval time system as the
     8 *  base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated.
    109 *
    1110 *  @author Ross Harman, MHPCC
    1211 *
    13  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2004-11-04 01:04:57 $
     12 *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-11-09 00:36:13 $
    1514 *
    1615 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2120#include <string.h>
    2221#include <math.h>
     22#include <ctype.h>
    2323
    2424#include "psTime.h"
    … …  
    3131#include "psMetadata.h"
    3232#include "psMetadataIO.h"
     33#include "psLookupTable.h"
    3334#include "psConstants.h"
    3435#include "psAstronomyErrors.h"
    3536
    36 #ifndef SER7_FILE
    37 #define SER7_FILE "../../data/ser7.dat"
    38 #pragma warning SER7_FILE was not defined in the makefile.
     37#ifndef TIME_CONFIG_FILE
     38#define TIME_CONFIG_FILE "../../config/time.config"
     39#pragma warning TIME_CONFIG_FILE was not defined in the makefile.
    3940#endif
    40 
    41 #ifndef TAIUTC_FILE
    42 #define TAIUTC_FILE "../../data/tai-utc.dat"
    43 #pragma warning TAIUTC_FILE was not defined in the makefile.
    44 #endif
    45 
    4641
    4742/** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */
    … …  
    7267                #define USEC_PER_DAY 86400000000.0
    7368
    74                 /** Max size of a line read */
    75                 #define LINESIZE 1024
    76 
    77 
    7869                /** Time metadata read from config file */
    7970                static psMetadata *timeMetadata = NULL;
    8071
    81 
    82 /** Read TAI/UTC data file.
    83  *
    84  *  Reads TAI/UTC data file.
    85  *
    86  *  @return  psImage*: Two-dimensional array with table data.
    87  */
    88 static psImage* readTaiUtcFile(
    89     const char *fileName                      ///< Name of file to read.
    90 );
    91 
    92 
    93 /** Read SER7 data file.
    94  *
    95  *  Reads SER7 data file.
    96  *
    97  *  @return  psImage*: Two-dimensional array with table data.
    98  */
    99 static psImage* readSer7File(
    100     const char *fileName                      ///< Name of file to read.
    101 );
    102 
    103 
    104 /** Read finals data file.
    105  *
    106  *  Reads finals data file.
    107  *
    108  *  @return  psImage*: Two-dimensional array with table data.
    109  */
    110 /* not yet used
    111 static psImage* readFinalsFile(
    112     char *fileName                      ///< Name of file to read.
    113 );
    114 */
    115 
    116 /** Read EOPC data file.
    117  *
    118  *  Reads EOPC data file.
    119  *
    120  *  @return  psImage*: Two-dimensional array with table data.
    121  */
    122 /* not yet used
    123 static psImage* readEopcFile(
    124     char *fileName                      ///< Name of file to read.
    125 );
    126 */
    127 
    128 /** Lookup and interpolate TAI-UTC data
    129  *
    130  *  Interpolates TAI-UTC data, given time in MJD format.
    131  *
    132  *  @return  double: Interpolated value.
    133  */
    134 static double lookupTaiUtcTable(
    135     const psTime *time                  ///< time to lookup.
    136 );
    137 
    138 /** Lookup and interpolate polar or UT1-UTC data
    139  *
    140  *  Interpolates TAI-UT1 or polar coordinate data, given time in MJD format.
    141  *
    142  *  @return  double: Interpolated value.
    143  */
    144 static double lookupSer7Table(
    145     const psTime *time,                 ///< time to lookup.
    146     psS32 col                             ///< Column to lookup.
    147 );
    148 
    149 static psImage* readTaiUtcFile(const char *fileName)
    150 {
    151     char line[LINESIZE];
    152     psS32 j = 0;
    153     psS32 maxLines = 1000;
    154     psF64 *ptr = NULL;
    155     psImage *table = NULL;
    156     FILE *fd = NULL;
    157 
    158 
    159     fd=fopen(fileName, "r");
    160     if(fd==NULL) {
    161         psError(PS_ERR_IO,true,
    162                 PS_ERRORTEXT_psTime_FILE_NOT_FOUND,
    163                 fileName,"Tai-Utc");
    164         return NULL;
    165     }
    166 
    167     // Table shouldn't be larger than 100 rows. Select columns are read based on indices.
    168     table  = psImageAlloc(4, maxLines, PS_TYPE_F64);
    169 
    170     while(fgets(line, LINESIZE, fd) != NULL) {
    171         if(j>=maxLines) {
    172             psError(PS_ERR_IO,true,
    173                     PS_ERRORTEXT_psTime_FILE_TOO_MANY_ROWS,
    174                     fileName,maxLines);
    175             return NULL;
    176             psFree(table);
    177         }
    178 
    179         ptr = table->data.F64[j++];
    180         *ptr = atof(line+16);                  // MJD
    181         *(ptr+1) = atof(line+38);              // TAI-UTC
    182         *(ptr+2) = atof(line+59);              // Contant #1
    183         *(ptr+3) = atof(line+69);              // Contant #2
    184     }
    185 
    186     *(psU32 *)&table->numRows = j;
    187 
    188     p_psMemSetPersistent(table,true);
    189     p_psMemSetPersistent(table->data.V,true);
    190     p_psMemSetPersistent(table->rawDataBuffer,true);
    191 
    192     return table;
    193 }
    194 
    195 static psImage* readSer7File(const char *fileName)
    196 {
    197     psBool beginRecord = false;
    198     char line[LINESIZE];
    199     psS32 j = 0;
    200     psS32 maxLines = 400;
    201     psImage *table = NULL;
    202     FILE *fd = NULL;
    203 
    204 
    205     fd=fopen(fileName, "r");
    206     if(fd==NULL) {
    207         psError(PS_ERR_IO,true,
    208                 PS_ERRORTEXT_psTime_FILE_NOT_FOUND,
    209                 fileName,"Ser7");
    210         return NULL ;
    211     }
    212 
    213     // Table shouldn't be larger than 400 rows. All columns in area of interest are read.
    214     table  = psImageAlloc(7, maxLines, PS_TYPE_F64);
    215 
    216     while(fgets(line, LINESIZE, fd) != NULL) {
    217 
    218         // Stop parsing after finding "These" in data file after read has began
    219         if(beginRecord) {
    220             if(strstr(line, "These")!=NULL) {
    221                 beginRecord = false;
    222                 break;
    223             }
    224 
    225             if(j>=maxLines) {
    226                 psError(PS_ERR_IO,true,
    227                         PS_ERRORTEXT_psTime_FILE_TOO_MANY_ROWS,
    228                         fileName,maxLines);
    229                 return NULL;
    230                 psFree(table);
    231             }
    232             psF64 *ptr = table->data.F64[j++];
    233             sscanf(line, "%lf %lf %lf %lf %lf %lf %lf", ptr, ptr+1, ptr+2, ptr+3, ptr+4, ptr+5, ptr+6);
    234         }
    235 
    236         // Start parsing after finding "x(arcsec)" in data file
    237         if(strstr(line, "x(arcsec)") != NULL) {
    238             beginRecord = true;
    239         }
    240     }
    241     *(psU32 *)&table->numRows = j;
    242 
    243     p_psMemSetPersistent(table,true);
    244     p_psMemSetPersistent(table->data.V,true);
    245     p_psMemSetPersistent(table->rawDataBuffer,true);
    246 
    247     return table;
    248 }
    249 
    250 
    251 static double lookupTaiUtcTable(const psTime *time)
    252 {
    253     psS32 hiIdx = 0;
    254     psS32 loIdx = 0;
    255     psS32 numRows = 0;
    256     double jd = 0.0;
    257     double mjd = 0.0;
    258     double out = 0.0;
    259     double denom = 0.0;
    260     double const1 = 0.0;
    261     double const2 = 0.0;
    262     double const3 = 0.0;
    263     static psImage* table = NULL;
    264 
    265     if (table == NULL) {
    266         table = readTaiUtcFile(TAIUTC_FILE);
    267         if (table == NULL) {
    268             psAbort(__func__,"Failed to read tai-utc.dat file");
    269         }
    270     }
    271 
    272     // Number of rows is the number of rows of data read from the data file
    273     numRows = table->numRows;
    274 
    275     // Determine Julian and modified Julian dates used in table lookup and time delta calculation
    276     jd = psTimeToJD(time);
    277     mjd = psTimeToMJD(time);
    278 
    279     // Variable in should be in MJD format
    280     if(jd < table->data.F64[0][0]) {                                        // MJD time before start of table
    281         return 0.0;
    282     } else if(jd > table->data.F64[numRows-1][0]) {                         // MJD time after end of table
    283         return table->data.F64[numRows-1][1];
    284     } else {                                                                // All other times..interpolate
    285         while(jd > table->data.F64[hiIdx][0]) {
    286             hiIdx++;
    287             if(hiIdx >= numRows) {
    288                 psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    289                         PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE,
    290                         jd, "TAI-UTC");
    291                 return 0.0;
    292             }
    293         }
    294         loIdx = hiIdx--;
    295         if(loIdx < 0) {
    296             psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    297                     PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE,
    298                     jd, "TAI-UTC");
     72/** Static function prototypes */
     73static char *cleanString(char *inString, int sLen);
     74static char* getToken(char **inString, char *delimiter, psParseErrorType *status);
     75psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables);
     76
     77/** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null
     78 *  terminated copy of the original input string. */
     79static char *cleanString(char *inString, int sLen)
     80{
     81    char *ptrB = NULL;
     82    char *ptrE = NULL;
     83    char *cleaned = NULL;
     84
     85
     86    ptrB = inString;
     87
     88    /* Skip over leading # or whitespace */
     89    while (isspace(*ptrB) || *ptrB=='#') {
     90        ptrB++;
     91    }
     92
     93    /* Skip over trailing whitespace, null terminators, and # characters */
     94    ptrE = inString + sLen;
     95    while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') {
     96        ptrE--;
     97    }
     98
     99    // Length, sLen, does not include '\0'
     100    sLen = ptrE - ptrB + 1;
     101
     102    // Adds '\0' to end of string and +1 to sLen
     103    cleaned = psStringNCopy(ptrB, sLen);
     104
     105    return cleaned;
     106}
     107
     108/** Returns cleaned token based on delimiter, but not including delimiter. Also changes the pointer location
     109 * the beginning of the string. Tokens are newly allocated null terminated strings. */
     110static char* getToken(char **inString, char *delimiter, psParseErrorType *status)
     111{
     112    char *cleanToken = NULL;
     113    int sLen = 0;
     114
     115
     116    // Skip over leading whitespace
     117    while(isspace(**inString)) {
     118        (*inString)++;
     119    }
     120
     121    // Length of token, not including delimiter
     122    sLen = strcspn(*inString, delimiter);
     123    if(sLen) {
     124
     125        // Create new, cleaned, and null terminated token
     126        cleanToken = cleanString(*inString, sLen);
     127
     128        // Move to end of token
     129        (*inString) += sLen;
     130    } else if(**inString!='\0' && sLen==0) {
     131        *status = PS_PARSE_ERROR_GENERAL;
     132    }
     133
     134    return cleanToken;
     135}
     136
     137/** Searches time tables in priority order and performs interpolation if input index value is within a table.
     138 * If the index value is out of range, the status is set accordingly. */
     139psF64 searchTables(psF64 index, psU64 column, psLookupStatusType *status, char *metadataTableNames[], psU32 nTables)
     140{
     141    psU32 i = 0;
     142    char *tableName = NULL;
     143    psF64 result = 0.0;
     144    psLookupTable *table = NULL;
     145    psMetadataItem *tableMetadataItem = NULL;
     146
     147
     148
     149    // Check time metadata. Function call reports errors.
     150    if(timeMetadata == NULL) {
     151        if(!psTimeInit(TIME_CONFIG_FILE))
    299152            return 0.0;
    300         }
    301 
    302         // since rows of table can not have same jd, this check is probably not needed.
    303         denom = table->data.F64[hiIdx][0]-table->data.F64[loIdx][0];
    304         if(fabs(denom) < FLT_EPSILON) {
    305             psError(PS_ERR_UNKNOWN,true,
    306                     PS_ERRORTEXT_psTime_TABLE_DUPLICATE_ROWS
    307                     "Ser7");
     153    }
     154
     155    // Search each table in priority order: ser7, eopc,finals
     156    for(i=0; i<nTables; i++) {
     157
     158        // Get table from metadata
     159        tableName = metadataTableNames[i];
     160        tableMetadataItem = psMetadataLookup(timeMetadata, tableName);
     161        if(tableMetadataItem == NULL) {
     162            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     163                    tableName);
    308164            return 0.0;
    309         } else {
    310             // Interpolate three constants required to calculate TAI-UTC
    311             const1 = table->data.F64[loIdx][1]+(table->data.F64[hiIdx][1]-table->data.F64[loIdx][1])
    312                      *(jd-table->data.F64[loIdx][0])/denom;
    313 
    314             const2 = table->data.F64[loIdx][2]+(table->data.F64[hiIdx][2]-table->data.F64[loIdx][2])
    315                      *(jd-table->data.F64[loIdx][0])/denom;
    316 
    317             const3 = table->data.F64[loIdx][3]+(table->data.F64[hiIdx][3]-table->data.F64[loIdx][3])
    318                      *(jd-table->data.F64[loIdx][0])/denom;
    319         }
    320     }
    321 
    322     out = const1 + (mjd - const2) * const3;
    323 
    324     return out;
    325 }
    326 
    327 static double lookupSer7Table(const psTime *time, psS32 col)
    328 {
    329     psS32 hiIdx = 0;
    330     psS32 loIdx = 0;
    331     psS32 numRows = 0;
    332     double out = 0.0;
    333     double denom = 0.0;
    334     double mjdUtc = 0.0;
    335     static psImage* table = NULL;
    336 
    337     if (table == NULL) {
    338         table = readSer7File(SER7_FILE);
    339         if (table == NULL) {
    340             psAbort(__func__,"Failed to read the ser7.dat file.");
    341         }
    342     }
    343 
    344     // Number of rows is the number of rows of data read from the data file
    345     numRows = table->numRows;
    346 
    347     // Interpolation value in should be in MJD format
    348     mjdUtc = psTimeToMJD(time);
    349     if(mjdUtc < table->data.F64[0][3]) {                                    // MJD time before start of table
    350         psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    351                 PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE,
    352                 mjdUtc, "Ser7");
    353 
    354     } else if(mjdUtc > table->data.F64[numRows-1][3]) {                     // MJD time after end of table
    355         psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    356                 PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE,
    357                 mjdUtc, "Ser7");
    358     } else {                                                                // All other times..interpolate
    359         while(mjdUtc > table->data.F64[hiIdx][3]) {
    360             hiIdx++;
    361             if(hiIdx >= numRows) {
    362                 psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    363                         PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE,
    364                         mjdUtc, "Ser7");
    365                 return 0.0;
    366             }
    367         }
    368 
    369         loIdx = hiIdx--;
    370         if(loIdx < 0) {
    371             psError(PS_ERR_BAD_PARAMETER_VALUE,true,
    372                     PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE,
    373                     mjdUtc, "Ser7");
    374         }
    375 
    376         denom = table->data.F64[hiIdx][3] - table->data.F64[loIdx][3];
    377         if(fabs(denom) < FLT_EPSILON) {
    378             psError(PS_ERR_UNKNOWN,true,
    379                     PS_ERRORTEXT_psTime_TABLE_DUPLICATE_ROWS,
    380                     "Ser7");
    381         } else {
    382             out = table->data.F64[loIdx][col]+(table->data.F64[hiIdx][col]-table->data.F64[loIdx][col])
    383                   *(mjdUtc-table->data.F64[loIdx][3])/denom;
    384         }
    385     }
    386 
    387     return out;
    388 }
    389 
    390 void psTimeInit(char *fileName)
    391 {
     165        }
     166        table = (psLookupTable*)tableMetadataItem->data.V;
     167        PS_PTR_CHECK_NULL(table,0.0);
     168
     169        // Attempt to interpolate table
     170        *status = PS_LOOKUP_SUCCESS;
     171        result = psLookupTableInterpolate(table, index, 2, status);
     172        if(*status == PS_LOOKUP_SUCCESS) {
     173            return result;
     174        } else if(*status == PS_LOOKUP_ERROR) {
     175            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED_NAME,
     176                    tableName);
     177            return 0.0;
     178        }
     179    }
     180
     181    return result;
     182}
     183
     184bool psTimeInit(char *fileName)
     185{
     186    bool foundTable = false;
    392187    char *tableDir = NULL;
    393188    char *tableNames = NULL;
    394     psU8 numTables = 0;
     189    char *namesPtr = NULL;
     190    char *metadataNamesPtr = NULL;
     191    char *tableName = NULL;
     192    char *fullTableName = NULL;
     193    psS32 i = 0;
     194    psS32 j = 0;
     195    psS32 numTables = 0;
    395196    psVector *tablesFrom = NULL;
    396197    psVector *tablesTo = NULL;
    397     psF64 beforeXp = 0.0;
    398     psF64 beforeYp = 0.0;
    399     psF64 beforeDut = 0.0;
    400198    psMetadataItem *metadataItem = NULL;
     199    psLookupTable *table = NULL;
     200    psLookupStatusType status = PS_LOOKUP_SUCCESS;
     201    char *metadataTableNames[4] = {"ser7", "eopc",  "finals", "tai"};
    401202
    402203
    403204    // Read time config file
    404     psMetadataParseConfig(&timeMetadata, fileName, true);
     205    if(psMetadataParseConfig(&timeMetadata, fileName, true) < 0) {
     206        return false;
     207    }
    405208
    406209    // Get path to time data files
    407210    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.dir");
    408     tableDir = (char*)metadataItem->data.V;
     211    if(metadataItem == NULL) {
     212        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     213                "psLib.time.tables.dir");
     214        return false;
     215    }
     216    tableDir = psStringCopy(metadataItem->data.V);
    409217
    410218    // Get number of tables
    411219    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.n");
    412     numTables = (psU8)metadataItem->data.S32;
    413 
    414     // Get table names
     220    if(metadataItem == NULL) {
     221        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     222                "psLib.time.tables.n");
     223        return false;
     224    }
     225    numTables = (psS32)metadataItem->data.S32;
     226
     227    // Table file names
    415228    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.files");
    416     tableNames = (char*)metadataItem->data.V;
    417 
    418     // Get ranges of tables
     229    if(metadataItem == NULL) {
     230        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     231                "psLib.time.tables.files");
     232        return false;
     233    }
     234    tableNames = psStringCopy(metadataItem->data.V);
     235
     236    // Get lower range of tables
    419237    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.from");
    420     tablesFrom = (psVector*)metadataItem->data.V;
     238    if(metadataItem == NULL) {
     239        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     240                "psLib.time.tables.from");
     241        return false;
     242    }
     243    tablesFrom = psVectorCopy(tablesFrom, metadataItem->data.V, PS_TYPE_F64);
     244    if(tablesFrom->n != numTables) {
     245        return false;
     246    }
     247
     248    // Get upper range of tables
    421249    metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.to");
    422     tablesTo = (psVector*)metadataItem->data.V;
    423 
    424     // Out of range values
    425     metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp");
    426     beforeXp =(psF64)metadataItem->data.F64;
    427     metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp");
    428     beforeYp =(psF64)metadataItem->data.F64;
    429     metadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut");
    430     beforeDut = (psF64)metadataItem->data.F64;
     250    if(metadataItem == NULL) {
     251        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     252                "psLib.time.tables.to");
     253        return false;
     254    }
     255    tablesTo = psVectorCopy(tablesTo, metadataItem->data.V, PS_TYPE_F64);
     256    if(tablesTo->n != numTables) {
     257        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesTo->n, numTables);
     258        return false;
     259    }
     260
     261    // Read time tables
     262    namesPtr = tableNames;
     263    while((tableName=getToken(&namesPtr, " ", &status)) != NULL) {
     264        fullTableName = (char*)psAlloc(strlen(tableDir)+strlen(tableName)+1+1); // Add one for last dir "/"
     265        fullTableName[0]='\0'; // Old strings may come back from psAlloc, so set init pos to EOL
     266        strcat(fullTableName, tableDir);
     267        strcat(fullTableName, "/");
     268        strcat(fullTableName, tableName);
     269
     270        if(i < numTables) {
     271            table = psLookupTableAlloc(fullTableName, tablesFrom->data.F64[i], tablesTo->data.F64[i]);
     272            psLookupTableRead(table);
     273        } else {
     274            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables);
     275        }
     276
     277        // Place tables into metadata slightly altered names as keys to create consistent naming conventions
     278        foundTable = false;
     279        for(j=0; j<numTables; j++) {
     280            metadataNamesPtr = strstr(tableName, metadataTableNames[j]);
     281            if(metadataNamesPtr != NULL) {
     282                psMetadataAdd(timeMetadata, PS_LIST_TAIL, strcat(metadataTableNames[j], "Table"),
     283                              PS_META_LOOKUPTABLE, NULL, table);
     284                foundTable = true;
     285            } else if(foundTable==false && j==numTables-1) {
     286                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, j, numTables);
     287            }
     288        }
     289
     290        psFree(fullTableName);
     291        psFree(tableName);
     292        i++;
     293    }
     294
     295    if(numTables > i+1) {
     296        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables);
     297    }
     298
     299    psFree(tableDir);
     300    psFree(tableNames);
     301    psFree(tablesFrom);
     302    psFree(tablesTo);
     303
     304    return true;
     305}
     306
     307void psTimeFinalize(void)
     308{
     309    psFree(timeMetadata);
    431310}
    432311
    … …  
    552431    // Convert Universal Time (UTC) to  UT1
    553432    ut1UtcDelta = psTimeAlloc(PS_TIME_UTC);
    554 
    555     ut1UtcDbl = psTimeGetUT1Delta(utcTime)*1e6;
     433    ut1UtcDbl = psTimeGetUT1Delta(taiTime)*1e6;
    556434    ut1UtcDelta->usec = (psU32)fabs(ut1UtcDbl);
    557435    if(ut1UtcDbl<0) {
    … …  
    561439    }
    562440
    563 
    564441    // Calculate UT1 as Julian Centuries since J2000.0
    565442    jdUt1Days = psTimeToJD(ut1Time);
    … …  
    580457    fracDays = fmod(mjdUt1Days, 1.0);
    581458
    582     // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications
     459    // Calculate Greenwich Mean Sidereal Time (GMST) in radians. Equation set up to minimize multiplications.
    583460    gmstRad = fracDays*TWOPI+(const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R;
    584461
    … …  
    602479double psTimeGetUT1Delta(const psTime *time)
    603480{
     481    psU32 nTables = 3;
     482    psF64 mjd = 0.0;
     483    psF64 result = 0.0;
     484    psF64 dut2ut1 = 0.0;
     485    psF64 t = 0.0;
     486    psVector *dut = NULL;
     487    psMetadataItem *tableMetadataItem = NULL;
     488    psLookupStatusType status = PS_LOOKUP_SUCCESS;
     489    char *metadataTableNames[3] = {"ser7", "eopc",  "finals"};
     490
    604491
    605492    // Error checks
    … …  
    607494    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);
    608495
    609     return lookupSer7Table(time, 6);
     496    if(time->type != PS_TIME_TAI) {
     497        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type);
     498        return 0.0;
     499    }
     500
     501    // Attempt to find value through table lookup and interpolation
     502    mjd = psTimeToMJD(time);
     503    result = searchTables(mjd, 0, &status, metadataTableNames, nTables);
     504
     505    // Value could not be found through table lookup and interpolation
     506    if(status == PS_LOOKUP_PAST_TOP) {
     507
     508        // Date too early for tables. Get default time delta value from metadata, and issue warning.
     509        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "UT1-UTC");
     510
     511        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut");
     512        if(tableMetadataItem == NULL) {
     513            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.dut");
     514            return 0.0;
     515        }
     516        result = tableMetadataItem->data.F64;
     517
     518    } else if(status == PS_LOOKUP_PAST_BOTTOM) {
     519        /* Date too late for tables. Issue warning and use following formulae for predicting
     520           ahead of the most recent available table entry.
     521             ut1-utc = [0] + [1]*(MJD - [2]) - (ut2-ut1)
     522             [0, 1, 2] = @psLib.time.predict.dut
     523             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)
     524             t = 2000.0 + (MJD - 51544.03)/365.2422
     525        */
     526        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC");
     527
     528        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.dut");
     529        if(tableMetadataItem == NULL) {
     530            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.dut");
     531            return 0.0;
     532        }
     533        dut = (psVector*)tableMetadataItem->data.V;
     534        PS_PTR_CHECK_NULL(dut,0.0);
     535
     536        t = 2000.0 + (mjd - 51544.03)/365.2422;
     537        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);
     538        result = dut->data.F64[0] + dut->data.F64[1]*(mjd - dut->data.F64[2]) - dut2ut1;
     539
     540    } else if(status != PS_LOOKUP_SUCCESS) {
     541        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
     542        return 0.0;
     543    }
     544
     545    return result;
    610546}
    611547
    612548struct psSphere* psTimeGetPoleCoords(const psTime* time)
    613549{
    614     double x = 0.0;
    615     double y = 0.0;
     550    psU32 nTables = 3;
     551    psF64 x = 0.0;
     552    psF64 y = 0.0;
     553    psF64 mjd = 0.0;
     554    psF64 a = 0.0;
     555    psF64 c = 0.0;
     556    psF64 mjdPred = 0.0;
    616557    struct psSphere* output = NULL;
    617 
     558    psLookupStatusType xStatus = PS_LOOKUP_SUCCESS;
     559    psLookupStatusType yStatus = PS_LOOKUP_SUCCESS;
     560    psMetadataItem *tableMetadataItem = NULL;
     561    char *metadataTableNames[3] = {"ser7", "eopc",  "finals"};
     562    psVector *xp = NULL;
     563    psVector *yp = NULL;
    618564
    619565    // Error checks
    … …  
    621567    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NULL);
    622568
    623     x = lookupSer7Table((psTime*)time, 4);
    624     y = lookupSer7Table((psTime*)time, 5);
    625 
     569    if(time->type != PS_TIME_TAI)
     570    {
     571        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type);
     572        return NULL;
     573    }
     574
     575    // Attempt to find value through table lookup and interpolation
     576    mjd = psTimeToMJD(time);
     577    x = searchTables(mjd, 0, &xStatus, metadataTableNames, nTables);
     578    y = searchTables(mjd, 0, &yStatus, metadataTableNames, nTables);
     579
     580    // Value could not be found through table lookup and interpolation
     581    if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP)
     582    {
     583
     584        // Date too earlier for tables. Get default polar coodinate values from metadata, and issue warning.
     585        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "polar motion");
     586
     587        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp");
     588        if(tableMetadataItem == NULL) {
     589            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.xp");
     590            return NULL;
     591        }
     592        x = tableMetadataItem->data.F64;
     593
     594        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp");
     595        if(tableMetadataItem == NULL) {
     596            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.yp");
     597            return NULL;
     598        }
     599        y = tableMetadataItem->data.F64;
     600
     601    } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM)
     602    {
     603
     604        /* Date too late for tables. Issue warning and use following formulae for predicting
     605           ahead of the most recent available table entry.
     606              x = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c
     607              [0], [1], [2], [3] = @psLib.time.predict.xp
     608              y = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c
     609              [0], [1], [2], [3] = @psLib.time.predict.yp
     610              a = 2*pi*(mjd - pslib.time.predict.mjd)/365.25
     611              c = 2*pi*(mjd - pslib.time.predict.mjd)/435.0
     612        */
     613        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "polar motion");
     614
     615
     616        // Get predicted MJD
     617        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd");
     618        if(tableMetadataItem == NULL) {
     619            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     620                    "psLib.time.predict.mjd");
     621            return NULL;
     622        }
     623        mjdPred = tableMetadataItem->data.F64;
     624
     625        // Get xp
     626        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp");
     627        if(tableMetadataItem == NULL) {
     628            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp");
     629            return NULL;
     630        }
     631        xp = (psVector*)tableMetadataItem->data.V;
     632        PS_PTR_CHECK_NULL(xp,NULL);
     633
     634        // Get yp
     635        tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp");
     636        if(tableMetadataItem == NULL) {
     637            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp");
     638            return NULL;
     639        }
     640        yp = (psVector*)tableMetadataItem->data.V;
     641        PS_PTR_CHECK_NULL(yp,NULL);
     642
     643        // Calculate "a" and "c" constants
     644        a = TWOPI*(mjd - mjdPred)/365.25;
     645        c = TWOPI*(mjd - mjdPred)/435.0;
     646
     647        // Calculate x and y polar coordinates
     648        x = xp->data.F64[0] +
     649            xp->data.F64[1]*cos(a) +
     650            xp->data.F64[2]*sin(a) +
     651            xp->data.F64[3]*cos(c) +
     652            xp->data.F64[4]*sin(c);
     653
     654        y = yp->data.F64[0] +
     655            yp->data.F64[1]*cos(a) +
     656            yp->data.F64[2]*sin(a) +
     657            yp->data.F64[3]*cos(c) +
     658            yp->data.F64[4]*sin(c);
     659
     660    } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS)
     661    {
     662        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
     663        return NULL;
     664    }
     665
     666    // Create output sphere and convert arcsec to radians (i.e. x/60/60*M_PI/180)
    626667    output = psAlloc(sizeof(psSphere));
    627     output->r = x * M_PI / 648000.0; // x in arcsec converted to radians, i.e., x/60/60*M_PI/180
     668    output->r = x * M_PI / 648000.0;
    628669    output->d = y * M_PI / 648000.0;
    629670
    … …  
    633674double psTimeGetTAIDelta(const psTime *time)
    634675{
     676    psU64 i = 0;
     677    psF64 jd = 0.0;
     678    psF64 mjd = 0.0;
     679    psF64 out = 0.0;
     680    psF64 const1 = 0.0;
     681    psF64 const2 = 0.0;
     682    psF64 const3 = 0.0;
     683    psLookupTable* table = NULL;
     684    psMetadataItem *tableMetadataItem = NULL;
     685    psVector *stats = NULL;
     686    psVector *results = NULL;
    635687
    636688    // Error checks
    … …  
    638690    PS_INT_CHECK_RANGE(time->usec,0,1e6-1,NAN);
    639691
    640     return lookupTaiUtcTable(time);
    641 }
     692    // Check time metadata
     693    if(timeMetadata == NULL) {
     694        if(!psTimeInit(TIME_CONFIG_FILE))
     695            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND, "psTime.config");
     696        return 0.0;
     697    }
     698
     699    // Get table from metadata
     700    tableMetadataItem = psMetadataLookup(timeMetadata, "taiTable");
     701    if(tableMetadataItem == NULL) {
     702        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "taiTable");
     703        return 0.0;
     704    }
     705    table = (psLookupTable*)tableMetadataItem->data.V;
     706    PS_PTR_CHECK_NULL(table,0);
     707
     708    // Determine Julian and modified Julian dates used in table lookup and time delta calculation
     709    jd = psTimeToJD(time);
     710    mjd = psTimeToMJD(time);
     711
     712    stats = psVectorAlloc(table->numCols, PS_TYPE_U32);
     713    results = psLookupTableInterpolateAll(table, jd, stats);
     714
     715    // Check for successful interpolation
     716    for(i=0; i<stats->n; i++) {
     717        if(stats->data.U32[i]!=PS_LOOKUP_PAST_BOTTOM && stats->data.U32[i]!=PS_LOOKUP_PAST_TOP) {
     718            psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED);
     719        }
     720    }
     721
     722    const1 = results->data.F64[0];
     723    const2 = results->data.F64[1];
     724    const3 = results->data.F64[2];
     725    out = const1 + (mjd - const2) * const3;
     726
     727    psFree(stats);
     728    psFree(results);
     729
     730    return out;
     731}
     732
    642733
    643734psS64 psTimeLeapseconds(const psTime *time1, const psTime *time2)
Note: See TracChangeset for help on using the changeset viewer.