Index: trunk/psLib/src/astro/psSphereOps.c =================================================================== --- trunk/psLib/src/astro/psSphereOps.c (revision 6039) +++ trunk/psLib/src/astro/psSphereOps.c (revision 6309) @@ -8,6 +8,6 @@ * @author Dave Robbins, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-18 23:49:06 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-02 23:19:58 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -44,6 +44,6 @@ { psSphereRot r,s,t; - - // directly from ADD -- there must be a better way?! + //The following represents: r =a rotation about the z-axis by alphaP, s =a rotation about + // the y-axis by deltaP, and t =a rotation about the z-axis by phiP. r.q0=0; r.q1=0; @@ -71,8 +71,10 @@ bool psMemCheckSphereRot(psPtr ptr) { + //See if the ptr corresponds to a psSphereRot* return ( psMemGetDeallocator(ptr) == (psFreeFunc)sphereRotFree ); } - +//This function is really a second allocate function for psSphereRot's that uses quaternion +// component inputs instead of angle inputs as in psSphereRotAlloc. psSphereRot* psSphereRotQuat(double q0, double q1, @@ -80,7 +82,10 @@ double q3) { + //allocate space for a new sphere rotation and set deallocator psSphereRot* rot = (psSphereRot*)psAlloc(sizeof(psSphereRot)); psMemSetDeallocator(rot, (psFreeFunc)sphereRotFree); + //The magnitude of a rotation quaternion should = 1 so we normalize here in case the + // inputs have not been. double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); rot->q0 = q0 / len; @@ -92,6 +97,8 @@ } -psSphereRot* psSphereRotConjugate(psSphereRot *out, const psSphereRot *in) -{ +psSphereRot* psSphereRotConjugate(psSphereRot *out, + const psSphereRot *in) +{ + //if input sphere rotation is NULL, return NULL if (in == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, @@ -99,8 +106,11 @@ return NULL; } + //if output sphere rotation is NULL, allocate a new sphere rotation to return if (out == NULL) { out = (psSphereRot* ) psAlloc(sizeof(psSphereRot)); psMemSetDeallocator(out, (psFreeFunc)sphereRotFree); } + //Since q3 is the magnitude and q0,q1,q2 are direction, the conjugate rotation + // is formed by -q0,-q1,-q2,+q3. Note that this is the same as q0,q1,q2,-q3. out->q0 = -in->q0; out->q1 = -in->q1; @@ -111,20 +121,16 @@ } -/****************************************************************************** -XXX: We convert Right Ascension angles to the range 0:PI. Is that acceptable? -XXX: Should we do something for Declination as well? - *****************************************************************************/ psSphere* psSphereRotApply(psSphere* out, const psSphereRot* transform, const psSphere* coord) { + //Make sure that input coordinates and rotation are not NULL PS_ASSERT_PTR_NON_NULL(transform, NULL); PS_ASSERT_PTR_NON_NULL(coord, NULL); - + //If output coordinates is NULL allocate a new psSphere if (out == NULL) { out = psSphereAlloc(); } - // double r = -coord->r; - // apply the transform by creating a new psSphereRot from the input coord + //apply the transform by creating a new psSphereRot from the input coord // and combining it with the input transform (see ADD) double cosD = cos(coord->d); @@ -134,27 +140,18 @@ sin(coord->d), 0.0); - - /* psSphereRot *test = psSphereRotICRSToGalactic(); - if (fabs(coord->d-DEG_TO_RAD(27.1283)) < 0.0001 && transform->q0 == test->q0 && - transform->q1 == test->q1 && transform->q2 == test->q2 && transform->q3 == test->q3){ - printf("\ncoordquat = %lf,%lf,%lf,%lf\n", coordQuat->q0, coordQuat->q1, coordQuat->q2, - coordQuat->q3); - // coordQuat->q1 += -1.0; - test->q3 = - test->q3; - *(psSphereRot**)&transform = test; - } - psFree(test); - */ - // Inserted by PAP + //Get the conjugate of the input rotation. Need to calculate r*p*R to apply rotation + // where R = r-conjugate. psSphereRot *conjugate = psSphereRotConjugate(NULL, transform); psSphereRot *temp = psSphereRotCombine(NULL, transform, coordQuat); psSphereRot *result = psSphereRotCombine(NULL, temp, conjugate); + //From ADD we can find the new sphere coordinates, + // r is calculated as tan^-1 (q1/q0), d is sin^-1 (q2) 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) { + //Simply for convention, we make sure all output r-parameters are positive in the range + // of 0 to 2pi + if (out->r < -0.000001) { out->r += 2.0 * M_PI; } @@ -165,51 +162,4 @@ psFree(coordQuat); return out; - // Pau. - - #if 0 - // psSphereRot* coordQuatConjugate = psSphereRotQuat( - // coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3); - // coordQuat = psSphereRotInvert(coordQuat); - - // calculate q=(rp)r' - // coordQuat = psSphereRotCombine(coordQuat, transform, coordQuat); - // coordQuat = psSphereRotCombine(coordQuat, coordQuat, coordQuatConjugate); - // N.B., we can recycle coordQuat right away due to the implementation of - // psSphereRotCombine; it puts the input values in a local variable first - - // out->r = atan2(coordQuat->q1,coordQuat->q0); - // out->d = asin(coordQuat->q2); - - - // psSphereRot *inv = psSphereRotInvert(transform); - psSphereRot *inv = (psSphereRot*)psAlloc(sizeof(psSphereRot)); - *inv = *transform; - // psSphereRot *result = psSphereRotCombine(NULL, transform, coordQuat); - // inv = psSphereRotInvert(inv); - // psSphereRotCombine(result, result, inv); - // psSphereRot *result = psSphereRotCombine(NULL, inv, coordQuat); - // psSphereRotCombine(result, result, transform); - // psSphereRot *result = psSphereRotCombine(NULL, coordQuat, transform); - // psSphereRotCombine(result, inv, result); - - psSphereRot *result = psSphereRotCombine(NULL, coordQuat, transform); - inv = psSphereRotInvert(inv); - // *(psSphereRot**)&transform = psSphereRotInvert( *(psSphereRot**)&transform); - result = psSphereRotCombine(result, inv, result); - // *(psSphereRot**)&transform = psSphereRotInvert( *(psSphereRot**)&transform); - - // out->r = atan2(result->q0, result->q1); - out->r = atan2(result->q1, result->q0); - out->d = asin(result->q2); - out->rErr = 0.0; - out->dErr = 0.0; - psFree(inv); - psFree(result); - - psFree(coordQuat); - // psFree(coordQuatConjugate); - - return out; - #endif } @@ -218,7 +168,8 @@ const psSphereRot* rot2) { + //Make sure that input rotations are not NULL PS_ASSERT_PTR_NON_NULL(rot1, NULL); PS_ASSERT_PTR_NON_NULL(rot2, NULL); - + //If output rotation is NULL, allocate a new psSphereRot to return if (out == NULL) { out = (psSphereRot* ) psAlloc(sizeof(psSphereRot)); @@ -235,12 +186,8 @@ 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; + //Combine rot1 & rot2. Formulas here came from ADD. out->q0 = a3*b0 + a0*b3 + a1*b2 - a2*b1; out->q1 = a3*b1 - a0*b2 + a1*b3 + a2*b0; out->q2 = a3*b2 + a0*b1 - a1*b0 + a2*b3; - out->q3 = b3*a3 - b2*a2 - b1*a1 - b0*a0; @@ -252,4 +199,6 @@ double phiP) { + //This function should produce identical results to psSphereRotConjugate in most + // or possibly all situations. Creates the inverse rotation from the input angles. return (psSphereRotAlloc(-phiP, -deltaP, -alphaP)); } @@ -258,11 +207,8 @@ { 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 ) { @@ -270,8 +216,7 @@ return NULL; } - // Calculate number of Julian centuries since 1900 T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; - + //Formulas for phiP, deltaP, alphaP came from ADD. psF64 phiP = - DEG_TO_RAD(270.0); psF64 deltaP = - (DEG_TO_RAD(23.0) + @@ -289,11 +234,8 @@ { 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 ) { @@ -301,8 +243,10 @@ return NULL; } - // Calculate number of Julian centuries since 1900 T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; + //Formulas for phiP, deltaP, alphaP came from ADD. Notice that the formulas are the + // same as for EclipticToICRS except that alpha=-phiP, deltaP=-deltaP, phiP=-alphaP to + // produce the inverse rotation as in psSphereRotInverse. psF64 alphaP = DEG_TO_RAD(270.0); psF64 deltaP = DEG_TO_RAD(23.0) + @@ -317,11 +261,8 @@ } -// 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); - */ + //Formulas for alphaP, deltaP, phiP came from the ADD for ICRSToGalactic. Notice that + // this is the reason the inverse gets allocated and returned here. psF64 alphaP = DEG_TO_RAD(180.0-192.85948); psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825); @@ -332,4 +273,5 @@ psSphereRot* psSphereRotICRSToGalactic(void) { + //Formulas for alphaP, deltaP, phiP came from ADD. psF64 alphaP = DEG_TO_RAD(180.0-192.85948); psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825); @@ -339,19 +281,5 @@ } -/****************************************************************************** -The basic idea is to project both positions onto the linear plane, with -position1 at the center, then calculate the linear offset between those -projections. - -XXX: Do I need to check for unacceptable transformation parameters? Maybe, - if the points are on the North/South Pole, etc? - -XXX: Do I need to somehow scale this projection? - -XXX: Does PS_LINEAR mode make sense? The result must be returned in psSphere - regardless of the mode. - -XXX: How to compound errors? - *****************************************************************************/ +//Calculates the difference between coordinates in position1 & 2 and returns this offset. psSphere* psSphereGetOffset(const psSphere* position1, const psSphere* position2, @@ -359,7 +287,7 @@ psSphereOffsetUnit unit) { + //Make sure that input coordinates are not NULL & that mode & unit are valid PS_ASSERT_PTR_NON_NULL(position1, NULL); PS_ASSERT_PTR_NON_NULL(position2, NULL); - // Check positions near 90 degree and issue warnings if necessary if (position1->d >= DEG_TO_RAD(90.0)) { @@ -373,5 +301,4 @@ return NULL; } - // Allocate return structure psSphere* tmp = psSphereAlloc(); @@ -380,4 +307,6 @@ // onto tangent plane, set point projected into psSphere structure x->r y->d if (mode == PS_LINEAR) { + //The basic idea is to project both positions onto the linear plane, with position1 + // at the center, then calculate the linear offset between those projections. psProjection* proj = psProjectionAlloc(position1->r, position1->d, @@ -385,12 +314,9 @@ 1.0, PS_PROJ_TAN); - // Perform projection onto tangent plane psPlane* lin = psProject(position2, proj); - // Set return values tmp->r = lin->x; tmp->d = lin->y; - // Free data structures allocated psFree(proj); @@ -410,5 +336,5 @@ tmp->dErr = 0.0; - // Convert to desired units + // Convert output to desired units if (unit == PS_ARCSEC) { tmp->r = RAD_TO_SEC(tmp->r); @@ -428,7 +354,7 @@ return NULL; } + // Invalid mode } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN, @@ -442,15 +368,5 @@ } -/****************************************************************************** -XXX: Do we need to check for unacceptable transformation parameters? Maybe, - if the points are on the North/South Pole, etc? - -XXX: Do we need to somehow scale this projection? - -XXX: I copied the algorithm from the ADD exactly. - -XXX: Should we compound errors? - *****************************************************************************/ - +//Applies the specified offset to the input position coordinates and returns the new position. psSphere* psSphereSetOffset(const psSphere* position, const psSphere* offset, @@ -458,4 +374,5 @@ psSphereOffsetUnit unit) { + //Make sure that input coordinates are not NULL & that mode & unit are valid PS_ASSERT_PTR_NON_NULL(position, NULL); PS_ASSERT_PTR_NON_NULL(offset, NULL); @@ -469,10 +386,8 @@ // new sphere coordinate if (mode == PS_LINEAR) { - // Allocate plane coordinate and set coordinate psPlane* lin = psPlaneAlloc(); lin->x = offset->r; lin->y = offset->d; - // Allocate and set projection structure psProjection* proj = psProjectionAlloc(position->r, @@ -481,8 +396,6 @@ 1.0, PS_PROJ_TAN); - // Project tangent plane coord to spherical coord tmp = psDeproject(lin, proj); - // Free data structures used psFree(proj); @@ -492,5 +405,4 @@ // to the position and wrap to 0 to 2pi } else if (mode == PS_SPHERICAL) { - // Convert offset unit to radians if (unit == PS_ARCSEC) { @@ -526,5 +438,4 @@ // Invalid mode report error } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,