Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 13838) +++ /trunk/psLib/src/astro/psCoord.c (revision 13839) @@ -1,18 +1,18 @@ /** @file psCoord.c -* -* @brief Contains basic coordinate transformation definitions and operations -* -* This file defines the basic types for astronomical coordinate -* transformation -* -* @ingroup CoordinateTransform -* -* @author GLG, MHPCC -* -* @version $Revision: 1.138 $ $Name: not supported by cvs2svn $ -* @date $Date: 2007-04-04 21:10:03 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * + * @brief Contains basic coordinate transformation definitions and operations + * + * This file defines the basic types for astronomical coordinate + * transformation + * + * @ingroup CoordinateTransform + * + * @author GLG, MHPCC + * + * @version $Revision: 1.139 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-15 01:03:22 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ /******************************************************************************/ /* INCLUDE FILES */ @@ -72,5 +72,5 @@ XXX: below is the code using the standard matrix representation. note that this inversion requires x->nX == 1, y->nY == 1 and x->nY <= 1, y->nX <= 1 - *****************************************************************************/ +*****************************************************************************/ psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) { @@ -79,7 +79,7 @@ PS_ASSERT_PTR_NON_NULL(transform->y, NULL); if ((transform->x->nX < 1) || - (transform->x->nY < 1) || - (transform->y->nX < 1) || - (transform->y->nY < 1)) { + (transform->x->nY < 1) || + (transform->y->nX < 1) || + (transform->y->nY < 1)) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible."); return(NULL); @@ -150,5 +150,5 @@ Why isn't this called p_psIsPlaneTransformLinear()? - *****************************************************************************/ +*****************************************************************************/ bool p_psIsProjectionLinear(psPlaneTransform *transform) { @@ -161,6 +161,6 @@ if (transform->x->coeff[i][j] != 0.0) { if (!(((i == 0) && (j == 0)) || - ((i == 0) && (j == 1)) || - ((i == 1) && (j == 0)))) { + ((i == 0) && (j == 1)) || + ((i == 1) && (j == 0)))) { return(false); } @@ -173,6 +173,6 @@ if (transform->y->coeff[i][j] != 0.0) { if (!(((i == 0) && (j == 0)) || - ((i == 0) && (j == 1)) || - ((i == 1) && (j == 0)))) { + ((i == 0) && (j == 1)) || + ((i == 1) && (j == 0)))) { return(false); } @@ -304,5 +304,5 @@ This transformation takes into account parameters beyond an objects spatial coordinates: term3 and term4 (magnitude and color). - *****************************************************************************/ +*****************************************************************************/ psPlane* psPlaneDistortApply( psPlane* out, @@ -379,26 +379,27 @@ psProjectionClass class = PS_PROJECTION_CLASS_NONE; switch (projection->type) { - case PS_PROJ_LIN: - case PS_PROJ_PLY: - case PS_PROJ_WRP: + case PS_PROJ_LIN: + case PS_PROJ_PLY: + case PS_PROJ_WRP: class = PS_PROJECTION_CLASS_CARTESIAN; break; - case PS_PROJ_TAN: - case PS_PROJ_DIS: - case PS_PROJ_SIN: - case PS_PROJ_STG: - case PS_PROJ_ZEA: - case PS_PROJ_ZPL: + case PS_PROJ_TAN: + case PS_PROJ_TNX: + case PS_PROJ_DIS: + case PS_PROJ_SIN: + case PS_PROJ_STG: + case PS_PROJ_ZEA: + case PS_PROJ_ZPL: class = PS_PROJECTION_CLASS_ZENITHAL; break; - case PS_PROJ_AIT: - case PS_PROJ_PAR: - case PS_PROJ_GLS: - case PS_PROJ_MER: - case PS_PROJ_CAR: + case PS_PROJ_AIT: + case PS_PROJ_PAR: + case PS_PROJ_GLS: + case PS_PROJ_MER: + case PS_PROJ_CAR: class = PS_PROJECTION_CLASS_PSEUDOCYL; break; - default: + default: psAbort("invalid projection"); break; @@ -414,10 +415,10 @@ switch (class) { - case PS_PROJECTION_CLASS_CARTESIAN: + case PS_PROJECTION_CLASS_CARTESIAN: out->x = (coord->r - projection->R); out->y = (coord->d - projection->D); break; - case PS_PROJECTION_CLASS_ZENITHAL: + case PS_PROJECTION_CLASS_ZENITHAL: sinDp = sin(projection->D); cosDp = cos(projection->D); @@ -427,53 +428,53 @@ cosDelta = cos(coord->d); - # if ELIXIR_CODE +# if ELIXIR_CODE sinTheta = cosDelta*cosAlpha*cosDp + sinDelta*sinDp; - sinPhiCT = cosDelta*sinAlpha; - cosPhiCT = cosDelta*cosAlpha*sinDp - sinDelta*cosDp; - # else + sinPhiCT = cosDelta*sinAlpha; + cosPhiCT = cosDelta*cosAlpha*sinDp - sinDelta*cosDp; +# else sinTheta = cosDelta*cosAlpha*cosDp + sinDelta*sinDp; - sinPhiCT = -cosDelta*sinAlpha; - cosPhiCT = -cosDelta*cosAlpha*sinDp + sinDelta*cosDp; - # endif - // Perform the specified projection - switch (projection->type) - { - case PS_PROJ_TAN: - case PS_PROJ_DIS: - // Gnomonic projection - out->x = +sinPhiCT / sinTheta; - out->y = -cosPhiCT / sinTheta; - break; - case PS_PROJ_SIN: - // Othrographic projection - out->x = +sinPhiCT; - out->y = -cosPhiCT; - break; - case PS_PROJ_STG: - // Othrographic projection - out->x = +2*sinPhiCT / (1 + sinTheta); - out->y = -2*cosPhiCT / (1 + sinTheta); - break; - case PS_PROJ_ZEA: - case PS_PROJ_ZPL: - // Zenithal Equal-Area - zeta = M_SQRT2 / sqrt (1 + sinTheta); - out->x = +zeta * sinPhiCT; - out->y = -zeta * cosPhiCT; - break; - default: - psAbort("invalid projection"); - break; - } - break; - - case PS_PROJECTION_CLASS_PSEUDOCYL: + sinPhiCT = -cosDelta*sinAlpha; + cosPhiCT = -cosDelta*cosAlpha*sinDp + sinDelta*cosDp; +# endif + // Perform the specified projection + switch (projection->type) { + case PS_PROJ_TAN: + case PS_PROJ_TNX: // XXX this is not quite right + case PS_PROJ_DIS: + // Gnomonic projection + out->x = +sinPhiCT / sinTheta; + out->y = -cosPhiCT / sinTheta; + break; + case PS_PROJ_SIN: + // Othrographic projection + out->x = +sinPhiCT; + out->y = -cosPhiCT; + break; + case PS_PROJ_STG: + // Othrographic projection + out->x = +2*sinPhiCT / (1 + sinTheta); + out->y = -2*cosPhiCT / (1 + sinTheta); + break; + case PS_PROJ_ZEA: + case PS_PROJ_ZPL: + // Zenithal Equal-Area + zeta = M_SQRT2 / sqrt (1 + sinTheta); + out->x = +zeta * sinPhiCT; + out->y = -zeta * cosPhiCT; + break; + default: + psAbort("invalid projection"); + break; + } + break; + + case PS_PROJECTION_CLASS_PSEUDOCYL: phi = coord->r - projection->R; theta = coord->d - projection->D; switch (projection->type) { - case PS_PROJ_AIT: + case PS_PROJ_AIT: // Hammer-Aitoff projection zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); @@ -481,19 +482,19 @@ out->y = zeta*sin(theta); break; - case PS_PROJ_GLS: + case PS_PROJ_GLS: // projection name? out->x = phi * cos(theta); out->x = theta; break; - case PS_PROJ_PAR: + case PS_PROJ_PAR: // Parabolic projection out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); out->y = M_PI*sin(theta/3.0); break; - case PS_PROJ_CAR: - case PS_PROJ_MER: + case PS_PROJ_CAR: + case PS_PROJ_MER: psAbort("projection not yet implemented"); break; - default: + default: psAbort("invalid projection"); break; @@ -501,5 +502,5 @@ break; - default: + default: psAbort("invalid projection"); break; @@ -535,26 +536,27 @@ psProjectionClass class = PS_PROJECTION_CLASS_NONE; switch (projection->type) { - case PS_PROJ_LIN: - case PS_PROJ_PLY: - case PS_PROJ_WRP: + case PS_PROJ_LIN: + case PS_PROJ_PLY: + case PS_PROJ_WRP: class = PS_PROJECTION_CLASS_CARTESIAN; break; - case PS_PROJ_TAN: - case PS_PROJ_DIS: - case PS_PROJ_SIN: - case PS_PROJ_STG: - case PS_PROJ_ZEA: - case PS_PROJ_ZPL: + case PS_PROJ_TAN: + case PS_PROJ_TNX: + case PS_PROJ_DIS: + case PS_PROJ_SIN: + case PS_PROJ_STG: + case PS_PROJ_ZEA: + case PS_PROJ_ZPL: class = PS_PROJECTION_CLASS_ZENITHAL; break; - case PS_PROJ_AIT: - case PS_PROJ_PAR: - case PS_PROJ_GLS: - case PS_PROJ_CAR: - case PS_PROJ_MER: + case PS_PROJ_AIT: + case PS_PROJ_PAR: + case PS_PROJ_GLS: + case PS_PROJ_CAR: + case PS_PROJ_MER: class = PS_PROJECTION_CLASS_PSEUDOCYL; break; - default: + default: psAbort("invalid projection"); break; @@ -575,10 +577,10 @@ // Perform inverse projection switch (class) { - case PS_PROJECTION_CLASS_CARTESIAN: + case PS_PROJECTION_CLASS_CARTESIAN: out->r = coord->x*projection->Xs + projection->R; out->d = coord->y*projection->Ys + projection->D; break; - case PS_PROJECTION_CLASS_ZENITHAL: + case PS_PROJECTION_CLASS_ZENITHAL: // Remove plate scales x = coord->x*projection->Xs; @@ -588,8 +590,8 @@ cosPhi = (R == 0) ? 1.0 : -y / R; - switch (projection->type) - { - case PS_PROJ_TAN: - case PS_PROJ_DIS: + switch (projection->type) { + case PS_PROJ_TAN: + case PS_PROJ_TNX:// XXX this is not quite right + case PS_PROJ_DIS: // Gnonomic deprojection rho = sqrt (1 + R*R); @@ -597,14 +599,14 @@ cosTheta = R / rho; break; - case PS_PROJ_SIN: + case PS_PROJ_SIN: // Orhtographic deprojection sinTheta = sqrt (1 - R*R); cosTheta = R; break; - case PS_PROJ_STG: + case PS_PROJ_STG: psAbort("STG not defined"); break; - case PS_PROJ_ZEA: - case PS_PROJ_ZPL: + case PS_PROJ_ZEA: + case PS_PROJ_ZPL: if (R > 2) return NULL; @@ -612,5 +614,5 @@ cosTheta = sqrt (1 - PS_SQR(sinTheta)); break; - default: + default: psAbort("invalid projection"); break; @@ -622,15 +624,15 @@ // Convert from projection spherical coordinates // psLib versions: - # if ELIXIR_CODE +# if ELIXIR_CODE // XXX the elixir version : does the ADD have a sign error? psF64 delta = asin(sinTheta*sinDp - cosTheta*cosPhi*cosDp); psF64 sinAlpha = +cosTheta*sinPhi; psF64 cosAlpha = +cosTheta*cosPhi*sinDp + sinTheta*cosDp; - # else - - psF64 delta = asin(sinTheta*sinDp + cosTheta*cosPhi*cosDp); +# else + + psF64 delta = asin(sinTheta*sinDp + cosTheta*cosPhi*cosDp); psF64 sinAlpha = -cosTheta*sinPhi; psF64 cosAlpha = -cosTheta*cosPhi*sinDp + sinTheta*cosDp; - # endif +# endif out->d = delta; @@ -638,8 +640,8 @@ break; - case PS_PROJECTION_CLASS_PSEUDOCYL: + case PS_PROJECTION_CLASS_PSEUDOCYL: switch (projection->type) { - case PS_PROJ_AIT: + case PS_PROJ_AIT: // Hammer-Aitoff deprojection // XXX EAM : need range check on z^2 : must be > 0 @@ -652,5 +654,5 @@ theta = asin(y*rho); break; - case PS_PROJ_PAR: + case PS_PROJ_PAR: // Parabolic deprojection rho = y/M_PI; @@ -658,9 +660,9 @@ theta = 3.0*asin(rho); break; - case PS_PROJ_GLS: + case PS_PROJ_GLS: phi = x/cos(y); theta = y; break; - default: + default: psAbort("invalid projection"); break; @@ -669,5 +671,5 @@ out->d = theta + projection->D; break; - default: + default: psAbort("invalid projection"); break; @@ -682,5 +684,5 @@ multiplies them. Basically, for each non-zero coeff in the trans1 coeff[][] array, you must multiply by all non-zero coeffs in trans2. - *****************************************************************************/ +*****************************************************************************/ static psPolynomial2D *multiplyDPoly2D( psPolynomial2D *trans1, @@ -713,5 +715,5 @@ /***************************************************************************** psPlaneTransformCombine(out, trans1, trans2) - *****************************************************************************/ +*****************************************************************************/ psPlaneTransform *psPlaneTransformCombine( psPlaneTransform *out, @@ -757,7 +759,7 @@ } else { if ((out->x->nX == orderX) && - (out->x->nY == orderY) && - (out->y->nX == orderX) && - (out->y->nY == orderY)) { + (out->x->nY == orderY) && + (out->y->nX == orderX) && + (out->y->nY == orderY)) { myPT = out; // @@ -885,5 +887,5 @@ XXX: This code ignores nRejIter and sigmaClip. We must call the ClipFit routines instead. - *****************************************************************************/ +*****************************************************************************/ bool psPlaneTransformFit( psPlaneTransform *trans, @@ -932,5 +934,5 @@ psPlaneTransformInvert(out, in, region, nSamples) - *****************************************************************************/ +*****************************************************************************/ psPlaneTransform *psPlaneTransformInvert( psPlaneTransform *out, @@ -972,5 +974,5 @@ } else { if ((out->x->nX == order) && (out->x->nY == order) && - (out->y->nX == order) && (out->y->nX == order)) { + (out->y->nX == order) && (out->y->nX == order)) { myPT = out; } else { @@ -1030,5 +1032,5 @@ const psPlaneTransform *transformation, const psPlane *coord -) + ) { PS_ASSERT_PTR_NON_NULL(transformation, NULL); @@ -1121,5 +1123,5 @@ fxnVal = psPlaneTransformApply(fxnVal, inToOut, coord); if (fabs(fxnVal->x - coord->x) <= fabs(deriv->x) && - fabs(fxnVal->y - coord->y) <= fabs(deriv->y)) { + fabs(fxnVal->y - coord->y) <= fabs(deriv->y)) { int x = (int)(ceil(fabs(deriv->x))); int y = (int)(ceil(fabs(deriv->y))); @@ -1205,42 +1207,45 @@ switch (type) { - case PS_PROJ_LIN: + case PS_PROJ_LIN: psStringAppend (&name, "%s-LIN", prefix); return name; - case PS_PROJ_PLY: + case PS_PROJ_PLY: psStringAppend (&name, "%s-PLY", prefix); return name; - case PS_PROJ_WRP: + case PS_PROJ_WRP: psStringAppend (&name, "%s-WRP", prefix); return name; - case PS_PROJ_TAN: + case PS_PROJ_TAN: psStringAppend (&name, "%s-TAN", prefix); return name; - case PS_PROJ_DIS: + case PS_PROJ_TNX: + psStringAppend (&name, "%s-TNX", prefix); + return name; + case PS_PROJ_DIS: psStringAppend (&name, "%s-DIS", prefix); return name; - case PS_PROJ_SIN: + case PS_PROJ_SIN: psStringAppend (&name, "%s-SIN", prefix); return name; - case PS_PROJ_STG: + case PS_PROJ_STG: psStringAppend (&name, "%s-STG", prefix); return name; - case PS_PROJ_AIT: + case PS_PROJ_AIT: psStringAppend (&name, "%s-AIT", prefix); return name; - case PS_PROJ_PAR: + case PS_PROJ_PAR: psStringAppend (&name, "%s-PAR", prefix); return name; - case PS_PROJ_GLS: + case PS_PROJ_GLS: psStringAppend (&name, "%s-GLS", prefix); return name; - case PS_PROJ_CAR: + case PS_PROJ_CAR: psStringAppend (&name, "%s-CAR", prefix); return name; - case PS_PROJ_MER: + case PS_PROJ_MER: psStringAppend (&name, "%s-MER", prefix); return name; - default: + default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown projection type: %d\n", type); return NULL; @@ -1257,4 +1262,5 @@ if (!strcmp (&name[4], "-WRP")) return PS_PROJ_WRP; if (!strcmp (&name[4], "-TAN")) return PS_PROJ_TAN; + if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX; if (!strcmp (&name[4], "-DIS")) return PS_PROJ_DIS; if (!strcmp (&name[4], "-SIN")) return PS_PROJ_SIN; Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 13838) +++ /trunk/psLib/src/astro/psCoord.h (revision 13839) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-04-04 21:10:03 $ + * @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-15 01:03:22 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -127,4 +127,5 @@ PS_PROJ_SIN, ///< Sine projection PS_PROJ_STG, ///< Sine projection + PS_PROJ_TNX, ///< Sine projection PS_PROJ_ZEA, ///< Sine projection PS_PROJ_ZPL, ///< Sine projection