Changeset 4540 for trunk/psLib/src/astronomy/psCoord.c

Timestamp:

Jul 12, 2005, 9:12:01 AM (21 years ago)

Author:

desonia

Message:

massive restructuring of codebase.

File:

• trunk/psLib/src/astronomy/psCoord.c (modified) (13 diffs)

trunk/psLib/src/astronomy/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.78 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-06-27 20:38:11 $
+*  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-07-12 19:12:00 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -30 +30 @@
 #include <math.h>
 #include <float.h>
-/******************************************************************************/
-/*  DEFINE STATEMENTS                                                         */
-/******************************************************************************/
 
 // Modified Julian Day 01/01/1900 00:00:00
@@ -40 +37 @@
 #define JULIAN_CENTURY 36525.0
 
-/******************************************************************************/
-/*  TYPE DEFINITIONS                                                          */
-/******************************************************************************/
-
-// None
-
-/*****************************************************************************/
-/*  GLOBAL VARIABLES                                                         */
-/*****************************************************************************/
-
-// None
-
-/*****************************************************************************/
-/*  FILE STATIC VARIABLES                                                    */
-/*****************************************************************************/
-
-// None
-
-/*****************************************************************************/
-/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
-/*****************************************************************************/
 static void planeFree(psPlane *p)
 {
     // There are non dynamic allocated items
+}
+
+static void sphereFree(psSphere *s)
+{
+    // There are non dynamic allocated items
+}
+
+static void planeTransformFree(psPlaneTransform *pt)
+{
+    psFree(pt->x);
+    psFree(pt->y);
+}
+
+static void planeDistortFree(psPlaneDistort *pt)
+{
+    psFree(pt->x);
+    psFree(pt->y);
+}
+
+static void projectionFree(psProjection *p)
+{
+    // There are no dynamically allocated items
 }
 
@@ -201 +195 @@
 }
 
-
-static void sphereFree(psSphere *s)
-{
-    // There are non dynamic allocated items
-}
-
 psSphere* psSphereAlloc(void)
 {
@@ -215 +203 @@
 }
 
-static void planeTransformFree(psPlaneTransform *pt)
-{
-    psFree(pt->x);
-    psFree(pt->y);
+psSphereRot* psSphereRotAlloc(double alphaP,
+                              double deltaP,
+                              double phiP)
+{
+    psSphereRot* rot = psAlloc(sizeof(psSphereRot));
+
+    double cosDelta = cos(deltaP);
+    double halfPhi = phiP / 2.0;
+    double sinHalfPhi = sin(halfPhi);
+
+    // equations are directly from ADD
+    double vx = cosDelta*cos(alphaP);
+    double vy = cosDelta*sin(alphaP);
+    double vz = sin(deltaP);
+
+    rot->q0 = vx*sinHalfPhi;
+    rot->q1 = vy*sinHalfPhi;
+    rot->q2 = vz*sinHalfPhi;
+    rot->q3 = cos(halfPhi);
+
+    return rot;
+}
+
+psSphereRot* psSphereRotQuat(double q0,
+                             double q1,
+                             double q2,
+                             double q3)
+{
+    psSphereRot* rot = psAlloc(sizeof(psSphereRot));
+
+    double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
+    rot->q0 = q0 / len;
+    rot->q1 = q1 / len;
+    rot->q2 = q2 / len;
+    rot->q3 = q3 / len;
+
+    return rot;
 }
 
@@ -260 +281 @@
 
     return (out);
-}
-
-static void planeDistortFree(psPlaneDistort *pt)
-{
-    psFree(pt->x);
-    psFree(pt->y);
 }
 
@@ -319 +334 @@
 
 /******************************************************************************
-alpha is LONGITUDE
-delta is LATITUDE
- 
-    alphaP: Take the target pole in the source system; calculate its LONGITUDE
-     in the target system.  That longitude is alphaP.
-    DeltaP: Take the target pole in the source system; calculate its LATITUDE
-     in the target system.  That longitude is deltaP.
-    phiP:   This is the LONGITUDE of the ascending node in the target system.
- *****************************************************************************/
-psSphereTransform* psSphereTransformAlloc(psF64 alphaP,
-        psF64 deltaP,
-        psF64 phiP)
-{
-    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
-
-    tmp->cosDeltaP = cos(deltaP);
-    tmp->sinDeltaP = sin(deltaP);
-    tmp->alphaP = alphaP;
-    tmp->phiP = phiP;
-
-    return (tmp);
-}
-
-/******************************************************************************
-XXX: Private Function.
- 
-piNormalize(): take an input angle in radians and convert it to the range 0:2*PI.
- *****************************************************************************/
-psF32 piNormalize(psF32 angle)
-{
-    while (angle < FLT_EPSILON) {
-        angle+=M_PI*2;
-    }
-
-    while (angle >= (M_PI*2)) {
-        angle-=M_PI*2;
-    }
-    return(angle);
-}
-
-/******************************************************************************
 XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
 XXX: Should we do something for Declination as well?
  *****************************************************************************/
-psSphere* psSphereTransformApply(psSphere* out,
-                                 const psSphereTransform* transform,
-                                 const psSphere* coord)
+psSphere* psSphereRotApply(psSphere* out,
+                           const psSphereRot* transform,
+                           const psSphere* coord)
 {
     PS_ASSERT_PTR_NON_NULL(transform, NULL);
@@ -371 +345 @@
 
     if (out == NULL) {
-        out = (psSphere* ) psAlloc(sizeof(psSphere));
-    }
-
-    psF64 alpha = coord->r;
-    psF64 delta = coord->d;
-    psF64 alphaMinusAlphaP = alpha - transform->alphaP;
-
-    psF64 eq55 = (sin(delta) * transform->cosDeltaP) -
-                 (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));
-    psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +
-                 (sin(delta) * transform->sinDeltaP);
-    psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);
-
-    psF64 theta = asin(eq55);
-    psF64 phi = atan2(eq56, eq57) + transform->phiP;
-    out->r = piNormalize(phi);
-    out->d = theta;
-
-    return(out);
-}
-
-psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time)
+        out = psSphereAlloc();
+    }
+
+
+    // apply the transform by creating a new psSphereRot from the input coord
+    // and combining it with the input transform (see ADD)
+    psSphereRot* coordRot = psSphereRotAlloc(coord->r, coord->d, 0);
+    coordRot->q3 = 0.0;
+    coordRot = psSphereRotCombine(coordRot, transform, coordRot);
+    // N.B., we can recycle coordRot right away due to the implementation of
+    // psSphereRotCombine puts the values of coordRot in a local variable first
+
+    out->r = atan2(coordRot->q1,coordRot->q0);
+    out->d = atan2(coordRot->q2,sqrt(coordRot->q1*coordRot->q1+coordRot->q0*coordRot->q0));
+
+    return out;
+}
+
+psSphereRot* psSphereRotCombine(psSphereRot* out,
+                                const psSphereRot* rot1,
+                                const psSphereRot* rot2)
+{
+    PS_ASSERT_PTR_NON_NULL(rot1, NULL);
+    PS_ASSERT_PTR_NON_NULL(rot2, NULL);
+
+    if (out == NULL) {
+        out = (psSphereRot* ) psAlloc(sizeof(psSphereRot));
+    }
+
+    double a0 = rot1->q0;
+    double a1 = rot1->q1;
+    double a2 = rot1->q2;
+    double a3 = rot1->q3;
+    double b0 = rot2->q0;
+    double b1 = rot2->q1;
+    double b2 = rot2->q2;
+    double b3 = rot2->q3;
+
+    // following came from ADD
+    out->q0 = b3*a0 + b2*a1 - b1*a2 + b0*a3;
+    out->q1 = b3*a1 - b2*a0 + b1*a3 + b0*a2;
+    out->q2 = b3*a2 + b2*a3 + b1*a0 - b0*a1;
+    out->q3 = b3*a3 - b3*a2 - b1*a1 - b0*a0;
+
+    return out;
+}
+
+psSphereRot *psSphereRotInvert(psSphereRot *rot)
+{
+    PS_ASSERT_PTR_NON_NULL(rot, NULL);
+
+    double norm = sqrt(rot->q0*rot->q0 + rot->q1*rot->q1 + rot->q2*rot->q2 + rot->q3*rot->q3);
+    rot->q1 = -rot->q1 / norm;
+    rot->q2 = -rot->q2 / norm;
+    rot->q3 = -rot->q3 / norm;
+
+    return rot;
+}
+
+psSphereRot* psSphereRotICRSToEcliptic(const psTime *time)
 {
     psF64 T;
@@ -421 +433 @@
 
     // Don't neglect the minus sign on deltaP (bug 244):
-    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
-}
-
-
-psSphereTransform* psSphereTransformEclipticToICRS(psTime *time)
+    return (psSphereRotAlloc(alphaP, deltaP, phiP));
+}
+
+
+psSphereRot* psSphereRotEclipticToICRS(const psTime *time)
 {
     psF64 T;
@@ -453 +465 @@
     psF64 phiP = 0.0;
 
-    return (psSphereTransformAlloc(alphaP, -deltaP, phiP));
+    return (psSphereRotAlloc(alphaP, -deltaP, phiP));
 }
 
 // XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS()
-psSphereTransform* psSphereTransformGalacticToICRS(void)
+psSphereRot* psSphereRotGalacticToICRS(void)
 {
     psF64 alphaP = DEG_TO_RAD(32.93192);
@@ -463 +475 @@
     psF64 phiP = DEG_TO_RAD(282.85948);
 
-    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
-}
-
-psSphereTransform* psSphereTransformICRSToGalactic(void)
+    return (psSphereRotAlloc(alphaP, deltaP, phiP));
+}
+
+psSphereRot* psSphereRotICRSToGalactic(void)
 {
     psF64 alphaP = DEG_TO_RAD(282.85948);
@@ -472 +484 @@
     psF64 phiP = DEG_TO_RAD(32.93192);
 
-    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
-}
-
-void projectionFree(psProjection *p)
-{
-    // There are no dynamically allocated items
+    return (psSphereRotAlloc(alphaP, deltaP, phiP));
 }
 
@@ -837 +844 @@
 
     // Create transform with proper constants
-    psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP);
+    psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP);
 
     // Apply transform to coordinates
-    psSphere *out = psSphereTransformApply(NULL, tmpST, coords);
+    psSphere *out = psSphereRotApply(NULL, tmpST, coords);
 
     psFree(tmpST);