IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6036


Ignore:
Timestamp:
Jan 18, 2006, 10:59:32 AM (21 years ago)
Author:
drobbin
Message:

Small update to TEO/polarTide. New fxn Time_TideUT1Corr done and tested.

Location:
trunk/psLib
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psEarthOrientation.c

    r6030 r6036  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2006-01-18 00:41:29 $
     10 *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2006-01-18 20:59:31 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    846846    if (tidalCorr != NULL && tidalCorr->s != 0.0) {
    847847        int nsec = in->nsec + (int)(tidalCorr->s * 1e9);
    848         if (nsec < 0.0) {
     848        if (nsec < 0) {
    849849            in->sec += -1;
    850             in->nsec = (int)(1e9) - nsec;
     850            in->nsec = (int)(1e9) + nsec;
    851851        } else {
    852852            in->nsec = nsec;
    … …  
    10091009
    10101010    // Calculate number of Julian centuries since 2000
    1011     //XXX: NOT SURE IF THIS IS CORRECT FOR THIS SITUATION
    1012     //    double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
    10131011    double RJD = MJD;
    10141012

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

    r5684 r6036  
    1010 *  @author Ross Harman, MHPCC
    1111 *
    12  *  @version $Revision: 1.76 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-12-05 22:00:48 $
     12 *  @version $Revision: 1.77 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-01-18 20:59:31 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    9696/** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null
    9797 *  terminated copy of the original input string. */
    98 static char *cleanString(char *inString, int sLen)
     98static char *cleanString(char *inString,
     99                         int sLen)
    99100{
    100101    char *ptrB = NULL;
    … …  
    126127/** Returns cleaned token based on delimiter, but not including delimiter. Also changes the pointer location
    127128 * the beginning of the string. Tokens are newly allocated null terminated strings. */
    128 static char* getToken(char **inString, char *delimiter, psParseErrorType *status)
     129static char* getToken(char **inString,
     130                      char *delimiter,
     131                      psParseErrorType *status)
    129132{
    130133    char *cleanToken = NULL;
    … …  
    169172// Searches time tables in priority order and performs interpolation if input index value is within a table.
    170173// If the index value is out of range, the status is set accordingly.
    171 psF64 p_psTimeSearchTables(psF64 index, psU64 column, char *metadataTableNames[],
    172                            psU32 nTables, psLookupStatusType* status)
     174psF64 p_psTimeSearchTables(psF64 index,
     175                           psU64 column,
     176                           char *metadataTableNames[],
     177                           psU32 nTables,
     178                           psLookupStatusType* status)
    173179{
    174180    char*            tableName          = NULL;
    … …  
    660666}
    661667
    662 psTime* psTimeConvert(psTime *time, psTimeType type)
     668psTime* psTimeConvert(psTime *time,
     669                      psTimeType type)
    663670{
    664671    // Error checks
    … …  
    745752}
    746753
    747 double psTimeToLMST(psTime *time, double longitude)
     754double psTimeToLMST(psTime *time,
     755                    double longitude)
    748756{
    749757    psF64  jdTdtDays    =  0.0;
    … …  
    818826}
    819827
    820 double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin)
     828double psTimeGetUT1Delta(const psTime *time,
     829                         psTimeBulletin bulletin)
    821830{
    822831    psU32              nTables               = 2;
    … …  
    899908
    900909    return result;
     910}
     911
     912static double DMOD(double x, double y)
     913{
     914    double value = x - y * trunc(x/y);
     915    return value;
     916}
     917
     918psTime *psTime_TideUT1Corr(const psTime *time)
     919{
     920    PS_ASSERT_PTR_NON_NULL(time, NULL);
     921    psTime *out = NULL;
     922
     923    // Convert psTime to MJD
     924    double MJD = psTimeToMJD(time);
     925    if (MJD == NAN) {
     926        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     927                "Time conversion to MJD failed.  Invalid input time.\n");
     928        return NULL;
     929    }
     930
     931    // Calculate number of Julian centuries since 2000
     932    double RJD = MJD;
     933
     934    //Formula comes from fortran reference
     935    //DMOD in fortran ref. = double remainder -> x - y * trunc(x/y)
     936    double T, L, LPRIME, CAPF, CAPD, OMEGA, THETA, CORZ;
     937    double ARG1, ARG2, ARG3, ARG4, ARG5, ARG6, ARG7, ARG8;
     938    double T2, T3, T4;
     939    T = (RJD - 51544.5) / 36525.0;
     940    T2 = T*T;
     941    T3 = T*T*T;
     942    T4 = T*T*T*T;
     943    L = -0.0002447 * T4 + 0.051635 * T3 + 31.8792 * T2 + 1717915923.2178 * T + 485868.249036;
     944    L = DMOD(L, 1296000.0);
     945    LPRIME = -0.00001149 * T4 - 0.000136 * T3 - 0.5532 * T2 + 129596581.0481 * T + 1287104.79305;
     946    LPRIME = DMOD(LPRIME, 1296000.0);
     947    CAPF = 0.00000417 * T4 - 0.001037 * T3 - 12.7512 * T2 + 1739527262.8478 * T + 335779.526232;
     948    CAPF = DMOD(CAPF, 1296000.0);
     949    CAPD = -0.00003169 * T4 + 0.006593 * T3 - 6.3706 * T2 + 1602961601.209 * T + 1072260.70369;
     950    CAPD = DMOD(CAPD, 1296000.0);
     951    OMEGA = -0.00005939 * T4 + 0.007702 * T3 + 7.4722 * T2 - 6962890.2665 * T + 450160.398036;
     952    OMEGA = DMOD(OMEGA, 1296000.0);
     953    THETA = (67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * T + 0.093104 * T2 -
     954             6.2e-6 * T3) * 15.0 + 648000.0;
     955    ARG7 = DMOD((-L - 2.0 * CAPF - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI)
     956           - M_PI / 2.0;
     957    ARG1 = DMOD((-2.0 * CAPF - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI) - M_PI / 2.0;
     958    ARG2 = DMOD((-2.0 * CAPF + 2.0 * CAPD - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI)
     959           - M_PI / 2.0;
     960    ARG3 = DMOD(THETA * M_PI / 648000.0, 2.0 * M_PI) - M_PI / 2.0;
     961    ARG4 = DMOD((-L - 2.0 * CAPF - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     962    ARG5 = DMOD((-2.0 * CAPF - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     963    ARG6 = DMOD((-2.0 * CAPF + 2.0 * CAPD - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0,
     964                2.0 * M_PI);
     965    ARG8 = DMOD((2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     966    CORZ =  0.0245 * sin(ARG7) + 0.0503 * cos(ARG7)
     967            +0.1210 * sin(ARG1) + 0.1605 * cos(ARG1)
     968            +0.0286 * sin(ARG2) + 0.0516 * cos(ARG2)
     969            +0.0864 * sin(ARG3) + 0.1771 * cos(ARG3)
     970            -0.0380 * sin(ARG4) - 0.0154 * cos(ARG4)
     971            -0.1617 * sin(ARG5) - 0.0720 * cos(ARG5)
     972            -0.0759 * sin(ARG6) - 0.0004 * cos(ARG6)
     973            -0.0196 * sin(ARG8) - 0.0038 * cos(ARG8);
     974    CORZ = CORZ * 0.1e-3;
     975
     976    double timeCheck = (double)(time->sec) + (double)(1e-9*time->nsec);
     977    if ( (timeCheck + CORZ) < 0.0 ) {
     978        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     979                "Invalid time for Tide Correction.\n");
     980        return NULL;
     981    }
     982    out = psTimeAlloc(time->type);
     983    *out = *time;
     984    //    out->sec = time->sec;
     985    //    out->nsec = time->nsec;
     986    //    out->leapsecond = time->leapsecond;
     987    if (out->type != PS_TIME_UT1) {
     988        out = psTimeConvert(out, PS_TIME_UT1);
     989    }
     990    if (fabs(CORZ) > 1.0) {
     991        int sec = (int)CORZ;
     992        out->sec += sec;
     993        int nsec = (int)((CORZ - sec)*1e9);
     994        out->nsec += nsec;
     995    } else {
     996        int nsec = out->nsec + (int)(CORZ * 1e9);
     997        if (nsec < 0) {
     998            out->sec += -1;
     999            out->nsec = (int)(1e9) + nsec;
     1000        } else {
     1001            out->nsec = nsec;
     1002        }
     1003    }
     1004    return out;
    9011005}
    9021006
    … …  
    10901194}
    10911195
    1092 long psTimeLeapSecondDelta(const psTime *time1, const psTime *time2)
     1196long psTimeLeapSecondDelta(const psTime *time1,
     1197                           const psTime *time2)
    10931198{
    10941199    psS64 diff = 0;
    … …  
    13971502}
    13981503
    1399 psTime* psTimeFromTT(psS64 sec, psU32 nsec)
     1504psTime* psTimeFromTT(psS64 sec,
     1505                     psU32 nsec)
    14001506{
    14011507    psTime*      outTime  = NULL;
    … …  
    14151521}
    14161522
    1417 psTime* psTimeFromUTC(psS64 sec, psU32 nsec, bool leapsecond)
     1523psTime* psTimeFromUTC(psS64 sec,
     1524                      psU32 nsec,
     1525                      bool leapsecond)
    14181526{
    14191527    psTime*   outTime   = NULL;
    … …  
    15191627}
    15201628
    1521 psTime* psTimeMath(const psTime *time, double delta)
     1629psTime* psTimeMath(const psTime *time,
     1630                   double delta)
    15221631{
    15231632    psF64 sec = 0.0;
    … …  
    15601669}
    15611670
    1562 double psTimeDelta(const psTime *time1, const psTime *time2)
     1671double psTimeDelta(const psTime *time1,
     1672                   const psTime *time2)
    15631673{
    15641674    psF64 out = 0.0;

  • trunk/psLib/src/astro/psTime.h

    r5507 r6036  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-11-12 03:37:34 $
     13 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2006-01-18 20:59:31 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    7676 * Reads config and data files associated with various time conversions.
    7777 *
    78  * @return  bool: True for success, false for failure.
     78 * @return bool:    True for success, false for failure.
    7979 */
    8080psBool p_psTimeInit(
    … …  
    9494 * Frees time data to be held in memory until the end of successful program execution.
    9595 *
    96  * @return  void: void.
     96 * @return void:    void.
    9797 */
    9898psBool p_psTimeFinalize(void);
    … …  
    110110 * of the struct are set to zero.
    111111 *
    112  * @return  psTime*: Struct with empty time.
     112 * @return psTime*:    Struct with empty time.
    113113 */
    114114psTime* psTimeAlloc(
    … …  
    133133 * (PS_TIME_TAI or PS_TIME_UTC) in the argument.
    134134 *
    135  *  @return  psTime*: Struct with current time.
     135 *  @return psTime*:    Struct with current time.
    136136 */
    137137psTime* psTimeGetNow(
    … …  
    143143 *  Converts psTime to UTC, TAI, UT1, or TT time based on the psTimeType argument.
    144144 *
    145  *  @return  psTime*: Pointer to psTime.
     145 *  @return psTime*:    Pointer to psTime.
    146146 */
    147147psTime* psTimeConvert(
    … …  
    155155 *  in UTC format, then it is converted.
    156156 *
    157  *  @return  double: LST Time.
     157 *  @return double:    LST Time.
    158158 */
    159159double psTimeToLMST(
    … …  
    166166 *  This function is necessary to for various SLALIB functions.
    167167 *
    168  *  @return  double: Time difference.
     168 *  @return double:    Time difference.
    169169 */
    170170double psTimeGetUT1Delta(
    … …  
    173173);
    174174
     175/** Provides tidal corrections to UT1-UTC.
     176 *
     177 *  Uses the Ray model of Simon et al.
     178 *
     179 *  @return psTime*:    The corrected time in UT1.
     180 */
     181psTime *psTime_TideUT1Corr(
     182    const psTime *time                 ///< psTime to be corrected.
     183);
     184
    175185/** Determine TAI - UTC from table lookup.
    176186 *
    177187 *  This function is necessary to for various psTime functions.
    178188 *
    179  *  @return  psF64: Time difference.
     189 *  @return psF64:      Time difference.
    180190 */
    181191psF64 p_psTimeGetTAIDelta(
    … …  
    187197 *  Determines the orientation of the polar axis at the given time.
    188198 *
    189  *  @return  psSphere*: Spherical coordinates of Earth's polar axias.
     199 *  @return psSphere*:      Spherical coordinates of Earth's polar axias.
    190200 */
    191201psSphere* p_psTimeGetPoleCoords(
    … …  
    197207 *  Calculates the number of leapseconds between two times.
    198208 *
    199  *  @return  long: leapseconds added between given times
     209 *  @return long:  leapseconds added between given times
    200210 */
    201211long psTimeLeapSecondDelta(
    … …  
    208218 *  Determines if the specified UTC time is a valid leapsecond.
    209219 *
    210  *  @return  bool: valid leap second
     220 *  @return bool:  valid leap second
    211221 */
    212222bool psTimeIsLeapSecond(
    … …  
    219229 *  subtract leapseconds.
    220230 *
    221  *  @return  double: Julian Date (JD) time.
     231 *  @return double:    Julian Date (JD) time.
    222232 */
    223233double psTimeToJD(
    … …  
    229239 *  add or subtract leapseconds.
    230240 *
    231  *  @return  double: Modified Julian Days (MJD) time.
     241 *  @return double:    Modified Julian Days (MJD) time.
    232242 */
    233243double psTimeToMJD(
    … …  
    240250 *  This function does not add or subtract leapseconds.
    241251 *
    242  *  @return  psString:     Pointer null terminated array of chars in ISO time.
     252 *  @return psString:     Pointer null terminated array of chars in ISO time.
    243253 */
    244254psString psTimeToISO(
    … …  
    250260 *  Converts psTime to timeval time. This function does not add or subtract leapseconds.
    251261 *
    252  *  @return  timeval*: timeval struct time.
     262 *  @return timeval*:  timeval struct time.
    253263 */
    254264struct timeval* psTimeToTimeval(

  • trunk/psLib/test/astro/tst_psTime_04.c

    r5216 r6036  
    1111 *      Test F - Attempt to read data file with typo in number
    1212 *      Test G - Free data
     13 *      Test H - Attempt to use all timer functions
     14 *      Test I - Tidal Corrections to UT1-UTC
    1315 *
    1416 *  @author  Ross Harman, MHPCC
    1517 *  @author  Eric Van Alst, MHPCC
     18 *  @author  David Robbins, MHPCC
    1619 *
    17  *  @version $Revision: 1.4 $  $Name: not supported by cvs2svn $
    18  *  @date  $Date: 2005-10-01 02:22:20 $
     20 *  @version $Revision: 1.5 $  $Name: not supported by cvs2svn $
     21 *  @date  $Date: 2006-01-18 20:59:32 $
    1922 *
    2023 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3336static psS32 testTimeInit6(void);
    3437static psS32 testTimer1(void);
     38static psS32 testTideUT1Corr(void);
    3539
    3640testDescription tests[] = {
    … …  
    4246                              {testTimeInit6,-6,"p_psTimeInit",0,false},
    4347                              {testTimer1,-7,"psTimer",0,false},
     48                              {testTideUT1Corr,-8,"psTime_TideUT1Corr",0,false},
    4449                              {NULL}
    4550                          };
    … …  
    136141}
    137142
     143psS32 testTideUT1Corr(void)
     144{
     145    psTime *empty = NULL;
     146    psTime *tide = NULL;
     147    psTime *noTide = psTimeAlloc(PS_TIME_UTC);
     148    noTide->sec = 1049160600;
     149    noTide->nsec = 0;
     150    noTide->leapsecond = false;
     151
     152    empty = psTime_TideUT1Corr(tide);
     153    if (empty != NULL) {
     154        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     155                "psTime_TideUT1Corr failed to return NULL for NULL input time.\n");
     156        return 1;
     157    }
     158
     159    empty = psTime_TideUT1Corr(noTide);
     160    if (empty->sec != 1049160599 || empty->nsec != 656982272) {
     161        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     162                "psTime_TideUT1Corr failed to return correct values.\n");
     163        printf("\nsec = %ld, nsec = %u\n", empty->sec, empty->nsec);
     164        return 2;
     165    }
     166
     167    if (!p_psTimeFinalize()) {
     168        return 3;
     169    }
     170
     171    psFree(empty);
     172    psFree(noTide);
     173
     174    return 0;
     175}
     176
Note: See TracChangeset for help on using the changeset viewer.