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

Timestamp:

Jan 23, 2006, 10:04:31 AM (20 years ago)

Author:

drobbin

Message:

updated PrecessMethod. Fixed MetadataAddStr to error on NULL input. Time_TideUT1Corr?

File:

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

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-18 23:49:06 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 20:04:31 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -1462 +1462 @@
 }
 
-/******************************************************************************
-psSpherePrecess(coords, fromTime, toTime):
- 
-XXX: Use static memory for tmpST.
- *****************************************************************************/
-psSphere *psSpherePrecess(psSphere *coords,
-                          const psTime *fromTime,
-                          const psTime *toTime)
+psSphereRot *psSpherePrecess(const psTime *fromTime,
+                             const psTime *toTime,
+                             psPrecessMethod mode)
 {
     // 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);
-
+    if (fromTime == NULL && toTime == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                "Invalid time inputs.  fromTime & toTime cannot both be NULL.\n");
+        return NULL;
+    }
+    PS_ASSERT_INT_WITHIN_RANGE(mode, PS_PRECESS_ROUGH, PS_PRECESS_IAU2000A, NULL);
+    psF64 fromMJD, toMJD;
     // 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);
-
+    //If either input time is NULL, assume it to be J2000 -> from SDRS as of rev 18
+    if (fromTime == NULL) {
+        fromMJD = MJD_2000;
+    } else {
+        fromMJD = psTimeToMJD(fromTime);
+    }
+    if (toTime == NULL) {
+        toMJD = MJD_2000;
+    } else {
+        toMJD = psTimeToMJD(toTime);
+    }
+    psTime *from = NULL;
+    psTime *to = NULL;
+
+    if (mode == PS_PRECESS_ROUGH) {
+        //For PS_PRECESS_ROUGH, no time/earthpole corrections are used.  This is the
+        //lowest level of detail mode.
+        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);
+        return tmpST;
+    } else if (mode == PS_PRECESS_IAU2000A) {
+        //For IAU2000A mode, run psEOC_PrecessionModel and then psSphereRot_CEOtoGCRS for
+        //each time.  Then difference the resulting rotations by adding the inverse
+        //rotation corresponding to fromTime to the toTime rotation.
+
+        //Since the time inputs are allowed to be NULL, either convert the MJD time
+        //or copy to non-NULL time
+        if (fromTime == NULL) {
+            from = psTimeFromMJD(fromMJD);
+        } else {
+            from = p_psTimeCopy(fromTime);
+        }
+        if (toTime == NULL) {
+            to = psTimeFromMJD(toMJD);
+        } else {
+            to = p_psTimeCopy(toTime);
+        }
+
+        psEarthPole *fromEP = psEOC_PrecessionModel(from);
+        psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
+        psEarthPole *toEP = psEOC_PrecessionModel(to);
+        psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP);
+        psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot);
+        psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj);
+        psFree(from);
+        psFree(to);
+        psFree(fromEP);
+        psFree(fromRot);
+        psFree(toEP);
+        psFree(toRot);
+        psFree(fromConj);
+        return out;
+    } else if (mode == PS_PRECESS_COMPLETE_A) {
+        //For PS_PRECESS_COMPLETE_A the same procedure as IAU2000A is used but with
+        //additional earthpole corrections from psEOC_PrecessionCorr.  The corrections
+        //for COMPLETE_A come from the IERS Bulletin A.
+
+        //Since the time inputs are allowed to be NULL, either convert the MJD time
+        //or copy to non-NULL time
+        if (fromTime == NULL) {
+            from = psTimeFromMJD(fromMJD);
+        } else {
+            from = p_psTimeCopy(fromTime);
+        }
+        if (toTime == NULL) {
+            to = psTimeFromMJD(toMJD);
+        } else {
+            to = p_psTimeCopy(toTime);
+        }
+
+        psEarthPole *fromEP = psEOC_PrecessionModel(from);
+        psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_A);
+        fromEP->x += fromCorr->x;
+        fromEP->y += fromCorr->y;
+        fromEP->s += fromCorr->s;
+        psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
+        psEarthPole *toEP = psEOC_PrecessionModel(to);
+        psEarthPole *toCorr = psEOC_PrecessionCorr(to, PS_IERS_A);
+        toEP->x += toCorr->x;
+        toEP->y += toCorr->y;
+        toEP->s += toCorr->s;
+        psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP);
+        psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot);
+        psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj);
+        psFree(from);
+        psFree(to);
+        psFree(fromEP);
+        psFree(fromCorr);
+        psFree(fromRot);
+        psFree(toEP);
+        psFree(toCorr);
+        psFree(toRot);
+        psFree(fromConj);
+        return out;
+    } else {  //mode == PS_PRECESS_COMPLETE_B
+        //For PS_PRECESS_COMPLETE_B the same procedure as IAU2000A is used but with
+        //additional earthpole corrections from psEOC_PrecessionCorr.  The corrections
+        //for COMPLETE_B come from the IERS Bulletin B.
+
+        //Since the time inputs are allowed to be NULL, either convert the MJD time
+        //or copy to non-NULL time
+        if (fromTime == NULL) {
+            from = psTimeFromMJD(fromMJD);
+        } else {
+            from = p_psTimeCopy(fromTime);
+        }
+        if (toTime == NULL) {
+            to = psTimeFromMJD(toMJD);
+        } else {
+            to = p_psTimeCopy(toTime);
+        }
+
+        psEarthPole *fromEP = psEOC_PrecessionModel(from);
+        psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_B);
+        fromEP->x += fromCorr->x;
+        fromEP->y += fromCorr->y;
+        fromEP->s += fromCorr->s;
+        psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
+        psEarthPole *toEP = psEOC_PrecessionModel(to);
+        psEarthPole *toCorr = psEOC_PrecessionCorr(to, PS_IERS_B);
+        toEP->x += toCorr->x;
+        toEP->y += toCorr->y;
+        toEP->s += toCorr->s;
+        psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP);
+        psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot);
+        psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj);
+        psFree(from);
+        psFree(to);
+        psFree(fromEP);
+        psFree(fromCorr);
+        psFree(fromRot);
+        psFree(toEP);
+        psFree(toCorr);
+        psFree(toRot);
+        psFree(fromConj);
+        return out;
+    }
     // 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);
+    //    psSphere *out = psSphereRotApply(NULL, tmpST, coords);
+    //    if (out->r < -0.0001) {
+    //        out->r += 2.0 * M_PI;
+    //    }
 }