Changeset 5507 for trunk/psLib/src/astro/psEarthOrientation.c

Timestamp:

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

Author:

drobbin

Message:

Implemented PrecessionCorr function. Not verified, very possibly wrong.

File:

• trunk/psLib/src/astro/psEarthOrientation.c (modified) (6 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;