Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 2634) +++ /trunk/psLib/src/astro/psCoord.c (revision 2635) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-03 19:43:43 $ +* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -19,7 +19,4 @@ /******************************************************************************/ #include "psType.h" -#include "psImage.h" -#include "psArray.h" -#include "psList.h" #include "psCoord.h" #include "psMemory.h" @@ -35,4 +32,10 @@ #define PS_COT(X) (1.0 / atan(X)) +#define DEG_TO_RAD(DEGREES) ((DEGREES) * M_PI / 180.0) +#define MIN_TO_RAD(MINUTES) ((MINUTES) * M_PI / (180.0 * 60.0)) +#define SEC_TO_RAD(SECONDS) ((SECONDS) * M_PI / (180.0 * 60.0 * 60.0)) +#define RAD_TO_DEG(RADIANS) ((RADIANS) * 180 / M_PI) +#define RAD_TO_MIN(RADIANS) ((RADIANS) * 180 * 60.0 / M_PI) +#define RAD_TO_SEC(RADIANS) ((RADIANS) * 180 * 60.0 * 60.0 / M_PI) /******************************************************************************/ @@ -162,9 +165,16 @@ /****************************************************************************** -This function prototype has been modified since the SDRS. - *****************************************************************************/ -psSphereTransform* psSphereTransformAlloc(double alphaP, - double deltaP, - double phiP) +alpha is LONGITUDE +delta is LATITUDE + + alphaP: Take the target pole in the source system; calculate its LONGITUDE + in the target system. That longitude is alphaP. + DeltaP: Take the target pole in the source system; calculate its LATITUDE + in the target system. That longitude is deltaP. + phiP: This is the LONGITUDE of the ascending node in the target system. + *****************************************************************************/ +psSphereTransform* psSphereTransformAlloc(psF64 alphaP, + psF64 deltaP, + psF64 phiP) { psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); @@ -179,10 +189,8 @@ /****************************************************************************** -This algorithm comes from an email from Gene. I assume (x,y) corresponds to -(r,d) in the sphere coordinates. - -XXX: In Gene's email, there are different variables with similar names (y -and Y) and in the code, there's cos_y, cos_Y, and cos(y). Verify that there -are no typo's. +XXX: It's not clear if this code correctly implements what the equations in +the ADD imply. Bug 243 will, hopefully, decide that. + +XXX: Once this code is verified, delete Gene's code. *****************************************************************************/ psSphere* psSphereTransformApply(psSphere* out, @@ -193,32 +201,62 @@ PS_PTR_CHECK_NULL(coord, NULL); - double sinY = 0.0; - double cosY = 0.0; - double sinX = 0.0; - double cosX = 0.0; - double x = 0.0; - double y = 0.0; - double dx = 0.0; - if (out == NULL) { out = (psSphere* ) psAlloc(sizeof(psSphere)); } - x = coord->r; - y = coord->d; - dx = x - transform->phiP; - sinY = cos(y) * sin(dx) * transform->sinDeltaP + - sin(y) * transform->cosDeltaP; - cosY = sqrt(1.0 - sinY * sinY); - sinX = (cos(y) * sin(dx) * transform->cosDeltaP - - sin(y) * transform->sinDeltaP) / cos(y); - cosX = cos(y) * cos(dx) / cos(y); - - out->r = atan2(sinX, cosX) + transform->alphaP; - out->d = atan2(sinY, cosY); - - return (out); -} - + psF64 alpha = coord->r; + psF64 delta = coord->d; + psF64 alphaMinusAlphaP = alpha - transform->alphaP; + + psF64 eq55 = (sin(delta) * transform->cosDeltaP) - + (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP)); + + psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) + + (sin(delta) * transform->sinDeltaP); + + psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP); + + psF64 theta = asin(eq55); + psF64 phi = atan2(eq56, eq57) + transform->alphaP; + + out->r = phi; + out->d = theta; + + return(out); + + /* + This algorithm comes from an email from Gene. I assume (x,y) corresponds to + (r,d) in the sphere coordinates. + + XXX: In Gene's email, there are different variables with similar names (y + and Y) and in the code, there's cos_y, cos_Y, and cos(y). Verify that there + are no typo's. + + psF64 sinY = 0.0; + psF64 cosY = 0.0; + psF64 sinX = 0.0; + psF64 cosX = 0.0; + psF64 x = 0.0; + psF64 y = 0.0; + psF64 dx = 0.0; + + x = coord->r; + y = coord->d; + dx = x - transform->phiP; + sinY = cos(y) * sin(dx) * transform->sinDeltaP + + sin(y) * transform->cosDeltaP; + cosY = sqrt(1.0 - sinY * sinY); + sinX = (cos(y) * sin(dx) * transform->cosDeltaP - + sin(y) * transform->sinDeltaP) / cos(y); + cosX = cos(y) * cos(dx) / cos(y); + + out->r = atan2(sinX, cosX) + transform->alphaP; + out->d = atan2(sinY, cosY); + */ +} + +// XXX: the ADD changed. This code no longer conforms to the ADD. In fact, +// there is no algorithm for this in the ADD. +// XXX: This is bug 244 psSphereTransform* psSphereTransformICRSToEcliptic(psTime time) { @@ -226,9 +264,12 @@ struct tm *tmTime = psTimeToTM(&time); - double year = (double)(1900 + tmTime->tm_year); - double T = year / 100.0; - double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T; - double alphaP = 0.0; - double phiP = 0.0; + psF64 year = (psF64)(1900 + tmTime->tm_year); + psF64 T = year / 100.0; + psF64 phi = DEG_TO_RAD(-23.452294) + + (DEG_TO_RAD(0.013013) * T) + + (DEG_TO_RAD(0.000001639) * T * T) - + (DEG_TO_RAD(0.000000503) * T * T * T); + psF64 alphaP = 0.0; + psF64 phiP = 0.0; psFree(tmTime); @@ -236,4 +277,6 @@ } +// XXX: the ADD changed. This code no longer conforms to the old ADD. +// XXX: Waiting for IfA to sepcify if T should include days/months/etc. psSphereTransform* psSphereTransformEclipticToICRS(psTime time) { @@ -241,24 +284,41 @@ struct tm *tmTime = psTimeToTM(&time); - double year = (double)(1900 + tmTime->tm_year); - double T = year / 100.0; - double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T; - double alphaP = 0.0; - double phiP = 0.0; + psF64 year = (psF64)(1900 + tmTime->tm_year); + psF64 T = year / 100.0; + psF64 alphaP = 0.0; + psF64 deltaP = DEG_TO_RAD(23.0) + + SEC_TO_RAD(27.0) + + MIN_TO_RAD(8.0) - + (SEC_TO_RAD(46.845) * T) - + (SEC_TO_RAD(0.0059) * T * T) + + (SEC_TO_RAD(0.00181) * T * T * T); + psF64 phiP = 0.0; psFree(tmTime); - return (psSphereTransformAlloc(phi, alphaP, phiP)); -} - -// XXX: SHould these be in radians, not degrees? + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); +} + +// XXX: Should these be in radians, not degrees? +// XXX: the ADD changed. This code no longer conforms to the ADD. In fact, +// there is no algorithm for this in the ADD. +// XXX: This is bug 245 psSphereTransform* psSphereTransformICRSToGalatic(void) { - return (psSphereTransformAlloc(62.6, 282.25, 33.0)); -} - -// XXX: SHould these be in radians, not degrees? + psF64 alphaP = DEG_TO_RAD(-62.6); + psF64 deltaP = DEG_TO_RAD(33.0); + psF64 phiP = DEG_TO_RAD(282.25); + + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); +} + +// XXX: Should these be in radians, not degrees? +// XXX: the ADD changed. This code no longer conforms to the old ADD. psSphereTransform* psSphereTransformGalaticToICRS(void) { - return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); + psF64 alphaP = DEG_TO_RAD(282.85948); + psF64 deltaP = DEG_TO_RAD(62.87175); + psF64 phiP = DEG_TO_RAD(32.93192); + + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); } @@ -273,17 +333,15 @@ PS_PTR_CHECK_NULL(projection, NULL); - float R = 0.0; - float alpha = 0.0; + // float R = 0.0; + // float alpha = 0.0; psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane)); if (projection->type == PS_PROJ_TAN) { - R = PS_COT(coord->r) * (180.0 / M_PI); - tmp->x = R * sin(coord->d); - tmp->y = R * cos(coord->d); + tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r); + tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r); } 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); + tmp->x = cos(coord->d) * sin(coord->d); + tmp->y = -cos(coord->d) * cos(coord->d); } else if (projection->type == PS_PROJ_CAR) { @@ -293,12 +351,9 @@ } 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; - + tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r))); } 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); + psF64 tmpF64 = sqrt(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); } else if (projection->type == PS_PROJ_PAR) { @@ -337,7 +392,7 @@ if (projection->type == PS_PROJ_TAN) { + tmp->r = atan2(-coord->y, coord->x); R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = atan2(-coord->y, coord->x); - tmp->r = atan(180.0 / (R * M_PI)); + tmp->d = atan(1.0 / R); } else if (projection->type == PS_PROJ_SIN) { @@ -434,12 +489,12 @@ if (unit == PS_ARCSEC) { - tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI; - tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI; + tmp->r = RAD_TO_SEC(tmp->r); + tmp->d = RAD_TO_SEC(tmp->d); } else if (unit == PS_ARCMIN) { - tmp->r = (tmp->r * 180.0 * 60.0) / M_PI; - tmp->d = (tmp->d * 180.0 * 60.0) / M_PI; + tmp->r = RAD_TO_MIN(tmp->r); + tmp->d = RAD_TO_MIN(tmp->d); } else if (unit == PS_DEGREE) { - tmp->r = (tmp->r * 180.0) / M_PI; - tmp->d = (tmp->d * 180.0) / M_PI; + tmp->r = RAD_TO_DEG(tmp->r); + tmp->d = RAD_TO_DEG(tmp->d); } else if (unit == PS_RADIAN) {} else { @@ -480,6 +535,6 @@ psSphere* tmp; psProjection proj; - double tmpR = 0.0; - double tmpD = 0.0; + psF64 tmpR = 0.0; + psF64 tmpD = 0.0; if (mode == PS_LINEAR) { @@ -495,15 +550,14 @@ tmp = psDeproject(&lin, &proj); return (tmp); - } else if (mode == PS_SPHERICAL) { if (unit == PS_ARCSEC) { - tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0); - tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0); + tmpR = SEC_TO_RAD(offset->r); + tmpD = SEC_TO_RAD(offset->d); } else if (unit == PS_ARCMIN) { - tmpR = (M_PI * offset->r) / (180.0 * 60.0); - tmpD = (M_PI * offset->d) / (180.0 * 60.0); + tmpR = MIN_TO_RAD(offset->r); + tmpD = MIN_TO_RAD(offset->d); } else if (unit == PS_DEGREE) { - tmpR = (M_PI * offset->r) / (180.0); - tmpD = (M_PI * offset->d) / (180.0); + tmpR = DEG_TO_RAD(offset->r); + tmpD = DEG_TO_RAD(offset->d); } else if (unit == PS_RADIAN) { tmpR = offset->r; @@ -533,15 +587,43 @@ } -// XXX: What algorithm should be used? + +/****************************************************************************** +psSpherePrecess(coords, fromTime, toTime): + +XXX: What algorithm should be used? +XXX: Verify that this is correct. +XXX: Use static memory for tmpST. + *****************************************************************************/ psSphere *psSpherePrecess(psSphere *coords, const psTime *fromTime, const psTime *toTime) { - /* - psF64 fromMJD = fromTime->sec/86400.0 + fromTime->usec/86400000000.0 + 40587.0; - psF64 toMJD = toTime->sec/86400.0 + toTime->usec/86400000000.0 + 40587.0; - psF64 = (toMJD - fromMJD) / 36525.0; - */ - return(NULL); -} - + psF64 fromMJD = fromTime->sec/86400.0 + + fromTime->usec/86400000000.0 + + 40587.0; + psF64 toMJD = toTime->sec/86400.0 + + toTime->usec/86400000000.0 + + 40587.0; + psF64 T = (toMJD - fromMJD) / 36525.0; + + psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0000839) * T * T) + + (DEG_TO_RAD(0.000005) * T * T * T)); + + psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) - + (DEG_TO_RAD(0.0001185) * T * T) - + (DEG_TO_RAD(0.0000116) * T * T * T); + + psF64 phiP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0003041) * T * T) + + (DEG_TO_RAD(0.0000051) * T * T * T)); + + psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP); + + psSphere *out = psSphereTransformApply(NULL, tmpST, coords); + + psFree(tmpST); + + return(out); +} + Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 2634) +++ /trunk/psLib/src/astro/psCoord.h (revision 2635) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-02 21:12:51 $ +* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -107,7 +107,7 @@ typedef struct { + double alphaP; ///< Longitude of the target system pole in the source system double cosDeltaP; ///< Cosine of target pole latitude in the source system double sinDeltaP; ///< Sine of target pole latitude in the source system - double alphaP; ///< Longitude of the target system pole in the source system double phiP; ///< Longitude of the ascending node in the target system } Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 2634) +++ /trunk/psLib/src/astronomy/psCoord.c (revision 2635) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-03 19:43:43 $ +* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -19,7 +19,4 @@ /******************************************************************************/ #include "psType.h" -#include "psImage.h" -#include "psArray.h" -#include "psList.h" #include "psCoord.h" #include "psMemory.h" @@ -35,4 +32,10 @@ #define PS_COT(X) (1.0 / atan(X)) +#define DEG_TO_RAD(DEGREES) ((DEGREES) * M_PI / 180.0) +#define MIN_TO_RAD(MINUTES) ((MINUTES) * M_PI / (180.0 * 60.0)) +#define SEC_TO_RAD(SECONDS) ((SECONDS) * M_PI / (180.0 * 60.0 * 60.0)) +#define RAD_TO_DEG(RADIANS) ((RADIANS) * 180 / M_PI) +#define RAD_TO_MIN(RADIANS) ((RADIANS) * 180 * 60.0 / M_PI) +#define RAD_TO_SEC(RADIANS) ((RADIANS) * 180 * 60.0 * 60.0 / M_PI) /******************************************************************************/ @@ -162,9 +165,16 @@ /****************************************************************************** -This function prototype has been modified since the SDRS. - *****************************************************************************/ -psSphereTransform* psSphereTransformAlloc(double alphaP, - double deltaP, - double phiP) +alpha is LONGITUDE +delta is LATITUDE + + alphaP: Take the target pole in the source system; calculate its LONGITUDE + in the target system. That longitude is alphaP. + DeltaP: Take the target pole in the source system; calculate its LATITUDE + in the target system. That longitude is deltaP. + phiP: This is the LONGITUDE of the ascending node in the target system. + *****************************************************************************/ +psSphereTransform* psSphereTransformAlloc(psF64 alphaP, + psF64 deltaP, + psF64 phiP) { psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); @@ -179,10 +189,8 @@ /****************************************************************************** -This algorithm comes from an email from Gene. I assume (x,y) corresponds to -(r,d) in the sphere coordinates. - -XXX: In Gene's email, there are different variables with similar names (y -and Y) and in the code, there's cos_y, cos_Y, and cos(y). Verify that there -are no typo's. +XXX: It's not clear if this code correctly implements what the equations in +the ADD imply. Bug 243 will, hopefully, decide that. + +XXX: Once this code is verified, delete Gene's code. *****************************************************************************/ psSphere* psSphereTransformApply(psSphere* out, @@ -193,32 +201,62 @@ PS_PTR_CHECK_NULL(coord, NULL); - double sinY = 0.0; - double cosY = 0.0; - double sinX = 0.0; - double cosX = 0.0; - double x = 0.0; - double y = 0.0; - double dx = 0.0; - if (out == NULL) { out = (psSphere* ) psAlloc(sizeof(psSphere)); } - x = coord->r; - y = coord->d; - dx = x - transform->phiP; - sinY = cos(y) * sin(dx) * transform->sinDeltaP + - sin(y) * transform->cosDeltaP; - cosY = sqrt(1.0 - sinY * sinY); - sinX = (cos(y) * sin(dx) * transform->cosDeltaP - - sin(y) * transform->sinDeltaP) / cos(y); - cosX = cos(y) * cos(dx) / cos(y); - - out->r = atan2(sinX, cosX) + transform->alphaP; - out->d = atan2(sinY, cosY); - - return (out); -} - + psF64 alpha = coord->r; + psF64 delta = coord->d; + psF64 alphaMinusAlphaP = alpha - transform->alphaP; + + psF64 eq55 = (sin(delta) * transform->cosDeltaP) - + (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP)); + + psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) + + (sin(delta) * transform->sinDeltaP); + + psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP); + + psF64 theta = asin(eq55); + psF64 phi = atan2(eq56, eq57) + transform->alphaP; + + out->r = phi; + out->d = theta; + + return(out); + + /* + This algorithm comes from an email from Gene. I assume (x,y) corresponds to + (r,d) in the sphere coordinates. + + XXX: In Gene's email, there are different variables with similar names (y + and Y) and in the code, there's cos_y, cos_Y, and cos(y). Verify that there + are no typo's. + + psF64 sinY = 0.0; + psF64 cosY = 0.0; + psF64 sinX = 0.0; + psF64 cosX = 0.0; + psF64 x = 0.0; + psF64 y = 0.0; + psF64 dx = 0.0; + + x = coord->r; + y = coord->d; + dx = x - transform->phiP; + sinY = cos(y) * sin(dx) * transform->sinDeltaP + + sin(y) * transform->cosDeltaP; + cosY = sqrt(1.0 - sinY * sinY); + sinX = (cos(y) * sin(dx) * transform->cosDeltaP - + sin(y) * transform->sinDeltaP) / cos(y); + cosX = cos(y) * cos(dx) / cos(y); + + out->r = atan2(sinX, cosX) + transform->alphaP; + out->d = atan2(sinY, cosY); + */ +} + +// XXX: the ADD changed. This code no longer conforms to the ADD. In fact, +// there is no algorithm for this in the ADD. +// XXX: This is bug 244 psSphereTransform* psSphereTransformICRSToEcliptic(psTime time) { @@ -226,9 +264,12 @@ struct tm *tmTime = psTimeToTM(&time); - double year = (double)(1900 + tmTime->tm_year); - double T = year / 100.0; - double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T; - double alphaP = 0.0; - double phiP = 0.0; + psF64 year = (psF64)(1900 + tmTime->tm_year); + psF64 T = year / 100.0; + psF64 phi = DEG_TO_RAD(-23.452294) + + (DEG_TO_RAD(0.013013) * T) + + (DEG_TO_RAD(0.000001639) * T * T) - + (DEG_TO_RAD(0.000000503) * T * T * T); + psF64 alphaP = 0.0; + psF64 phiP = 0.0; psFree(tmTime); @@ -236,4 +277,6 @@ } +// XXX: the ADD changed. This code no longer conforms to the old ADD. +// XXX: Waiting for IfA to sepcify if T should include days/months/etc. psSphereTransform* psSphereTransformEclipticToICRS(psTime time) { @@ -241,24 +284,41 @@ struct tm *tmTime = psTimeToTM(&time); - double year = (double)(1900 + tmTime->tm_year); - double T = year / 100.0; - double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T; - double alphaP = 0.0; - double phiP = 0.0; + psF64 year = (psF64)(1900 + tmTime->tm_year); + psF64 T = year / 100.0; + psF64 alphaP = 0.0; + psF64 deltaP = DEG_TO_RAD(23.0) + + SEC_TO_RAD(27.0) + + MIN_TO_RAD(8.0) - + (SEC_TO_RAD(46.845) * T) - + (SEC_TO_RAD(0.0059) * T * T) + + (SEC_TO_RAD(0.00181) * T * T * T); + psF64 phiP = 0.0; psFree(tmTime); - return (psSphereTransformAlloc(phi, alphaP, phiP)); -} - -// XXX: SHould these be in radians, not degrees? + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); +} + +// XXX: Should these be in radians, not degrees? +// XXX: the ADD changed. This code no longer conforms to the ADD. In fact, +// there is no algorithm for this in the ADD. +// XXX: This is bug 245 psSphereTransform* psSphereTransformICRSToGalatic(void) { - return (psSphereTransformAlloc(62.6, 282.25, 33.0)); -} - -// XXX: SHould these be in radians, not degrees? + psF64 alphaP = DEG_TO_RAD(-62.6); + psF64 deltaP = DEG_TO_RAD(33.0); + psF64 phiP = DEG_TO_RAD(282.25); + + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); +} + +// XXX: Should these be in radians, not degrees? +// XXX: the ADD changed. This code no longer conforms to the old ADD. psSphereTransform* psSphereTransformGalaticToICRS(void) { - return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); + psF64 alphaP = DEG_TO_RAD(282.85948); + psF64 deltaP = DEG_TO_RAD(62.87175); + psF64 phiP = DEG_TO_RAD(32.93192); + + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); } @@ -273,17 +333,15 @@ PS_PTR_CHECK_NULL(projection, NULL); - float R = 0.0; - float alpha = 0.0; + // float R = 0.0; + // float alpha = 0.0; psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane)); if (projection->type == PS_PROJ_TAN) { - R = PS_COT(coord->r) * (180.0 / M_PI); - tmp->x = R * sin(coord->d); - tmp->y = R * cos(coord->d); + tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r); + tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r); } 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); + tmp->x = cos(coord->d) * sin(coord->d); + tmp->y = -cos(coord->d) * cos(coord->d); } else if (projection->type == PS_PROJ_CAR) { @@ -293,12 +351,9 @@ } 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; - + tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r))); } 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); + psF64 tmpF64 = sqrt(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); } else if (projection->type == PS_PROJ_PAR) { @@ -337,7 +392,7 @@ if (projection->type == PS_PROJ_TAN) { + tmp->r = atan2(-coord->y, coord->x); R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = atan2(-coord->y, coord->x); - tmp->r = atan(180.0 / (R * M_PI)); + tmp->d = atan(1.0 / R); } else if (projection->type == PS_PROJ_SIN) { @@ -434,12 +489,12 @@ if (unit == PS_ARCSEC) { - tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI; - tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI; + tmp->r = RAD_TO_SEC(tmp->r); + tmp->d = RAD_TO_SEC(tmp->d); } else if (unit == PS_ARCMIN) { - tmp->r = (tmp->r * 180.0 * 60.0) / M_PI; - tmp->d = (tmp->d * 180.0 * 60.0) / M_PI; + tmp->r = RAD_TO_MIN(tmp->r); + tmp->d = RAD_TO_MIN(tmp->d); } else if (unit == PS_DEGREE) { - tmp->r = (tmp->r * 180.0) / M_PI; - tmp->d = (tmp->d * 180.0) / M_PI; + tmp->r = RAD_TO_DEG(tmp->r); + tmp->d = RAD_TO_DEG(tmp->d); } else if (unit == PS_RADIAN) {} else { @@ -480,6 +535,6 @@ psSphere* tmp; psProjection proj; - double tmpR = 0.0; - double tmpD = 0.0; + psF64 tmpR = 0.0; + psF64 tmpD = 0.0; if (mode == PS_LINEAR) { @@ -495,15 +550,14 @@ tmp = psDeproject(&lin, &proj); return (tmp); - } else if (mode == PS_SPHERICAL) { if (unit == PS_ARCSEC) { - tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0); - tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0); + tmpR = SEC_TO_RAD(offset->r); + tmpD = SEC_TO_RAD(offset->d); } else if (unit == PS_ARCMIN) { - tmpR = (M_PI * offset->r) / (180.0 * 60.0); - tmpD = (M_PI * offset->d) / (180.0 * 60.0); + tmpR = MIN_TO_RAD(offset->r); + tmpD = MIN_TO_RAD(offset->d); } else if (unit == PS_DEGREE) { - tmpR = (M_PI * offset->r) / (180.0); - tmpD = (M_PI * offset->d) / (180.0); + tmpR = DEG_TO_RAD(offset->r); + tmpD = DEG_TO_RAD(offset->d); } else if (unit == PS_RADIAN) { tmpR = offset->r; @@ -533,15 +587,43 @@ } -// XXX: What algorithm should be used? + +/****************************************************************************** +psSpherePrecess(coords, fromTime, toTime): + +XXX: What algorithm should be used? +XXX: Verify that this is correct. +XXX: Use static memory for tmpST. + *****************************************************************************/ psSphere *psSpherePrecess(psSphere *coords, const psTime *fromTime, const psTime *toTime) { - /* - psF64 fromMJD = fromTime->sec/86400.0 + fromTime->usec/86400000000.0 + 40587.0; - psF64 toMJD = toTime->sec/86400.0 + toTime->usec/86400000000.0 + 40587.0; - psF64 = (toMJD - fromMJD) / 36525.0; - */ - return(NULL); -} - + psF64 fromMJD = fromTime->sec/86400.0 + + fromTime->usec/86400000000.0 + + 40587.0; + psF64 toMJD = toTime->sec/86400.0 + + toTime->usec/86400000000.0 + + 40587.0; + psF64 T = (toMJD - fromMJD) / 36525.0; + + psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0000839) * T * T) + + (DEG_TO_RAD(0.000005) * T * T * T)); + + psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) - + (DEG_TO_RAD(0.0001185) * T * T) - + (DEG_TO_RAD(0.0000116) * T * T * T); + + psF64 phiP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + (DEG_TO_RAD(0.0003041) * T * T) + + (DEG_TO_RAD(0.0000051) * T * T * T)); + + psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP); + + psSphere *out = psSphereTransformApply(NULL, tmpST, coords); + + psFree(tmpST); + + return(out); +} + Index: /trunk/psLib/src/astronomy/psCoord.h =================================================================== --- /trunk/psLib/src/astronomy/psCoord.h (revision 2634) +++ /trunk/psLib/src/astronomy/psCoord.h (revision 2635) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-02 21:12:51 $ +* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -107,7 +107,7 @@ typedef struct { + double alphaP; ///< Longitude of the target system pole in the source system double cosDeltaP; ///< Cosine of target pole latitude in the source system double sinDeltaP; ///< Sine of target pole latitude in the source system - double alphaP; ///< Longitude of the target system pole in the source system double phiP; ///< Longitude of the ascending node in the target system } Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 2634) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 2635) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-01 19:56:06 $ + * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,6 +24,4 @@ *****************************************************************************/ #define PS_ONE 1.0 - - /***************************************************************************** Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 2634) +++ /trunk/psLib/src/math/psConstants.h (revision 2635) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-01 19:56:06 $ + * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-06 19:49:36 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -24,6 +24,4 @@ *****************************************************************************/ #define PS_ONE 1.0 - - /*****************************************************************************