Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 1287) +++ trunk/psLib/src/astro/psCoord.c (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -22,5 +22,7 @@ #include "psCoord.h" #include "psMemory.h" - +#include "psAbort.h" +#include +#include psPlane *psPlaneTransformApply(psPlane *out, @@ -127,4 +129,117 @@ } +// XXX: Is this the correct way to calculate this? +float cot(float x) +{ + return(1.0 / atan(x)); +} + +float arg(float x, float y) +{ + if (x > 0) { + return(atan(y/x)); + } else if ((x==0) && (y==0)) { + return(0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(-0.5 * M_PI); + } else if ((x==0) && (y==0)) { + return(M_PI + atan(y/x)); + } else if ((x==0) && (y==0)) { + return(-M_PI + atan(y/x)); + } + + psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y); + return(0.0); +} + + +psPlane *psProject(const psSphere *coord, + const psProjection *projection) +{ + float R = 0.0; + float alpha = 0.0; + psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane)); + + if (projection->type == PS_PROJ_TAN) { + R = cot(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_SIN) { + R = cos(coord->r) * (180.0 / M_PI); + tmp->x = R * sin(coord->d); + tmp->y = R * cos(coord->d); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->x = coord->d; + tmp->y = coord->r; + + } else if (projection->type == PS_PROJ_MER) { + tmp->x = coord->d; + tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI; + + } else if (projection->type == PS_PROJ_AIT) { + alpha = 1.0 / ((180.0 / M_PI) * + sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5))); + + tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d); + tmp->y = alpha * sin(coord->d); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} + +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection) +{ + float R = 0.0; + float chu = 0.0; + float chu1 = 0.0; + float chu2 = 0.0; + psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere)); + + if (projection->type == PS_PROJ_TAN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = atan(180.0 / (R * M_PI)); + + } else if (projection->type == PS_PROJ_SIN) { + R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); + tmp->d = arg(-coord->y, coord->x); + tmp->r = acos((R * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_CAR) { + tmp->d = coord->x; + tmp->r = coord->y; + + } else if (projection->type == PS_PROJ_MER) { + tmp->d = coord->x; + tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0; + + } else if (projection->type == PS_PROJ_AIT) { + chu1 = (coord->x * M_PI) / 720.0; + chu1*= chu1; + chu2 = (coord->y * M_PI) / 360.0; + chu2*= chu2; + chu = sqrt(1.0 - chu1 - chu2); + tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, + (coord->x * chu * M_PI)/360.0); + tmp->r = asin((coord->y * chu * M_PI) / 180.0); + + } else if (projection->type == PS_PROJ_PAR) { + psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n"); + + } else if (projection->type == PS_PROJ_GLS) { + psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n"); + } + + return(tmp); +} psSphere *psSphereGetOffset(const psSphere *restrict position1, @@ -133,8 +248,11 @@ psSphereOffsetUnit unit) { - return(NULL); -} - - + psAbort(__func__, "This function has not be dedfined.\n"); + return(NULL); +} + + +// XXX: I copied the algorithm from the ADD exactly. Arguments mode and unit +// are ignored. psSphere *psSphereSetOffset(const psSphere *restrict position, const psSphere *restrict offset, @@ -142,4 +260,24 @@ psSphereOffsetUnit unit) { - return(NULL); -} + psPlane lin; + psSphere *tmp; + psProjection proj; + + if (mode == PS_LINEAR) { + proj.R = position->r; + proj.D = position->d; + proj.Xs = 0.0; + proj.Ys = 0.0; + proj.type = PS_PROJ_TAN; + + lin.x = offset->r; + lin.y = offset->d; + + tmp = psDeproject(&lin, &proj); + return(tmp); + } else if (mode == PS_SPHERICAL) { + psAbort(__func__, "Sperical offset modes are not defined.\n"); + } + psAbort(__func__, "Unrecognized offset mode\n"); + return(NULL); +} Index: trunk/psLib/src/astro/psCoord.h =================================================================== --- trunk/psLib/src/astro/psCoord.h (revision 1287) +++ trunk/psLib/src/astro/psCoord.h (revision 1293) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-23 03:13:39 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-24 02:42:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,6 +27,6 @@ typedef struct { - double x; ///< x position - double y; ///< y position + double x; ///< x position + double y; ///< y position double xErr; ///< Error in x position double yErr; ///< Error in y position @@ -36,6 +36,6 @@ typedef struct { - double r; ///< RA - double d; ///< Dec + double r; ///< RA + double d; ///< Dec double rErr; ///< Error in RA double dErr; ///< Error in Dec @@ -63,22 +63,24 @@ double sinNPlat; ///< sin of North Pole latitude double cosNPlat; ///< cos of North Pole latitude - double sinZP; ///< sin of Forst PT os Ares lon - double cosZP; ///< cos of Forst PT os Ares lon + double sinZP; ///< sin of Forst PT os Ares lon + double cosZP; ///< cos of Forst PT os Ares lon } psSphereTransform; typedef enum { - PS_PROJ_TAN, ///< Tangent projection - PS_PROJ_SIN, ///< Sine projection - PS_PROJ_AIT, ///< Aitoff projection - PS_PROJ_PAR, ///< Par projection - PS_PROJ_GLS, ///< GLS projection - PS_PROJ_NTYPE ///< Number of types; must be last. + PS_PROJ_TAN, ///< Tangent projection + PS_PROJ_SIN, ///< 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_NTYPE ///< Number of types; must be last. } psProjectionType; typedef struct { - double R; ///< Coordinates of projection center - double D; ///< Coordinates of projection center + double R; ///< Coordinates of projection center + double D; ///< Coordinates of projection center double Xs; ///< plate-scale in X direction double Ys; ///< plate-scale in Y direction @@ -128,6 +130,6 @@ const psProjection *projection); -psSphere *psProject(const psPlane *coord, - const psProjection *projection); +psSphere *psDeproject(const psPlane *coord, + const psProjection *projection); psSphere *psSphereGetOffset(const psSphere *restrict position1,