IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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;
Note: See TracChangeset for help on using the changeset viewer.