Changeset 2650 for trunk/psLib/src/astro/psCoord.c

Timestamp:

Dec 7, 2004, 12:28:29 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/astro/psCoord.c (modified) (6 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-12-06 19:49:36 $
+*  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-12-07 22:28:28 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -189 +189 @@
 
 /******************************************************************************
-XXX: It's not clear if this code correctly implements what the equations in
-the ADD imply.  Bug 243 will, hopefully, decide that.
- 
-XXX: Once this code is verified, delete Gene's code.
  *****************************************************************************/
 psSphere* psSphereTransformApply(psSphere* out,
@@ -224 +220 @@
 
     return(out);
-
-    /*
-    This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
-    (r,d) in the sphere coordinates.
-     
-    XXX: In Gene's email, there are different variables with similar names (y
-    and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
-    are no typo's.
-     
-        psF64 sinY = 0.0;
-        psF64 cosY = 0.0;
-        psF64 sinX = 0.0;
-        psF64 cosX = 0.0;
-        psF64 x = 0.0;
-        psF64 y = 0.0;
-        psF64 dx = 0.0;
-     
-        x = coord->r;
-        y = coord->d;
-        dx = x - transform->phiP;
-        sinY = cos(y) * sin(dx) * transform->sinDeltaP +
-               sin(y) * transform->cosDeltaP;
-        cosY = sqrt(1.0 - sinY * sinY);
-        sinX = (cos(y) * sin(dx) * transform->cosDeltaP -
-                sin(y) * transform->sinDeltaP) / cos(y);
-        cosX = cos(y) * cos(dx) / cos(y);
-     
-        out->r = atan2(sinX, cosX) + transform->alphaP;
-        out->d = atan2(sinY, cosY);
-    */
-}
-
-// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
-// there is no algorithm for this in the ADD.
-// XXX: This is bug 244
-psSphereTransform* psSphereTransformICRSToEcliptic(psTime time)
-{
-    PS_PTR_CHECK_NULL(&time, NULL);
-
-    struct tm *tmTime = psTimeToTM(&time);
-    psF64 year = (psF64)(1900 + tmTime->tm_year);
-    psF64 T = year / 100.0;
-    psF64 phi = DEG_TO_RAD(-23.452294) +
-                (DEG_TO_RAD(0.013013) * T) +
-                (DEG_TO_RAD(0.000001639) * T * T) -
-                (DEG_TO_RAD(0.000000503) * T * T * T);
-    psF64 alphaP = 0.0;
-    psF64 phiP = 0.0;
-
-    psFree(tmTime);
-    return (psSphereTransformAlloc(phi, alphaP, phiP));
-}
-
-// XXX: the ADD changed.  This code no longer conforms to the old ADD.
-// XXX: Waiting for IfA to sepcify if T should include days/months/etc.
-psSphereTransform* psSphereTransformEclipticToICRS(psTime time)
-{
-    PS_PTR_CHECK_NULL(&time, NULL);
-
-    struct tm *tmTime = psTimeToTM(&time);
-    psF64 year = (psF64)(1900 + tmTime->tm_year);
+}
+
+psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time)
+{
+    PS_PTR_CHECK_NULL(time, NULL);
+    psF64 MJD = psTimeToMJD(time);
+    struct tm *tmTime = psTimeToTM(time);
+    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
     psF64 T = year / 100.0;
     psF64 alphaP = 0.0;
@@ -296 +239 @@
 
     psFree(tmTime);
+
+    // Don't neglect the minus sign on deltaP (bug 244):
+    return (psSphereTransformAlloc(alphaP, -deltaP, phiP));
+}
+
+
+psSphereTransform* psSphereTransformEclipticToICRS(psTime *time)
+{
+    PS_PTR_CHECK_NULL(time, NULL);
+    psF64 MJD = psTimeToMJD(time);
+    struct tm *tmTime = psTimeToTM(time);
+    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
+    psF64 T = year / 100.0;
+    psF64 alphaP = 0.0;
+    psF64 deltaP = DEG_TO_RAD(23.0) +
+                   SEC_TO_RAD(27.0) +
+                   MIN_TO_RAD(8.0) -
+                   (SEC_TO_RAD(46.845) * T) -
+                   (SEC_TO_RAD(0.0059) * T * T) +
+                   (SEC_TO_RAD(0.00181) * T * T * T);
+    psF64 phiP = 0.0;
+
+    psFree(tmTime);
     return (psSphereTransformAlloc(alphaP, deltaP, phiP));
 }
 
-// XXX: Should these be in radians, not degrees?
-// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
-// there is no algorithm for this in the ADD.
-// XXX: This is bug 245
+// XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalaticToICRS()
 psSphereTransform* psSphereTransformICRSToGalatic(void)
 {
-    psF64 alphaP = DEG_TO_RAD(-62.6);
-    psF64 deltaP = DEG_TO_RAD(33.0);
-    psF64 phiP = DEG_TO_RAD(282.25);
+    psF64 alphaP = DEG_TO_RAD(32.93192);
+    psF64 deltaP = DEG_TO_RAD(-62.87175);
+    psF64 phiP = DEG_TO_RAD(282.85948);
 
     return (psSphereTransformAlloc(alphaP, deltaP, phiP));
 }
 
-// XXX: Should these be in radians, not degrees?
-// XXX: the ADD changed.  This code no longer conforms to the old ADD.
 psSphereTransform* psSphereTransformGalaticToICRS(void)
 {
@@ -588 +549 @@
 
 
+
 /******************************************************************************
 psSpherePrecess(coords, fromTime, toTime):
@@ -599 +561 @@
                           const psTime *toTime)
 {
-    psF64 fromMJD = fromTime->sec/86400.0 +
-                    fromTime->usec/86400000000.0 +
-                    40587.0;
-    psF64 toMJD = toTime->sec/86400.0 +
-                  toTime->usec/86400000000.0 +
-                  40587.0;
+    psF64 fromMJD = psTimeToMJD(fromTime);
+    psF64 toMJD = psTimeToMJD(toTime);
     psF64 T = (toMJD - fromMJD) / 36525.0;