Changeset 5437 for trunk/psLib/src/astro/psSphereOps.c

Timestamp:

Oct 20, 2005, 4:14:02 PM (21 years ago)

Author:

drobbin

Message:

Implemented, debugged, fixed psSphereRot functions/tests. 1 problem remains + testing.

File:

• trunk/psLib/src/astro/psSphereOps.c (modified) (9 diffs)

trunk/psLib/src/astro/psSphereOps.c

@@ -8 +8 @@
  *  @author Dave Robbins, MHPCC
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-18 22:14:47 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-21 02:14:02 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -62 +62 @@
 
     // calculate t*s*r.
-    psSphereRot* result = psSphereRotCombine(NULL,&t,&s);
-    psSphereRotCombine(result,result,&r);
+    psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
+    psSphereRot* result = psSphereRotCombine(NULL, temp, &r);
+    psFree(temp);
 
     return result;
@@ -89 +90 @@
 
     return rot;
+}
+
+psSphereRot* psSphereRotConjugate(psSphereRot *out, const psSphereRot *in)
+{
+    if (in == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_NULL, true,
+                "psSphereRot input cannot be NULL.\n");
+        return NULL;
+    }
+    if (out == NULL) {
+        out = (psSphereRot* ) psAlloc(sizeof(psSphereRot));
+        psMemSetDeallocator(out, (psFreeFunc)sphereRotFree);
+    }
+    out->q0 = -in->q0;
+    out->q1 = -in->q1;
+    out->q2 = -in->q2;
+    out->q3 = in->q3;
+
+    return out;
 }
 
@@ -115 +135 @@
                                  0.0);
 
+    // Inserted by PAP
+    psSphereRot *conjugate = psSphereRotConjugate(NULL, transform);
+    psSphereRot *temp = psSphereRotCombine(NULL, transform, coordQuat);
+    psSphereRot *result = psSphereRotCombine(NULL, temp, conjugate);
+    out->r = atan2(result->q1, result->q0);
+    //    out->r = atan2(result->q1, result->q0);
+    out->d = asin(result->q2);
+    out->rErr = 0.0;
+    out->dErr = 0.0;
+
+    if (out->r < -0.0001) {
+        out->r += 2.0 * M_PI;
+    }
+
+    psFree(conjugate);
+    psFree(temp);
+    psFree(result);
+    psFree(coordQuat);
+    return out;
+    // Pau.
+
+    #if 0
     //    psSphereRot* coordQuatConjugate = psSphereRotQuat(
     //                                       coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3);
@@ -158 +200 @@
 
     return out;
+    #endif
 }
 
@@ -194 +237 @@
 }
 
-psSphereRot *psSphereRotInvert(psSphereRot *rot)
-{
-    if (rot == NULL) {
-        return NULL;
-        // XXX: Error?
-    }
-
-    rot->q0 = -rot->q0;
-    rot->q1 = -rot->q1;
-    rot->q2 = -rot->q2;
-
-    //    rot->q3 = rot->q3;
-    return rot;
-}
-
+psSphereRot *psSphereRotInvert(double alphaP,
+                               double deltaP,
+                               double phiP)
+{
+    return (psSphereRotAlloc(-phiP, -deltaP, -alphaP));
+}
 
 psSphereRot* psSphereRotEclipticToICRS(const psTime *time)
+{
+    psF64 T;
+
+    // Check for null parameter
+    PS_ASSERT_PTR_NON_NULL(time, NULL);
+
+    // Convert psTime to MJD
+    psF64 MJD = psTimeToMJD(time);
+
+    // Check the specified MJD is greater than 1900
+    if ( MJD < MJD_1900 ) {
+        psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD);
+        return NULL;
+    }
+
+    // Calculate number of Julian centuries since 1900
+    T = ( MJD - MJD_1900 ) / JULIAN_CENTURY;
+
+    psF64 phiP = - DEG_TO_RAD(270.0);
+    psF64 deltaP = - (DEG_TO_RAD(23.0) +
+                      MIN_TO_RAD(27.0) +
+                      SEC_TO_RAD(8.26) -
+                      (SEC_TO_RAD(46.845) * T) -
+                      (SEC_TO_RAD(0.0059) * T * T) +
+                      (SEC_TO_RAD(0.00181) * T * T * T));
+    psF64 alphaP = - DEG_TO_RAD(90.0);
+
+    return (psSphereRotAlloc(alphaP, deltaP, phiP));
+}
+
+psSphereRot* psSphereRotICRSToEcliptic(const psTime *time)
 {
     psF64 T;
@@ -241 +306 @@
 }
 
-psSphereRot* psSphereRotICRSToEcliptic(const psTime *time)
-{
-    return psSphereRotInvert(psSphereRotEclipticToICRS(time));
-}
-
 // XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS()
 psSphereRot* psSphereRotGalacticToICRS(void)
 {
+    /*    psF64 phiP = - DEG_TO_RAD(180.0-192.85948);
+        psF64 deltaP = - DEG_TO_RAD(90.0 - 27.12825);
+        psF64 alphaP = - DEG_TO_RAD(90.0+32.93192);
+    */
     psF64 alphaP = DEG_TO_RAD(180.0-192.85948);
-    psF64 deltaP = DEG_TO_RAD(90.0-62.87175);
+    psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825);
     psF64 phiP = DEG_TO_RAD(90.0+32.93192);
+    return (psSphereRotAlloc(-phiP,-deltaP,-alphaP));
+}
+
+psSphereRot* psSphereRotICRSToGalactic(void)
+{
+    psF64 alphaP = DEG_TO_RAD(180.0-192.85948);
+    psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825);
+    psF64 phiP = DEG_TO_RAD(90.0+32.93192);
 
     return (psSphereRotAlloc(alphaP, deltaP, phiP));
-}
-
-psSphereRot* psSphereRotICRSToGalactic(void)
-{
-    return psSphereRotInvert(psSphereRotGalacticToICRS());
 }
 
@@ -479 +546 @@
     // 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) +
+    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));
@@ -497 +564 @@
     // Apply transform to coordinates
     psSphere *out = psSphereRotApply(NULL, tmpST, coords);
+    if (out->r < -0.0001) {
+        out->r += 2.0 * M_PI;
+    }
 
     psFree(tmpST);