Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 10822) +++ trunk/psLib/src/astro/psCoord.c (revision 10823) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.127 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-10 17:29:29 $ +* @version $Revision: 1.128 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-22 21:19:47 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -26,4 +26,5 @@ #include "psLogMsg.h" #include "psTrace.h" +#include "psAbort.h" #include "psMatrix.h" @@ -347,5 +348,10 @@ } - +typedef enum { + PS_PROJ_NONE, + PS_PROJ_CARTESIAN, + PS_PROJ_ZENITHAL, + PS_PROJ_PSEUDOCYL, +}psProjectionClass; psPlane* psProject( @@ -364,5 +370,31 @@ psF64 sinDp, cosDp, sinAlpha, cosAlpha, sinDelta, cosDelta, zeta; - bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); + psProjectionClass class = PS_PROJ_NONE; + switch (projection->type) { + case PS_PROJ_LIN: + case PS_PROJ_PLY: + case PS_PROJ_WRP: + class = PS_PROJ_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: + class = PS_PROJ_ZENITHAL; + break; + case PS_PROJ_AIT: + case PS_PROJ_PAR: + case PS_PROJ_GLS: + case PS_PROJ_MER: + case PS_PROJ_CAR: + class = PS_PROJ_PSEUDOCYL; + break; + + default: + psAbort ("psCoord", "invalid projection"); + break; + } // Allocate return value @@ -374,5 +406,11 @@ } - if (zenithal) { + switch (class) { + case PS_PROJ_CARTESIAN: + out->x = (coord->r - projection->R); + out->y = (coord->d - projection->D); + break; + + case PS_PROJ_ZENITHAL: sinDp = sin(projection->D); cosDp = cos(projection->D); @@ -389,42 +427,74 @@ # else - sinTheta = sinDelta*sinDp + cosDelta*cosDp*cosAlpha; - cosPhiCT = sinDelta*cosDp - cosDelta*sinDp*cosAlpha; + sinTheta = cosDelta*cosAlpha*cosDp + sinDelta*sinDp; sinPhiCT = -cosDelta*sinAlpha; + cosPhiCT = -cosDelta*cosAlpha*sinDp + sinDelta*cosDp; # endif - -} else { - phi = coord->r - projection->R; - theta = coord->d - projection->D; -} - -// Perform the specified projection -switch (projection->type) { - case PS_PROJ_TAN: - // 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_AIT: - // Hammer-Aitoff projection - zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); - out->x = 2.0*zeta*cos(theta)*sin(phi/2.0); - out->y = zeta*sin(theta); - break; - 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); + // 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 ("psCoord", "invalid projection"); + break; + } + break; + + case PS_PROJ_PSEUDOCYL: + phi = coord->r - projection->R; + theta = coord->d - projection->D; + switch (projection->type) + { + case PS_PROJ_AIT: + // Hammer-Aitoff projection + zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); + out->x = 2.0*zeta*cos(theta)*sin(phi/2.0); + out->y = zeta*sin(theta); + break; + case PS_PROJ_GLS: + // projection name? + out->x = phi * cos(theta); + out->x = theta; + break; + 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: + psAbort ("psCoord", "projection not yet implemented"); + break; + default: + psAbort ("psCoord", "invalid projection"); + break; + } + break; + default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("The projection type, %d, is unknown."), - projection->type); - psFree(out); - return NULL; + psAbort ("psCoord", "invalid projection"); + break; } @@ -446,4 +516,5 @@ psF64 rho = 0.0; + psF64 rho2 = 0.0; psF64 sinTheta = 0.0; psF64 cosTheta = 0.0; @@ -453,4 +524,33 @@ psF64 theta = 0.0; psF64 phi = 0.0; + psF64 x, y, R; + + psProjectionClass class = PS_PROJ_NONE; + switch (projection->type) { + case PS_PROJ_LIN: + case PS_PROJ_PLY: + case PS_PROJ_WRP: + class = PS_PROJ_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: + class = PS_PROJ_ZENITHAL; + break; + case PS_PROJ_AIT: + case PS_PROJ_PAR: + case PS_PROJ_GLS: + case PS_PROJ_CAR: + case PS_PROJ_MER: + class = PS_PROJ_PSEUDOCYL; + break; + + default: + psAbort ("psCoord", "invalid projection"); + break; + } // Allocate return sphere structure @@ -466,63 +566,48 @@ } - // Remove plate scales - psF64 x = coord->x*projection->Xs; - psF64 y = coord->y*projection->Ys; - psF64 R = sqrt(x*x + y*y); - - bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); - // Perform inverse projection - switch (projection->type) { - case PS_PROJ_TAN: - // Gnonomic deprojection - // the MHPCC version here was fine - # if 0 - + switch (class) { + case PS_PROJ_CARTESIAN: + out->r = coord->x*projection->Xs + projection->R; + out->d = coord->y*projection->Ys + projection->D; + break; + + case PS_PROJ_ZENITHAL: + // Remove plate scales + x = coord->x*projection->Xs; + y = coord->y*projection->Ys; + R = sqrt(x*x + y*y); sinPhi = (R == 0) ? 0.0 : +x / R; cosPhi = (R == 0) ? 1.0 : -y / R; - rho = sqrt (1 + R*R); - sinTheta = 1 / rho; - cosTheta = R / rho; - # else - + + switch (projection->type) + { + case PS_PROJ_TAN: + case PS_PROJ_DIS: + // Gnonomic deprojection rho = sqrt (1 + R*R); - sinTheta = 1 / rho; - cosTheta = R / rho; - sinPhi = (R == 0) ? 0.0 : +x / R; - cosPhi = (R == 0) ? 1.0 : -y / R; - # endif - - break; - case PS_PROJ_SIN: - // Orhtographic deprojection - sinTheta = sqrt (1 - R*R); - cosTheta = R; - sinPhi = (R == 0) ? 0.0 : +x / R; - cosPhi = (R == 0) ? 1.0 : -y / R; - break; - case PS_PROJ_AIT: - // Hammer-Aitoff deprojection - // XXX EAM : need range check on z^2 : must be > 0 - // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2 - rho = sqrt(1.0 - PS_SQR(x/4.0) - PS_SQR(y/2.0)); - phi = 2.0*atan2((2.0*rho*rho-1.0), x*rho); - theta = asin(y*rho); - break; - case PS_PROJ_PAR: - // Parabolic deprojection - rho = y/M_PI; - phi = x/(1.0 - 4.0*rho*rho); - theta = 3.0*asin(rho); - break; - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("The projection type, %d, is unknown."), - projection->type); - psFree(out); - return NULL; - } - - if (zenithal) { + sinTheta = 1 / rho; + cosTheta = R / rho; + break; + case PS_PROJ_SIN: + // Orhtographic deprojection + sinTheta = sqrt (1 - R*R); + cosTheta = R; + break; + case PS_PROJ_STG: + psAbort ("coords", "STG not defined"); + break; + case PS_PROJ_ZEA: + case PS_PROJ_ZPL: + if (R > 2) + return NULL; + sinTheta = 1 - 0.5*PS_SQR(R); + cosTheta = sqrt (1 - PS_SQR(sinTheta)); + break; + default: + psAbort ("psCoord", "invalid projection"); + break; + } + psF64 sinDp = sin(projection->D); psF64 cosDp = cos(projection->D); @@ -531,11 +616,11 @@ // psLib versions: # 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 = sinTheta*cosDp + cosTheta*cosPhi*sinDp; + // 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*cosDp*cosPhi); + psF64 delta = asin(sinTheta*sinDp + cosTheta*cosPhi*cosDp); psF64 sinAlpha = -cosTheta*sinPhi; psF64 cosAlpha = -cosTheta*cosPhi*sinDp + sinTheta*cosDp; @@ -544,7 +629,40 @@ out->d = delta; out->r = atan2(sinAlpha, cosAlpha) + projection->R; - } else { + break; + + case PS_PROJ_PSEUDOCYL: + switch (projection->type) + { + case PS_PROJ_AIT: + // Hammer-Aitoff deprojection + // XXX EAM : need range check on z^2 : must be > 0 + // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2 + rho2 = 1.0 - PS_SQR(0.25*x) - PS_SQR(0.5*y); + if (rho2 < 0) + return (NULL); + rho = sqrt(rho2); + phi = 2.0*atan2(0.5*x*rho, 2.0*rho2 - 1.0); + theta = asin(y*rho); + break; + case PS_PROJ_PAR: + // Parabolic deprojection + rho = y/M_PI; + phi = x/(1.0 - 4.0*rho*rho); + theta = 3.0*asin(rho); + break; + case PS_PROJ_GLS: + phi = x/cos(y); + theta = y; + break; + default: + psAbort ("psCoord", "invalid projection"); + break; + } out->r = phi + projection->R; out->d = theta + projection->D; + break; + default: + psAbort ("psCoord", "invalid projection"); + break; } @@ -1070,2 +1188,83 @@ return(sphere); } + +char * psProjectTypeToString (psProjectionType type, char *prefix) +{ + + char *name = NULL; + + switch (type) { + case PS_PROJ_LIN: + psStringAppend (&name, "%s-LIN", prefix); + return name; + case PS_PROJ_PLY: + psStringAppend (&name, "%s-PLY", prefix); + return name; + case PS_PROJ_WRP: + psStringAppend (&name, "%s-WRP", prefix); + return name; + case PS_PROJ_TAN: + psStringAppend (&name, "%s-TAN", prefix); + return name; + case PS_PROJ_DIS: + psStringAppend (&name, "%s-DIS", prefix); + return name; + case PS_PROJ_SIN: + psStringAppend (&name, "%s-SIN", prefix); + return name; + case PS_PROJ_STG: + psStringAppend (&name, "%s-STG", prefix); + return name; + case PS_PROJ_AIT: + psStringAppend (&name, "%s-AIT", prefix); + return name; + case PS_PROJ_PAR: + psStringAppend (&name, "%s-PAR", prefix); + return name; + case PS_PROJ_GLS: + psStringAppend (&name, "%s-GLS", prefix); + return name; + case PS_PROJ_CAR: + psStringAppend (&name, "%s-CAR", prefix); + return name; + case PS_PROJ_MER: + psStringAppend (&name, "%s-MER", prefix); + return name; + + default: + psLogMsg ("psLib.astrom", 2, "warning: unknown projection type %d\n", type); + return NULL; + } + return NULL; +} + +psProjectionType psProjectTypeFromString (char *name) +{ + + psProjectionType type = PS_PROJ_NTYPE; + if (!strcmp (&name[4], "-LIN")) + type = PS_PROJ_LIN; + if (!strcmp (&name[4], "-PLY")) + type = PS_PROJ_PLY; + if (!strcmp (&name[4], "-WRP")) + type = PS_PROJ_WRP; + if (!strcmp (&name[4], "-TAN")) + type = PS_PROJ_TAN; + if (!strcmp (&name[4], "-DIS")) + type = PS_PROJ_DIS; + if (!strcmp (&name[4], "-SIN")) + type = PS_PROJ_SIN; + if (!strcmp (&name[4], "-STG")) + type = PS_PROJ_STG; + if (!strcmp (&name[4], "-AIT")) + type = PS_PROJ_AIT; + if (!strcmp (&name[4], "-PAR")) + type = PS_PROJ_PAR; + if (!strcmp (&name[4], "-GLS")) + type = PS_PROJ_GLS; + if (!strcmp (&name[4], "-CAR")) + type = PS_PROJ_CAR; + if (!strcmp (&name[4], "-MER")) + type = PS_PROJ_MER; + return type; +} Index: trunk/psLib/src/astro/psCoord.h =================================================================== --- trunk/psLib/src/astro/psCoord.h (revision 10822) +++ trunk/psLib/src/astro/psCoord.h (revision 10823) @@ -11,6 +11,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.54 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-10 17:29:51 $ + * @version $Revision: 1.55 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-22 21:19:47 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -122,11 +122,18 @@ */ typedef enum { + PS_PROJ_LIN, ///< Linear projection + PS_PROJ_PLY, ///< Linear polynomial projection + PS_PROJ_WRP, ///< Linear polynomial projection PS_PROJ_TAN, ///< Tangent projection + PS_PROJ_DIS, ///< Sine projection PS_PROJ_SIN, ///< Sine projection + PS_PROJ_STG, ///< Sine projection + PS_PROJ_ZEA, ///< Sine projection + PS_PROJ_ZPL, ///< Sine projection PS_PROJ_AIT, ///< Aitoff projection PS_PROJ_PAR, ///< Par projection - // PS_PROJ_GLS, ///< GLS projection - // PS_PROJ_CAR, ///< CAR projection - // PS_PROJ_MER, ///< MER projection + PS_PROJ_GLS, ///< GLS projection + PS_PROJ_CAR, ///< CAR projection + PS_PROJ_MER, ///< MER projection PS_PROJ_NTYPE ///< Number of types; must be last. } psProjectionType; @@ -438,4 +445,8 @@ ); +char * psProjectTypeToString (psProjectionType type, char *prefix); +psProjectionType psProjectTypeFromString (char *name); + + #define PS_PRINT_PLANE_TRANSFORM(NAME) \ { \