Changeset 5507

Timestamp:

Nov 11, 2005, 5:37:34 PM (20 years ago)

Author:

drobbin

Message:

Implemented PrecessionCorr function. Not verified, very possibly wrong.

Location:

trunk/psLib

Files:

• src/astro/psEarthOrientation.c (modified) (6 diffs)
• src/astro/psTime.c (modified) (1 diff)
• src/astro/psTime.h (modified) (2 diffs)
• test/astro/tst_psEarthOrientation.c (modified) (3 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-10 00:13:50 $
+*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-12 03:37:34 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -487 +487 @@
 }
 
-
 psEarthPole *psEOC_PrecessionCorr(const psTime *time,
                                   psTimeBulletin bulletin)
 {
-    return NULL;
+    PS_ASSERT_PTR_NON_NULL(time,NULL);
+    PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL);
+
+    int tableColumn;
+    double mjd;
+    psLookupStatusType status;
+    char* tableName[1] = {"dailyTable"};
+    double x, y, s;
+    //    psMetadataItem *tableMetadataItem = NULL;
+    //   psVector *temp = NULL;
+    //    double calc = 0.0;
+    //    double t = 0.0;
+    //    psMetadata *timeMetadata = p_psTimeGetMetadata();
+
+    mjd = psTimeToMJD(time);
+
+    if (bulletin == PS_IERS_A) {
+        tableColumn = 1;
+    } else {
+        tableColumn = 4;
+    }
+
+    x = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
+    if ( status == PS_LOOKUP_ERROR ) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Time table search returned unsuccessful status.\n");
+        return NULL;
+    } else if (status != PS_LOOKUP_SUCCESS) {
+        psSphere *sphere = p_psTimeGetPoleCoords(time);
+        x = sphere->r;
+        psFree(sphere);
+    } else {
+        x = SEC_TO_RAD(x);
+    }
+
+
+    /*    }else if(status == PS_LOOKUP_PAST_TOP) {
+            // Date too early for tables. Get default time delta value from metadata, and issue warning.
+            psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"xp");
+     
+            // Lookup value from time metadata loaded from pslib.config
+            tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp");
+            if(tableMetadataItem == NULL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                        "psLib.time.before.xp");
+                return NULL;
+            }
+            x = tableMetadataItem->data.F64;
+     
+        } else if(status == PS_LOOKUP_PAST_BOTTOM) {
+            tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd");
+            if(tableMetadataItem == NULL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                        "psLib.time.predict.mjd");
+                return NULL;
+            }
+            mjdPred = tableMetadataItem->data.F64;
+     
+            // Get xp
+            tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp");
+            if(tableMetadataItem == NULL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp");
+                return NULL;
+            }
+            xp = (psVector*)tableMetadataItem->data.V;
+            PS_ASSERT_PTR_NON_NULL(xp,NULL);
+     
+            // Get yp
+            tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp");
+            if(tableMetadataItem == NULL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp");
+                return NULL;
+            }
+            yp = (psVector*)tableMetadataItem->data.V;
+            PS_ASSERT_PTR_NON_NULL(yp,NULL);
+     
+            // Calculate "a" and "c" constants
+            a = TWOPI*(mjd - mjdPred)/365.25;
+            c = TWOPI*(mjd - mjdPred)/435.0;
+     
+            // Calculate x and y polar coordinates
+            x = xp->data.F64[0] +
+                xp->data.F64[1]*cos(a) +
+                xp->data.F64[2]*sin(a) +
+                xp->data.F64[3]*cos(c) +
+                xp->data.F64[4]*sin(c);
+     
+     
+     
+     
+     
+        }
+    */
+
+
+    tableColumn++;
+    y = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
+    if ( status == PS_LOOKUP_ERROR ) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Time table search returned unsuccessful status.\n");
+        return NULL;
+    } else if (status != PS_LOOKUP_SUCCESS) {
+        psSphere *sphere = p_psTimeGetPoleCoords(time);
+        y = sphere->d;
+        psFree(sphere);
+    } else {
+        y = SEC_TO_RAD(y);
+    }
+
+    tableColumn++;
+    s = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
+    if ( status == PS_LOOKUP_ERROR ) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Time table search returned unsuccessful status.\n");
+        return NULL;
+    } else if (status != PS_LOOKUP_SUCCESS) {
+        s = psTimeGetUT1Delta(time, bulletin);
+    } else {
+        s = SEC_TO_RAD(s);
+    }
+
+    //XXX: Do I need to do something for nutation terms w/periods of less than 2 days???
+    psEarthPole *out = psEarthPoleAlloc();
+    out->x = x;
+    out->y = y;
+    out->s = s;
+
+    //    psFree(tableMetadataItem);
+    //    psFree(timeMetadata);
+    return out;
 }
 
@@ -566 +694 @@
 }
 
-
 psSphereRot* psSphereRot_TEOtoCEO(const psTime *time)
 {
@@ -589 +716 @@
                                   psTimeBulletin bulletin)
 {
-
-    return NULL;
+    PS_ASSERT_PTR_NON_NULL(time, NULL);
+
+    psEarthPole *out = psEarthPoleAlloc();
+    psSphere *in = NULL;
+    double s;
+
+    //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate
+    //polar motion coordinates.
+    in = p_psTimeGetPoleCoords(time);
+    if (in == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n");
+        return NULL;
+    }
+    out->x = in->r;
+    out->y = in->d;
+
+    s = psTimeGetUT1Delta(time, bulletin);
+    out->s = s;
+
+    //XXX: Apply polar tide correction here???
+    psEarthPole *correction = psEOC_PolarTideCorr(time);
+    out->x += correction->x;
+    out->y += correction->y;
+    out->s += correction->s;
+
+
+    psFree(in);
+    psFree(correction);
+    return out;
 }
 
@@ -711 +866 @@
         diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
         diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
-     
+
         maxi = 0;
         for (int i = 1; i < 4; ++i) {
@@ -718 +873 @@
             }
         }
-     
+
         double p = 0.5 * sqrt(diag_sum[maxi]);
         recip = 1.0 / (4.0 * p);
-     
+
         if (maxi == 0) {
             out->q0 = p;

trunk/psLib/src/astro/psTime.c

@@ -10 +10 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-05 03:51:41 $
+ *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-12 03:37:34 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

trunk/psLib/src/astro/psTime.h

@@ -11 +11 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-03 04:31:22 $
+ *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-12 03:37:34 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -375 +375 @@
 char* p_psGetConfigFileName();
 
-/**
- *
- *
- *
+/** Searches the IERS time tables for a specified entry location.
+ *
+ *  Returns the interpolated double precision (arcsec) value at the specified entry
+ *  location.  Inputs to specify are the time index in mjd, the column number
+ *  corresponding to Xp, Yp, or Sp (UT1-UTC) in IERS A or B, the time table names,
+ *  and the number of time tables.
+ *
+ *  @return psF64:          Resulting table entry for specified parameters.
  */
 psF64 p_psTimeSearchTables(
-    psF64 index,                       ///<
-    psU64 column,                      ///<
-    char *metadataTableNames[],        ///<
-    psU32 nTables,                     ///<
-    psLookupStatusType* status         ///<
+    psF64 index,                       ///< time index for which to search
+    psU64 column,                      ///< column number of specified index
+    char *metadataTableNames[],        ///< names of IERS tables to search
+    psU32 nTables,                     ///< number of IERS tables to search
+    psLookupStatusType* status         ///< status of table search
 );
 

trunk/psLib/test/astro/tst_psEarthOrientation.c

@@ -5 +5 @@
 *  @author d-Rob, MHPCC
 *
-*  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-10 00:13:51 $
+*  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-12 03:37:34 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -139 +139 @@
 psS32 testEOCPrecession(void)
 {
-
+    psTime *empty = NULL;
+    psTime *time = psTimeAlloc(PS_TIME_TAI);
+    time->sec = 1131579114;
+    time->nsec = 498489000;
+    time->leapsecond = false;
+
+
+    //Tests for Precession Correction function//
+    //Return NULL for NULL time input
+    psEarthPole *pcorr = NULL;
+    pcorr = psEOC_PrecessionCorr(empty, PS_IERS_A);
+    if (pcorr != NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psEOC_PrecessionCorr failed to return NULL for NULL time input.\n");
+        return 4;
+    }
+
+    //Return NULL for Invalid IERS table
+    pcorr = psEOC_PrecessionCorr(time, 3);
+    if (pcorr != NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psEOC_PrecessionCorr failed to return NULL for incorrect IERS table.\n");
+        return 5;
+    }
+    psFree(pcorr);
+
+    //Check values from IERS table A
+    pcorr = psEOC_PrecessionCorr(time, PS_IERS_A);
+    if ( pcorr == NULL ) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, false,
+                "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
+        return 6;
+    } else {
+        printf("Precession Correction output (IERSA) = x,y,s = %.8g, %.8g, %.8g\n",
+               pcorr->x, pcorr->y, pcorr->s);
+    }
+    psFree(pcorr);
+
+    //Check values from IERS table B
+    pcorr = psEOC_PrecessionCorr(time, PS_IERS_B);
+    if ( pcorr == NULL ) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, false,
+                "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
+        return 7;
+    } else {
+        printf("Precession Correction output (IERSB) = x,y,s = %.8g, %.8g, %.8g\n",
+               pcorr->x, pcorr->y, pcorr->s);
+    }
+    psFree(pcorr);
+
+    psFree(time);
     return 0;
 }
@@ -167 +217 @@
     }
 
+    psEarthPole *polarMotion = NULL;
+    polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B);
+    if (polarMotion != NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
+        return 3;
+    }
+    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
+    if (polarMotion == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, false,
+                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
+        return 4;
+    }
+
+    psSphere *interm = p_psTimeGetPoleCoords(in);
+    printf(" -- sphere = r=%.8g, d=%.8g\n", interm->r, interm->d);
     printf(" -- Eop = x=%.8g, y=%.8g, s=%.8g\n", eop->x, eop->y, eop->s);
+    printf(" -- PolarMotion = x=%.8g, y=%.8g, s=%.8g\n", polarMotion->x,
+           polarMotion->y, polarMotion->s);
     psFree(in);
     psFree(eop);
+    psFree(polarMotion);
+    psFree(interm);
     return 0;
 }