Index: trunk/psLib/src/astronomy/psCoord.c =================================================================== --- trunk/psLib/src/astronomy/psCoord.c (revision 3359) +++ trunk/psLib/src/astronomy/psCoord.c (revision 3450) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.57 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-03-02 03:39:58 $ +* @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-03-18 21:34:43 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -399,6 +399,4 @@ } -// XXX: Must test psProjectionAlloc() and projectionFree(). -// XXX: Must rewrite code and tests to use these functions. psProjection* psProjectionAlloc( psF64 R, @@ -419,10 +417,4 @@ } - -/****************************************************************************** -XXX: Waiting for the definition of the PS_PROJ_PAR projection. -XXX: Waiting for the definition of the PS_PROJ_GLS projection. -XXX: Must apply scaling at the end. - *****************************************************************************/ psPlane* psProject(const psSphere* coord, const psProjection* projection) @@ -431,83 +423,50 @@ PS_PTR_CHECK_NULL(projection, NULL); - psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane)); + psF64 theta = 0.0; + psF64 phi = 0.0; + + // Allocate return value + psPlane* out = psPlaneAlloc(); + + // Convert to projection spherical coordinate system + theta = asin( sin(coord->d)*sin(projection->D) + + cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R)); + phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R), + sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) ); + + // Perform the specified projection + // Gnomonic projection if (projection->type == PS_PROJ_TAN) { - - /* - // ******************************************** - // From the ADD - // ******************************************** - // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); - // return(NULL); - psF32 sinTheta; - psF32 cosThetaCosPhi; - psF32 cosThetaSinPhi; - - sinTheta = (sin(coord->d) * sin(projection->D)) + - (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R))); - cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) - - (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R))); - cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R); - - tmp->x = -cosThetaSinPhi / sinTheta; - tmp->y = cosThetaCosPhi / sinTheta; - - */ - // ******************************************** - // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html - // delta_0 and phi_1 are the projection centers, not the point being projected. - // (delta_0, phi_1) == (projection->R, projection->D) - // ******************************************** - psF32 cosC = (sin(projection->D) * sin(coord->d)) + - (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R)); - tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC; - tmp->y = ((cos(projection->D) * sin(coord->d)) - - (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) / - cosC; - + out->x = (cos(theta)*sin(phi))/sin(theta); + out->y = (-1.0*cos(theta)*cos(phi))/sin(theta); + // Othrographic projection } else if (projection->type == PS_PROJ_SIN) { - // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work."); - // return(NULL); - - tmp->x = cos(coord->d) * sin(coord->d); - tmp->y = -cos(coord->d) * 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(DEG_TO_RAD(45.0) + (0.5 * coord->r))); - } else if (projection->type == PS_PROJ_AIT) { - psF64 tmpF64 = PS_SQRT_F32(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d)))); - tmpF64 = 1.0 / tmpF64; - tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d); - tmp->y = sin(coord->r) * tmpF64; - - } else if (projection->type == PS_PROJ_PAR) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED, - "PS_PROJ_PAR"); - - } else if (projection->type == PS_PROJ_GLS) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED, - "PS_PROJ_GLS"); + out->x = cos(theta)*sin(phi); + out->y = -1.0*cos(theta)*cos(phi); + // Hammer-Aitoff projection + } else if ( projection->type == PS_PROJ_AIT) { + psF64 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); + // Parabolic projection + } else if ( projection->type == PS_PROJ_PAR) { + out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); + out->y = PS_PI*sin(theta/3.0); } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, projection->type); - } - - return (tmp); -} - -/****************************************************************************** -XXX: Waiting for algorithms. -XXX: Waiting for the definition of the PS_PROJ_PAR projection. -XXX: Waiting for the definition of the PS_PROJ_GLS projection. -XXX: Must apply scaling at the beginning. - *****************************************************************************/ + psFree(out); + return NULL; + } + + // Apply plate scales + out->x *= projection->Xs; + out->y *= projection->Ys; + + // Return output + return out; +} + psSphere* psDeproject(const psPlane* coord, const psProjection* projection) @@ -515,90 +474,52 @@ PS_PTR_CHECK_NULL(coord, NULL); PS_PTR_CHECK_NULL(projection, NULL); - // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work."); - // return(NULL); - 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) { - /* - // ******************************************** - // From the ADD - // ******************************************** - psF32 phi; - psF32 theta; - psF32 sinDelta; - psF32 cosDeltaCosAlphaMinusAlphaP; - psF32 cosDeltaSinAlphaMinusAlphaP; - - phi = atan2(-coord->y, coord->x); - R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); - theta = atan(1.0 / R); - - sinDelta = (sin(theta) * sin(projection->D)) + - (cos(theta) * cos(projection->D) * cos(phi)); - cosDeltaCosAlphaMinusAlphaP = - (sin(theta) * cos(projection->D)) - - (cos(theta) * sin(projection->D) * cos(phi)); - cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi); - tmp->d = asin(sinDelta); - tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R; - */ - // ******************************************** - // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html - // delta_0 and phi_1 are the projection centers, not the point being projected. - // (delta_0, phi_1) == (projection->R, projection->D) - // XXX: figure out how to use the two-argument atan2() here. - // ******************************************** - psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); - psF32 C = atan(row); - tmp->d = asin((cos(C) * sin(projection->D)) + - ((coord->y * sin(C) * cos(projection->D)) / row)); - psF32 tmpAtan = atan((coord->x * sin(C)) / - ((row * cos(projection->D) * cos(C)) - - (coord->y * sin(projection->D) * sin(C)))); - tmp->r = projection->R + tmpAtan; - - - } else if (projection->type == PS_PROJ_SIN) { - R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = atan2(-coord->y, coord->x); - tmp->r = acos((R * PS_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 * PS_PI / 180.0)))) - 180.0; - - } else if (projection->type == PS_PROJ_AIT) { - chu1 = (coord->x * PS_PI) / 720.0; - chu1 *= chu1; - chu2 = (coord->y * PS_PI) / 360.0; - chu2 *= chu2; - chu = PS_SQRT_F32(1.0 - chu1 - chu2); - tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * PS_PI) / 360.0); - tmp->r = asin((coord->y * chu * PS_PI) / 180.0); - - } else if (projection->type == PS_PROJ_PAR) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED, - "PS_PROJ_PAR"); - - } else if (projection->type == PS_PROJ_GLS) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED, - "PS_PROJ_GLS"); + + psF64 theta = 0.0; + psF64 phi = 0.0; + + // Allocate return sphere structure + psSphere* out = psSphereAlloc(); + + // Remove plate scales + psF64 x = coord->x/projection->Xs; + psF64 y = coord->y/projection->Ys; + + // Perform inverse projection + // Gnonomic deprojection + if ( projection->type == PS_PROJ_TAN) { + phi = atan(-1.0*x/y); + theta = atan(1.0/sqrt(x*x+y*y)); + // Orhtographic deprojection + } else if ( projection->type == PS_PROJ_SIN) { + phi = atan((-1.0*x)/y); + theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y)); + // Hammer-Aitoff deprojection + } else if ( projection->type == PS_PROJ_AIT) { + psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0))); + phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) ); + theta = asin(y*z); + // Parabolic deprojection + } else if ( projection->type == PS_PROJ_PAR) { + psF64 rho = y/PS_PI; + phi = x/(1.0 - 4.0*rho*rho); + theta = 3.0*asin(rho); + // Invalid deprojection type } else { psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, projection->type); - } - - return (tmp); + psFree(out); + return NULL; + } + + // Convert from projection spherical coordinates + out->d = asin( sin(theta)*sin(projection->D) + + cos(theta)*cos(projection->D)*cos(phi) ); + out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi), + sin(theta)*cos(projection->D) - + cos(theta)*sin(projection->D)*cos(phi) ); + + // Return sphere coordinate + return out; }