Changeset 6184 for trunk/psLib/src/astro

Timestamp:

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

Author:

drobbin

Message:

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

Location:

trunk/psLib/src/astro

Files:

• psEarthOrientation.c (modified) (2 diffs)
• psEarthOrientation.h (modified) (3 diffs)
• psTime.c (modified) (4 diffs)
• psTime.h (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;
+    //    }
 }

trunk/psLib/src/astro/psEarthOrientation.h

@@ -9 +9 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-18 23:49:06 $
+*  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-23 20:04:31 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -40 +40 @@
 typedef enum {
     PS_PRECESS_ROUGH,                  ///< roughest, lowest level of detail
-    PS_PRECESS_COMPLETE,               ///< complete level of detail
+    PS_PRECESS_COMPLETE_A,             ///< complete level of detail using IERS A
+    PS_PRECESS_COMPLETE_B,             ///< complete level of detail using IERS B
     PS_PRECESS_IAU2000A,               ///< highest level of detail
 }
@@ -179 +180 @@
  *  between two times.  The equinoxes shall be Julian equinoxes.
  *
- *  @return psSphere* the resulting spherical rotation
+ *  If NULL is provided for either time, it is assumed to have the reference
+ *  equinox value of J2000.  The mode argument is used to specify the level of
+ *  detail used in the calculation.
+ *
+ *  @return psSphere*:       the resulting spherical rotation
  */
-psSphere* psSpherePrecess(
-    psSphere *coords,                  ///< coordinates (modified in-place)
+psSphereRot* psSpherePrecess(
+    //    psSphere *coords,                  ///< coordinates (modified in-place)
     const psTime *fromTime,            ///< equinox of coords input
-    const psTime *toTime               ///< equinox of coords output
+    const psTime *toTime,              ///< equinox of coords output
+    psPrecessMethod mode               ///< level of detail to use
 );
 

trunk/psLib/src/astro/psTime.c

@@ -10 +10 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.77 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-18 20:59:31 $
+ *  @version $Revision: 1.78 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 20:04:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -920 +920 @@
     PS_ASSERT_PTR_NON_NULL(time, NULL);
     psTime *out = NULL;
+    //Also see psEOC_PolarTideCorr for more info.
 
     // Convert psTime to MJD
@@ -982 +983 @@
     out = psTimeAlloc(time->type);
     *out = *time;
-    //    out->sec = time->sec;
-    //    out->nsec = time->nsec;
-    //    out->leapsecond = time->leapsecond;
     if (out->type != PS_TIME_UT1) {
         out = psTimeConvert(out, PS_TIME_UT1);
     }
+    //see if corrections include seconds or just nano-seconds
+    //nano-seconds must be converted (scaled) to integer values & total cannot be negative
     if (fabs(CORZ) > 1.0) {
         int sec = (int)CORZ;
@@ -1806 +1806 @@
 }
 
+psTime *p_psTimeCopy(const psTime *inTime)
+{
+    PS_ASSERT_PTR_NON_NULL(inTime, NULL);
+    psTime *outTime = psTimeAlloc(inTime->type);
+    *outTime = *inTime;
+    return outTime;
+}
+

trunk/psLib/src/astro/psTime.h

@@ -11 +11 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-18 20:59:31 $
+ *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 20:04:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -432 +432 @@
 psF64 psTimerStop(void);
 
+/** Private function for copying a psTime.
+ *
+ *  @return psTime*:        New copy of existing psTime.
+ */
+psTime *p_psTimeCopy(
+    const psTime *inTime               ///< input time to copy.
+);
+
 /// @}