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

Timestamp:

Jan 18, 2006, 1:49:06 PM (21 years ago)

Author:

drobbin

Message:

chgd psMetadataType's to psDataType's. Moved SpherePrecess to EOC. checked in verified tests.

File:

• trunk/psLib/src/astro/psEarthOrientation.c (modified) (8 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-18 20:59:31 $
+ *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-18 23:49:06 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -470 +470 @@
     // Convert psTime to MJD
     double MJD = psTimeToMJD(time);
+    if (MJD == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to MJD failed.  Invalid input time.\n");
+        return NULL;
+    }
 
     // Calculate number of Julian centuries since 2000
@@ -630 +635 @@
 
     double MJD = psTimeToMJD(time);
-
-    if (MJD < 41684.0 || MJD > 53334.0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "Invalid time input.  Date, %lf, is out of range (41684-53334)\n", MJD);
-        return out;
-    }
-
+    if (MJD == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to MJD failed.  Invalid input time.\n");
+        return NULL;
+    }
+    /*
+        if (MJD < 41684.0 || MJD > 53334.0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Invalid time input.  Date, %lf, is out of range (41684-53334)\n", MJD);
+            return out;
+        }
+    */
     // Check if EOC data loaded
     if(!eocInitialized) {
@@ -854 +864 @@
     }
     double T = psTimeToJD(in);
+    if (T == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to JD failed.  Invalid input time.\n");
+        return NULL;
+    }
     T += -2451545.0;
     double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
@@ -874 +889 @@
 
     double MJD = psTimeToMJD(time);
-    /*
-        if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) ||
-                (MJD > 53858.0 && bulletin == PS_IERS_A) ) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    "Invalid time input.  Date, %lf, is out of range\n", MJD);
-            return out;
-        }
-    */
+    if (MJD == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to MJD failed.  Invalid input time.\n");
+        return NULL;
+    }
+
     // Check if EOC data loaded
     if(! eocInitialized) {
@@ -1007 +1020 @@
     // Convert psTime to MJD
     double MJD = psTimeToMJD(time);
+    if (MJD == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to MJD failed.  Invalid input time.\n");
+        return NULL;
+    }
 
     // Calculate number of Julian centuries since 2000
@@ -1086 +1104 @@
     // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(time, NULL);
-    if (time->type == PS_TIME_UT1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "Invalid time input.  Time cannot be of type UT1.\n");
-        return NULL;
-    }
-
+    /*    if (time->type == PS_TIME_UT1) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Invalid time input.  Time cannot be of type UT1.\n");
+            return NULL;
+        }
+    */
     // Convert psTime to MJD
     double MJD = psTimeToMJD(time);
+    if (MJD == NAN) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Time conversion to MJD failed.  Invalid input time.\n");
+        return NULL;
+    }
 
     // Calculate number of Julian centuries since 2000
@@ -1438 +1461 @@
     return out;
 }
+
+/******************************************************************************
+psSpherePrecess(coords, fromTime, toTime):
+ 
+XXX: Use static memory for tmpST.
+ *****************************************************************************/
+psSphere *psSpherePrecess(psSphere *coords,
+                          const psTime *fromTime,
+                          const psTime *toTime)
+{
+    // Check input for NULL pointers
+    PS_ASSERT_PTR_NON_NULL(coords, NULL);
+    PS_ASSERT_PTR_NON_NULL(fromTime, NULL);
+    PS_ASSERT_PTR_NON_NULL(toTime, NULL);
+
+    // Calculate Julian centuries
+    psF64 fromMJD = psTimeToMJD(fromTime);
+    psF64 toMJD = psTimeToMJD(toTime);
+    psF64 T = (toMJD - fromMJD) / JULIAN_CENTURY;
+
+    // Calculate conversion constants
+    //    psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
+    psF64 alphaP = DEG_TO_RAD(180.0) + ((DEG_TO_RAD(0.6406161) * T) +
+                                        (DEG_TO_RAD(0.0000839) * T * T) +
+                                        (DEG_TO_RAD(0.000005) * T * T * T));
+
+    psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) -
+                   (DEG_TO_RAD(0.0001185) * T * T) -
+                   (DEG_TO_RAD(0.0000116) * T * T * T);
+
+    //    psF64 phiP = DEG_TO_RAD(90.0) + ((DEG_TO_RAD(0.6406161) * T) +
+    psF64 phiP = DEG_TO_RAD(180.0) + ((DEG_TO_RAD(0.6406161) * T) +
+                                      (DEG_TO_RAD(0.0003041) * T * T) +
+                                      (DEG_TO_RAD(0.0000051) * T * T * T));
+
+    // Create transform with proper constants
+    psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP);
+
+    // Apply transform to coordinates
+    psSphere *out = psSphereRotApply(NULL, tmpST, coords);
+    if (out->r < -0.0001) {
+        out->r += 2.0 * M_PI;
+    }
+
+    psFree(tmpST);
+
+    return(out);
+}