Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 1531) +++ /trunk/psLib/src/astro/psCoord.c (revision 1532) @@ -10,10 +10,12 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-08-13 23:33:13 $ +* @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-08-13 23:43:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ - +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ #include "psType.h" #include "psImage.h" @@ -26,141 +28,39 @@ #include #include - -static float cot(float x); -static float arg(float x, float y); - -psPlane* psPlaneTransformApply(psPlane* out, - const psPlaneTransform* transform, - const psPlane* coords) -{ - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - out->x = transform->x->coeff[0][0] + - (transform->x->coeff[1][0] * coords->x) + - (transform->x->coeff[0][1] * coords->y); - - out->y = transform->y->coeff[0][0] + - (transform->y->coeff[1][0] * coords->x) + - (transform->y->coeff[0][1] * coords->y); - - return (out); -} - -// This transformation takes into account parameters beyond an objects -// spatial coordinates: term3 and term4 (magnitude and color). -psPlane* psPlaneDistortApply(psPlane* out, - const psPlaneDistort* transform, - const psPlane* coords, - float term3, - float term4) -{ - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - - out->x = transform->x->coeff[0][0][0][0] + - (transform->x->coeff[1][0][0][0] * coords->x) + - (transform->x->coeff[0][1][0][0] * coords->y) + - (transform->x->coeff[0][0][1][0] * term3) + - (transform->x->coeff[0][0][0][1] * term4); - - out->y = transform->y->coeff[0][0][0][0] + - (transform->y->coeff[1][0][0][0] * coords->x) + - (transform->y->coeff[0][1][0][0] * coords->y) + - (transform->y->coeff[0][0][1][0] * term3) + - (transform->y->coeff[0][0][0][1] * term4); - - return (out); -} - -// This function prototype has been modified since the SDRS. -psSphereTransform* psSphereTransformAlloc(double NPlat, - double Xo, - double xo) -{ - psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); - - tmp->sinPhi = sin(NPlat); - tmp->cosPhi = cos(NPlat); - tmp->Xo = Xo; - tmp->xo = xo; - - return (tmp); -} - -/****************************************************************************** -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. - *****************************************************************************/ -psSphere* psSphereTransformApply(psSphere* out, - const psSphereTransform* transform, - const psSphere* coord) -{ - 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->xo; - sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi; - cosY = sqrt(1.0 - sinY * sinY); - sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y); - cosX = cos(y) * cos(dx) / cos(y); - - out->r = atan2(sinX, cosX) + transform->Xo; - out->d = atan2(sinY, cosY); - - return (out); -} - -psSphereTransform* psSphereTransformICRStoEcliptic(psTime time) -{ - 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 Xo = 0.0; - double xo = 0.0; - - return (psSphereTransformAlloc(phi, Xo, xo)); -} - -psSphereTransform* psSphereTransformEcliptictoICRS(psTime time) -{ - 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 Xo = 0.0; - double xo = 0.0; - - return (psSphereTransformAlloc(phi, Xo, xo)); -} - -psSphereTransform* psSphereTransformICRStoGalatic(void) -{ - return (psSphereTransformAlloc(62.6, 282.25, 33.0)); -} - -psSphereTransform* psSphereTransformGalatictoICRS(void) -{ - return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); -} - -float cot(float x) +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +// None + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +static float p_psCot(float x); +static float p_psArg(float x, float y); + +/****************************************************************************** +XXX: Do this with a macro. + *****************************************************************************/ +float p_psCot(float x) { return (1.0 / atan(x)); @@ -170,6 +70,6 @@ XXX: Verify this arc tan function. *****************************************************************************/ -float arg(float x, - float y) +float p_psArg(float x, + float y) { if (x > 0) { @@ -185,6 +85,145 @@ } - psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y); + psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y); return (0.0); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ +psPlane* psPlaneTransformApply(psPlane* out, + const psPlaneTransform* transform, + const psPlane* coords) +{ + if (out == NULL) { + out = (psPlane* ) psAlloc(sizeof(psPlane)); + } + out->x = transform->x->coeff[0][0] + + (transform->x->coeff[1][0] * coords->x) + + (transform->x->coeff[0][1] * coords->y); + + out->y = transform->y->coeff[0][0] + + (transform->y->coeff[1][0] * coords->x) + + (transform->y->coeff[0][1] * coords->y); + + return (out); +} + +/****************************************************************************** +This transformation takes into account parameters beyond an objects spatial +coordinates: term3 and term4 (magnitude and color). + *****************************************************************************/ +psPlane* psPlaneDistortApply(psPlane* out, + const psPlaneDistort* transform, + const psPlane* coords, + float color, + float magnitude) +{ + if (out == NULL) { + out = (psPlane* ) psAlloc(sizeof(psPlane)); + } + + out->x = transform->x->coeff[0][0][0][0] + + (transform->x->coeff[1][0][0][0] * coords->x) + + (transform->x->coeff[0][1][0][0] * coords->y) + + (transform->x->coeff[0][0][1][0] * color) + + (transform->x->coeff[0][0][0][1] * magnitude); + + out->y = transform->y->coeff[0][0][0][0] + + (transform->y->coeff[1][0][0][0] * coords->x) + + (transform->y->coeff[0][1][0][0] * coords->y) + + (transform->y->coeff[0][0][1][0] * color) + + (transform->y->coeff[0][0][0][1] * magnitude); + + return (out); +} + +/****************************************************************************** +This function prototype has been modified since the SDRS. + *****************************************************************************/ +psSphereTransform* psSphereTransformAlloc(double NPlat, + double Xo, + double xo) +{ + psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); + + tmp->sinPhi = sin(NPlat); + tmp->cosPhi = cos(NPlat); + tmp->Xo = Xo; + tmp->xo = xo; + + return (tmp); +} + +/****************************************************************************** +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. + *****************************************************************************/ +psSphere* psSphereTransformApply(psSphere* out, + const psSphereTransform* transform, + const psSphere* coord) +{ + 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->xo; + sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi; + cosY = sqrt(1.0 - sinY * sinY); + sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y); + cosX = cos(y) * cos(dx) / cos(y); + + out->r = atan2(sinX, cosX) + transform->Xo; + out->d = atan2(sinY, cosY); + + return (out); +} + +psSphereTransform* psSphereTransformICRStoEcliptic(psTime time) +{ + 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 Xo = 0.0; + double xo = 0.0; + + return (psSphereTransformAlloc(phi, Xo, xo)); +} + +psSphereTransform* psSphereTransformEcliptictoICRS(psTime time) +{ + 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 Xo = 0.0; + double xo = 0.0; + + return (psSphereTransformAlloc(phi, Xo, xo)); +} + +psSphereTransform* psSphereTransformICRStoGalatic(void) +{ + return (psSphereTransformAlloc(62.6, 282.25, 33.0)); +} + +psSphereTransform* psSphereTransformGalatictoICRS(void) +{ + return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); } @@ -201,5 +240,5 @@ if (projection->type == PS_PROJ_TAN) { - R = cot(coord->r) * (180.0 / M_PI); + R = p_psCot(coord->r) * (180.0 / M_PI); tmp->x = R * sin(coord->d); tmp->y = R * cos(coord->d); @@ -249,10 +288,10 @@ if (projection->type == PS_PROJ_TAN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = arg(-coord->y, coord->x); + tmp->d = p_psArg(-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->d = p_psArg(-coord->y, coord->x); tmp->r = acos((R * M_PI) / 180.0); @@ -271,5 +310,5 @@ 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->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); tmp->r = asin((coord->y * chu * M_PI) / 180.0); @@ -314,4 +353,5 @@ lin = psProject(position2, &proj); tmp = psDeproject(lin, &proj); + psFree(lin); // XXX: Do we need to convert units in tmp? @@ -349,8 +389,8 @@ /****************************************************************************** -XXX: Do I need to check for unacceptable transformation parameters? Maybe, +XXX: Do we need to check for unacceptable transformation parameters? Maybe, if the points are on the North/South Pole, etc? -XXX: Do I need to somehow scale this projection? +XXX: Do we need to somehow scale this projection? XXX: I copied the algorithm from the ADD exactly. Index: /trunk/psLib/src/astronomy/psAstrometry.c =================================================================== --- /trunk/psLib/src/astronomy/psAstrometry.c (revision 1531) +++ /trunk/psLib/src/astronomy/psAstrometry.c (revision 1532) @@ -8,12 +8,13 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-08-13 22:41:24 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-08-13 23:43:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ - +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ #include - #include "psFunctions.h" #include "psAstrometry.h" @@ -21,74 +22,31 @@ #include "psAbort.h" #include "slalib.h" - -psExposure* psExposureAlloc(double ra, double dec, double hourAngle, - double zenith, double azimuth, double localTime, float date, - float rotAngle, float temperature, float pressure, float humidity, - float exposureTime) -{ - psExposure* exp = psAlloc(sizeof(psExposure)); - - *(double *)&exp->ra = ra; - *(double *)&exp->dec = dec; - *(double *)&exp->hourAngle = hourAngle; - *(double *)&exp->zenith = zenith; - *(double *)&exp->azimuth = azimuth; - *(double *)&exp->localTime = localTime; - *(float *)&exp->date = date; - *(float *)&exp->rotAngle = rotAngle; - *(float *)&exp->temperature = temperature; - *(float *)&exp->pressure = pressure; - *(float *)&exp->humidity = humidity; - *(float *)&exp->exposureTime = exposureTime; - - return exp; -} -#define TBD 0.0 -/***************************************************************************** -XXX: Several of the input params to sla_aoppa() are currently undefined. -We are awaiting futher direction from IfA on this. - *****************************************************************************/ -psGrommit* psGrommitAlloc(const psExposure* exp) -{ - if (exp == NULL) { - psAbort(__func__, "the 'exp' parameter is NULL\n"); - } - - double date = TBD; // "mjd" in psExposure will become a psTime - // from which it will be possible to get UTC. - double dut = 0.0; - double elongm = TBD; - double phim = TBD; - double hm = TBD; - double xp = 0.0; - double yp = 0.0; - double tdk = exp->temperature; - double pmb = exp->pressure; - double rh = exp->humidity; - double wl = exp->wavelength; - double tlr = TBD; - double *AOPRMS = NULL; - - slaAoppa(date, dut, elongm, phim, hm, xp, yp, - tdk, pmb, rh, wl, tlr, AOPRMS); - - psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit)); - *(double *)&grommit->latitude = AOPRMS[0]; - *(double *)&grommit->sinLat = AOPRMS[1]; - *(double *)&grommit->cosLat = AOPRMS[2]; - *(double *)&grommit->height = AOPRMS[3]; - *(double *)&grommit->abberationMag = AOPRMS[4]; - *(double *)&grommit->temperature = AOPRMS[5]; - *(double *)&grommit->pressure = AOPRMS[6]; - *(double *)&grommit->humidity = AOPRMS[7]; - *(double *)&grommit->wavelength = AOPRMS[8]; - *(double *)&grommit->lapseRate = AOPRMS[9]; - *(double *)&grommit->refractA = AOPRMS[10]; - *(double *)&grommit->refractB = AOPRMS[11]; - *(double *)&grommit->longitudeOffset = AOPRMS[12]; - *(double *)&grommit->siderealTime = AOPRMS[13]; - - return (grommit); -} +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +// None + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ void p_psGrommitFree(psGrommit* grommit) @@ -97,33 +55,9 @@ } -psCell* psCellInFPA(const psPlane* fpaCoord, - const psFPA* FPA) -{ - psChip* tmpChip = NULL; - psPlane* chipCoord = NULL; - psCell* outCell = NULL; - - if (fpaCoord == NULL) { - psAbort(__func__, "input parameter fpaCoord is NULL."); - } - if (FPA == NULL) { - psAbort(__func__, "input parameter FPA is NULL."); - } - - // Determine which chip contains the fpaCoords. - tmpChip = psChipInFPA(fpaCoord, FPA); - - // Convert to those chip coordinates. - chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, tmpChip); - - // Determine which cell contains those chip coordinates. - outCell = psCellInChip(chipCoord, tmpChip); - - psFree(tmpChip); - psFree(chipCoord); - - return (outCell); -} - +/***************************************************************************** +p_psCheckValidImageCoords(): this is a private function which simply +determines if the supplied x,y coordinates are in the range for the supplied +psImage. + *****************************************************************************/ int p_psCheckValidImageCoords(double x, double y, @@ -142,382 +76,4 @@ return (1); -} - -psChip* psChipInFPA(const psPlane* fpaCoord, - const psFPA* FPA) -{ - psArray* chips = FPA->chips; - int nChips = chips->n; - psPlane* chipCoord = NULL; - psCell *tmpCell = NULL; - - if (fpaCoord == NULL) { - psAbort(__func__, "input parameter fpaCoord is NULL."); - } - if (FPA == NULL) { - psAbort(__func__, "input parameter FPA is NULL."); - } - - // Loop through every chip in this FPA. Convert the original - // FPA coordinates to chip coordinates for that chip. Then, - // determine if any cells in that chip contain those chip - // coordinates. - for (int i = 0; i < nChips; i++) { - psChip* tmpChip = chips->data[i]; - chipCoord = psPlaneTransformApply(chipCoord, tmpChip->fromFPA, - fpaCoord); - - tmpCell = psCellInChip(chipCoord, tmpChip); - if (tmpCell != NULL) { - psFree(chipCoord); - return(tmpChip); - } - psFree(chipCoord); - } - - return (NULL); -} - -psCell* psCellInChip(const psPlane* chipCoord, - const psChip* chip) -{ - psPlane* cellCoord = NULL; - psArray* cells; - - // We return NULL if either of the input parameters is NULL. - if (chipCoord == NULL) { - psAbort(__func__, "the 'chipCoord' parameter is NULL\n"); - } - if (chip == NULL) { - psAbort(__func__, "the 'chip' parameter is NULL\n"); - } - - cells = chip->cells; - if (cells == NULL) { - return NULL; - } - - // We loop over each cell in the chip. We transform the chipCoord into - // a cellCoord for that cell and determine if that cellCoord is valid. - // If so, then we return that cell. - for (int i = 0; i < cells->n; i++) { - psCell* tmpCell = (psCell* ) cells->data[i]; - psArray* readouts = tmpCell->readouts; - - if (readouts != NULL) { - for (int j = 0; j < readouts->n; j++) { - psReadout* tmpReadout = readouts->data[j]; - - cellCoord = psPlaneTransformApply(cellCoord, - tmpCell->fromChip, - chipCoord); - - if (p_psCheckValidImageCoords(cellCoord->x, - cellCoord->y, - tmpReadout->image)) { - psFree(cellCoord); - return (tmpCell); - } - } - } - } - - psFree(cellCoord); - return (NULL); -} - -psPlane* psCoordCellToChip(psPlane* outCoord, - const psPlane* inCoord, - const psCell* cell) -{ - if (inCoord == NULL) { - psAbort(__func__, "input parameter inCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - return (psPlaneTransformApply(outCoord, cell->toChip, inCoord)); -} - -psPlane* psCoordChipToFPA(psPlane* outCoord, - const psPlane* inCoord, - const psChip* chip) -{ - if (inCoord == NULL) { - psAbort(__func__, "input parameter inCoord is NULL."); - } - if (chip == NULL) { - psAbort(__func__, "input parameter chip is NULL."); - } - - return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord)); -} - -psPlane* psCoordFPAToTP(psPlane* outCoord, - const psPlane* inCoord, - double color, - double magnitude, - const psFPA* fpa) -{ - if (inCoord == NULL) { - psAbort(__func__, "input parameter inCoord is NULL."); - } - if (fpa == NULL) { - psAbort(__func__, "input parameter fpa is NULL."); - } - - return(psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, - color, magnitude)); -} - -/***************************************************************************** -XXX: What about units for the (x,y) coords? - *****************************************************************************/ -psSphere* psCoordTPToSky(psSphere* outSphere, - const psPlane* tpCoord, - const psGrommit* grommit) -{ - double AOB = 0.0; - double ZOB = 0.0; - double HOB = 0.0; - - if (tpCoord == NULL) { - psAbort(__func__, "input parameter tpCoord is NULL."); - } - if (grommit == NULL) { - psAbort(__func__, "input parameter grommit is NULL."); - } - if (outSphere == NULL) { - outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); - } - - slaAopqk(tpCoord->x, tpCoord->y, (double *) grommit, - &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d); - - return (outSphere); -} - -psPlane* psCoordCellToFPA(psPlane* fpaCoord, - const psPlane* cellCoord, - const psCell* cell) -{ - if (cellCoord == NULL) { - psAbort(__func__, "input parameter cellCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord)); -} - -psSphere* psCoordCellToSky(psSphere* skyCoord, - const psPlane* cellCoord, - double color, - double magnitude, - const psCell* cell) -{ - if (cellCoord == NULL) { - psAbort(__func__, "input parameter cellCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - psPlane* fpaCoord = NULL; - psPlane* tpCoord = NULL; - psFPA* parFPA = (cell->parent)->parent; - psGrommit* tmpGrommit = NULL; - - // Convert the input cell coordinates to FPA coordinates. - fpaCoord = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord); - - // Convert the FPA coordinates to tangent plane Coordinates. - tpCoord = psPlaneDistortApply(tpCoord, parFPA->toTangentPlane, - fpaCoord, color, magnitude); - - // Generate a grommit for this FPA. - tmpGrommit = psGrommitAlloc(parFPA->exposure); - - // Convert the tangent plane Coordinates to sky coordinates. - skyCoord = psCoordTPToSky(skyCoord, tpCoord, tmpGrommit); - - psFree(fpaCoord); - psFree(tpCoord); - psFree(tmpGrommit); - - return(skyCoord); -} - -psSphere* psCoordCellToSkyQuick(psSphere* outSphere, - const psPlane* cellCoord, - const psCell* cell) -{ - if (cellCoord == NULL) { - psAbort(__func__, "input parameter cellCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - psPlane *tpCoord = NULL; - psChip *chip = cell->parent; - psFPA *FPA = chip->parent; - psProjectionType oldProjectionType; - - if (outSphere == NULL) { - outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); - } - - // Determine the tangent plane coordinates. - tpCoord = psPlaneTransformApply(tpCoord, cell->toTP, cellCoord); - - // Save the old projection type and set the new projection type to TAN. - oldProjectionType = FPA->projection->type; - FPA->projection->type = PS_PROJ_TAN; - - // Deproject the tangent plane coordinates a sphere. - outSphere = psDeproject(tpCoord, FPA->projection); - - // Restore old projection type. Free memory. - FPA->projection->type = oldProjectionType; - psFree(tpCoord); - - return (outSphere); -} - -/***************************************************************************** -XXX: What about units for the (x,y) coords? - *****************************************************************************/ -psPlane* psCoordSkyToTP(psPlane* tpCoord, - const psSphere* in, - const psGrommit* grommit) -{ - if (in == NULL) { - psAbort(__func__, "input parameter in is NULL."); - } - if (grommit == NULL) { - psAbort(__func__, "input parameter grommit is NULL."); - } - if (tpCoord == NULL) { - tpCoord = (psPlane* ) psAlloc(sizeof(psPlane)); - } - - slaOapqk("RA", in->r, in->d, (double *) grommit, &tpCoord->x, &tpCoord->y); - - return(tpCoord); -} - - -psPlane* psCoordTPToFPA(psPlane* fpaCoord, - const psPlane* tpCoord, - double color, - double magnitude, - const psFPA* fpa) -{ - if (tpCoord == NULL) { - psAbort(__func__, "input parameter tpCoord is NULL."); - } - if (fpa == NULL) { - psAbort(__func__, "input parameter fpa is NULL."); - } - - return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, - tpCoord, color, magnitude)); -} - -psPlane* psCoordFPAToChip(psPlane* chipCoord, - const psPlane* fpaCoord, - const psChip* chip) -{ - if (fpaCoord == NULL) { - psAbort(__func__, "input parameter fpaCoord is NULL."); - } - if (chip == NULL) { - psAbort(__func__, "input parameter chip is NULL."); - } - - psFPA *FPA = chip->parent; - if (FPA == NULL) { - psAbort(__func__, "chip->parent is NULL"); - } - - // Determine which chip contains these FPA coordinates. - psChip *newChip = psChipInFPA(fpaCoord, FPA); - if (newChip == NULL) { - return(NULL); - } - - chipCoord = psPlaneTransformApply(chipCoord, newChip->fromFPA, fpaCoord); - psFree(newChip); - return(chipCoord); -} - -psPlane* psCoordChipToCell(psPlane* cellCoord, - const psPlane* chipCoord, - const psCell* cell) -{ - if (chipCoord == NULL) { - psAbort(__func__, "input parameter chipCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - psChip *chip = cell->parent; - if (chip == NULL) { - psAbort(__func__, "cell->parent is NULL"); - } - - // Determine which cell contains these FPA coordinates. - psCell *newCell = psCellInChip(chipCoord, chip); - if (newCell == NULL) { - return(NULL); - } - - cellCoord = psPlaneTransformApply(cellCoord, newCell->fromChip, chipCoord); - psFree(newCell); - return(cellCoord); -} - -/***************************************************************************** -XXX: Should we do anything to determine which cell contains the coordinates? - *****************************************************************************/ -psPlane* psCoordSkyToCell(psPlane* cellCoord, - const psSphere* skyCoord, - double color, - double magnitude, - const psCell* cell) -{ - if (skyCoord == NULL) { - psAbort(__func__, "input parameter skyCoord is NULL."); - } - if (cell == NULL) { - psAbort(__func__, "input parameter cell is NULL."); - } - - psChip *parChip = cell->parent; - psFPA *parFPA = parChip->parent; - psGrommit* grommit = parFPA->grommit; - - // Convert the skyCoords to tangent plane coords. - psPlane *tpCoord = psCoordSkyToTP(tpCoord, skyCoord, grommit); - - // Convert the tangent plane coords to FPA coords. - psPlane *fpaCoord = psCoordTPToFPA(fpaCoord, tpCoord, color, - magnitude, parFPA); - - // Convert the FPA coords to chip coords. - psPlane *chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, parChip); - - // Convert the chip coords to cell coords. - cellCoord = psCoordChipToCell(cellCoord, chipCoord, cell); - - psFree(tpCoord); - psFree(fpaCoord); - psFree(chipCoord); - - return (cellCoord); } @@ -602,4 +158,484 @@ } +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ +psExposure* psExposureAlloc(double ra, double dec, double hourAngle, + double zenith, double azimuth, double localTime, float date, + float rotAngle, float temperature, float pressure, float humidity, + float exposureTime) +{ + psExposure* exp = psAlloc(sizeof(psExposure)); + + *(double *)&exp->ra = ra; + *(double *)&exp->dec = dec; + *(double *)&exp->hourAngle = hourAngle; + *(double *)&exp->zenith = zenith; + *(double *)&exp->azimuth = azimuth; + *(double *)&exp->localTime = localTime; + *(float *)&exp->date = date; + *(float *)&exp->rotAngle = rotAngle; + *(float *)&exp->temperature = temperature; + *(float *)&exp->pressure = pressure; + *(float *)&exp->humidity = humidity; + *(float *)&exp->exposureTime = exposureTime; + + return exp; +} +#define TBD 0.0 +/***************************************************************************** +XXX: Several of the input params to sla_aoppa() are currently undefined. +We are awaiting futher direction from IfA on this. + *****************************************************************************/ +psGrommit* psGrommitAlloc(const psExposure* exp) +{ + if (exp == NULL) { + psAbort(__func__, "the 'exp' parameter is NULL\n"); + } + + double date = TBD; // "mjd" in psExposure will become a psTime + // from which it will be possible to get UTC. + double dut = 0.0; + double elongm = TBD; + double phim = TBD; + double hm = TBD; + double xp = 0.0; + double yp = 0.0; + double tdk = exp->temperature; + double pmb = exp->pressure; + double rh = exp->humidity; + double wl = exp->wavelength; + double tlr = TBD; + double *AOPRMS = NULL; + + slaAoppa(date, dut, elongm, phim, hm, xp, yp, + tdk, pmb, rh, wl, tlr, AOPRMS); + + psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit)); + *(double *)&grommit->latitude = AOPRMS[0]; + *(double *)&grommit->sinLat = AOPRMS[1]; + *(double *)&grommit->cosLat = AOPRMS[2]; + *(double *)&grommit->height = AOPRMS[3]; + *(double *)&grommit->abberationMag = AOPRMS[4]; + *(double *)&grommit->temperature = AOPRMS[5]; + *(double *)&grommit->pressure = AOPRMS[6]; + *(double *)&grommit->humidity = AOPRMS[7]; + *(double *)&grommit->wavelength = AOPRMS[8]; + *(double *)&grommit->lapseRate = AOPRMS[9]; + *(double *)&grommit->refractA = AOPRMS[10]; + *(double *)&grommit->refractB = AOPRMS[11]; + *(double *)&grommit->longitudeOffset = AOPRMS[12]; + *(double *)&grommit->siderealTime = AOPRMS[13]; + + return (grommit); +} + +psCell* psCellInFPA(const psPlane* fpaCoord, + const psFPA* FPA) +{ + psChip* tmpChip = NULL; + psPlane* chipCoord = NULL; + psCell* outCell = NULL; + + if (fpaCoord == NULL) { + psAbort(__func__, "input parameter fpaCoord is NULL."); + } + if (FPA == NULL) { + psAbort(__func__, "input parameter FPA is NULL."); + } + + // Determine which chip contains the fpaCoords. + tmpChip = psChipInFPA(fpaCoord, FPA); + + // Convert to those chip coordinates. + chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, tmpChip); + + // Determine which cell contains those chip coordinates. + outCell = psCellInChip(chipCoord, tmpChip); + + psFree(tmpChip); + psFree(chipCoord); + + return (outCell); +} + +psChip* psChipInFPA(const psPlane* fpaCoord, + const psFPA* FPA) +{ + psArray* chips = FPA->chips; + int nChips = chips->n; + psPlane* chipCoord = NULL; + psCell *tmpCell = NULL; + + if (fpaCoord == NULL) { + psAbort(__func__, "input parameter fpaCoord is NULL."); + } + if (FPA == NULL) { + psAbort(__func__, "input parameter FPA is NULL."); + } + + // Loop through every chip in this FPA. Convert the original + // FPA coordinates to chip coordinates for that chip. Then, + // determine if any cells in that chip contain those chip + // coordinates. + for (int i = 0; i < nChips; i++) { + psChip* tmpChip = chips->data[i]; + chipCoord = psPlaneTransformApply(chipCoord, tmpChip->fromFPA, + fpaCoord); + + tmpCell = psCellInChip(chipCoord, tmpChip); + if (tmpCell != NULL) { + psFree(chipCoord); + return(tmpChip); + } + psFree(chipCoord); + } + + return (NULL); +} + +psCell* psCellInChip(const psPlane* chipCoord, + const psChip* chip) +{ + psPlane* cellCoord = NULL; + psArray* cells; + + // We return NULL if either of the input parameters is NULL. + if (chipCoord == NULL) { + psAbort(__func__, "the 'chipCoord' parameter is NULL\n"); + } + if (chip == NULL) { + psAbort(__func__, "the 'chip' parameter is NULL\n"); + } + + cells = chip->cells; + if (cells == NULL) { + return NULL; + } + + // We loop over each cell in the chip. We transform the chipCoord into + // a cellCoord for that cell and determine if that cellCoord is valid. + // If so, then we return that cell. + for (int i = 0; i < cells->n; i++) { + psCell* tmpCell = (psCell* ) cells->data[i]; + psArray* readouts = tmpCell->readouts; + + if (readouts != NULL) { + for (int j = 0; j < readouts->n; j++) { + psReadout* tmpReadout = readouts->data[j]; + + cellCoord = psPlaneTransformApply(cellCoord, + tmpCell->fromChip, + chipCoord); + + if (p_psCheckValidImageCoords(cellCoord->x, + cellCoord->y, + tmpReadout->image)) { + psFree(cellCoord); + return (tmpCell); + } + } + } + } + + psFree(cellCoord); + return (NULL); +} + +psPlane* psCoordCellToChip(psPlane* outCoord, + const psPlane* inCoord, + const psCell* cell) +{ + if (inCoord == NULL) { + psAbort(__func__, "input parameter inCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + return (psPlaneTransformApply(outCoord, cell->toChip, inCoord)); +} + +psPlane* psCoordChipToFPA(psPlane* outCoord, + const psPlane* inCoord, + const psChip* chip) +{ + if (inCoord == NULL) { + psAbort(__func__, "input parameter inCoord is NULL."); + } + if (chip == NULL) { + psAbort(__func__, "input parameter chip is NULL."); + } + + return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord)); +} + +psPlane* psCoordFPAToTP(psPlane* outCoord, + const psPlane* inCoord, + double color, + double magnitude, + const psFPA* fpa) +{ + if (inCoord == NULL) { + psAbort(__func__, "input parameter inCoord is NULL."); + } + if (fpa == NULL) { + psAbort(__func__, "input parameter fpa is NULL."); + } + + return(psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, + color, magnitude)); +} + +/***************************************************************************** +XXX: What about units for the (x,y) coords? + *****************************************************************************/ +psSphere* psCoordTPToSky(psSphere* outSphere, + const psPlane* tpCoord, + const psGrommit* grommit) +{ + double AOB = 0.0; + double ZOB = 0.0; + double HOB = 0.0; + + if (tpCoord == NULL) { + psAbort(__func__, "input parameter tpCoord is NULL."); + } + if (grommit == NULL) { + psAbort(__func__, "input parameter grommit is NULL."); + } + if (outSphere == NULL) { + outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); + } + + slaAopqk(tpCoord->x, tpCoord->y, (double *) grommit, + &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d); + + return (outSphere); +} + +psPlane* psCoordCellToFPA(psPlane* fpaCoord, + const psPlane* cellCoord, + const psCell* cell) +{ + if (cellCoord == NULL) { + psAbort(__func__, "input parameter cellCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord)); +} + +psSphere* psCoordCellToSky(psSphere* skyCoord, + const psPlane* cellCoord, + double color, + double magnitude, + const psCell* cell) +{ + if (cellCoord == NULL) { + psAbort(__func__, "input parameter cellCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + psPlane* fpaCoord = NULL; + psPlane* tpCoord = NULL; + psFPA* parFPA = (cell->parent)->parent; + psGrommit* tmpGrommit = NULL; + + // Convert the input cell coordinates to FPA coordinates. + fpaCoord = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord); + + // Convert the FPA coordinates to tangent plane Coordinates. + tpCoord = psPlaneDistortApply(tpCoord, parFPA->toTangentPlane, + fpaCoord, color, magnitude); + + // Generate a grommit for this FPA. + tmpGrommit = psGrommitAlloc(parFPA->exposure); + + // Convert the tangent plane Coordinates to sky coordinates. + skyCoord = psCoordTPToSky(skyCoord, tpCoord, tmpGrommit); + + psFree(fpaCoord); + psFree(tpCoord); + psFree(tmpGrommit); + + return(skyCoord); +} + +psSphere* psCoordCellToSkyQuick(psSphere* outSphere, + const psPlane* cellCoord, + const psCell* cell) +{ + if (cellCoord == NULL) { + psAbort(__func__, "input parameter cellCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + psPlane *tpCoord = NULL; + psChip *chip = cell->parent; + psFPA *FPA = chip->parent; + psProjectionType oldProjectionType; + + if (outSphere == NULL) { + outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); + } + + // Determine the tangent plane coordinates. + tpCoord = psPlaneTransformApply(tpCoord, cell->toTP, cellCoord); + + // Save the old projection type and set the new projection type to TAN. + oldProjectionType = FPA->projection->type; + FPA->projection->type = PS_PROJ_TAN; + + // Deproject the tangent plane coordinates a sphere. + outSphere = psDeproject(tpCoord, FPA->projection); + + // Restore old projection type. Free memory. + FPA->projection->type = oldProjectionType; + psFree(tpCoord); + + return (outSphere); +} + +/***************************************************************************** +XXX: What about units for the (x,y) coords? + *****************************************************************************/ +psPlane* psCoordSkyToTP(psPlane* tpCoord, + const psSphere* in, + const psGrommit* grommit) +{ + if (in == NULL) { + psAbort(__func__, "input parameter in is NULL."); + } + if (grommit == NULL) { + psAbort(__func__, "input parameter grommit is NULL."); + } + if (tpCoord == NULL) { + tpCoord = (psPlane* ) psAlloc(sizeof(psPlane)); + } + + slaOapqk("RA", in->r, in->d, (double *) grommit, &tpCoord->x, &tpCoord->y); + + return(tpCoord); +} + + +psPlane* psCoordTPToFPA(psPlane* fpaCoord, + const psPlane* tpCoord, + double color, + double magnitude, + const psFPA* fpa) +{ + if (tpCoord == NULL) { + psAbort(__func__, "input parameter tpCoord is NULL."); + } + if (fpa == NULL) { + psAbort(__func__, "input parameter fpa is NULL."); + } + + return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, + tpCoord, color, magnitude)); +} + +psPlane* psCoordFPAToChip(psPlane* chipCoord, + const psPlane* fpaCoord, + const psChip* chip) +{ + if (fpaCoord == NULL) { + psAbort(__func__, "input parameter fpaCoord is NULL."); + } + if (chip == NULL) { + psAbort(__func__, "input parameter chip is NULL."); + } + + psFPA *FPA = chip->parent; + if (FPA == NULL) { + psAbort(__func__, "chip->parent is NULL"); + } + + // Determine which chip contains these FPA coordinates. + psChip *newChip = psChipInFPA(fpaCoord, FPA); + if (newChip == NULL) { + return(NULL); + } + + chipCoord = psPlaneTransformApply(chipCoord, newChip->fromFPA, fpaCoord); + psFree(newChip); + return(chipCoord); +} + +psPlane* psCoordChipToCell(psPlane* cellCoord, + const psPlane* chipCoord, + const psCell* cell) +{ + if (chipCoord == NULL) { + psAbort(__func__, "input parameter chipCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + psChip *chip = cell->parent; + if (chip == NULL) { + psAbort(__func__, "cell->parent is NULL"); + } + + // Determine which cell contains these FPA coordinates. + psCell *newCell = psCellInChip(chipCoord, chip); + if (newCell == NULL) { + return(NULL); + } + + cellCoord = psPlaneTransformApply(cellCoord, newCell->fromChip, chipCoord); + psFree(newCell); + return(cellCoord); +} + +/***************************************************************************** +XXX: Should we do anything to determine which cell contains the coordinates? + *****************************************************************************/ +psPlane* psCoordSkyToCell(psPlane* cellCoord, + const psSphere* skyCoord, + double color, + double magnitude, + const psCell* cell) +{ + if (skyCoord == NULL) { + psAbort(__func__, "input parameter skyCoord is NULL."); + } + if (cell == NULL) { + psAbort(__func__, "input parameter cell is NULL."); + } + + psChip *parChip = cell->parent; + psFPA *parFPA = parChip->parent; + psGrommit* grommit = parFPA->grommit; + + // Convert the skyCoords to tangent plane coords. + psPlane *tpCoord = psCoordSkyToTP(tpCoord, skyCoord, grommit); + + // Convert the tangent plane coords to FPA coords. + psPlane *fpaCoord = psCoordTPToFPA(fpaCoord, tpCoord, color, + magnitude, parFPA); + + // Convert the FPA coords to chip coords. + psPlane *chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, parChip); + + // Convert the chip coords to cell coords. + cellCoord = psCoordChipToCell(cellCoord, chipCoord, cell); + + psFree(tpCoord); + psFree(fpaCoord); + psFree(chipCoord); + + return (cellCoord); +} + /***************************************************************************** XXX: Should we do anything to determine which cell contains the coordinates? Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 1531) +++ /trunk/psLib/src/astronomy/psCoord.c (revision 1532) @@ -10,10 +10,12 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-08-13 23:33:13 $ +* @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-08-13 23:43:29 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ - +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ #include "psType.h" #include "psImage.h" @@ -26,141 +28,39 @@ #include #include - -static float cot(float x); -static float arg(float x, float y); - -psPlane* psPlaneTransformApply(psPlane* out, - const psPlaneTransform* transform, - const psPlane* coords) -{ - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - out->x = transform->x->coeff[0][0] + - (transform->x->coeff[1][0] * coords->x) + - (transform->x->coeff[0][1] * coords->y); - - out->y = transform->y->coeff[0][0] + - (transform->y->coeff[1][0] * coords->x) + - (transform->y->coeff[0][1] * coords->y); - - return (out); -} - -// This transformation takes into account parameters beyond an objects -// spatial coordinates: term3 and term4 (magnitude and color). -psPlane* psPlaneDistortApply(psPlane* out, - const psPlaneDistort* transform, - const psPlane* coords, - float term3, - float term4) -{ - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - - out->x = transform->x->coeff[0][0][0][0] + - (transform->x->coeff[1][0][0][0] * coords->x) + - (transform->x->coeff[0][1][0][0] * coords->y) + - (transform->x->coeff[0][0][1][0] * term3) + - (transform->x->coeff[0][0][0][1] * term4); - - out->y = transform->y->coeff[0][0][0][0] + - (transform->y->coeff[1][0][0][0] * coords->x) + - (transform->y->coeff[0][1][0][0] * coords->y) + - (transform->y->coeff[0][0][1][0] * term3) + - (transform->y->coeff[0][0][0][1] * term4); - - return (out); -} - -// This function prototype has been modified since the SDRS. -psSphereTransform* psSphereTransformAlloc(double NPlat, - double Xo, - double xo) -{ - psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); - - tmp->sinPhi = sin(NPlat); - tmp->cosPhi = cos(NPlat); - tmp->Xo = Xo; - tmp->xo = xo; - - return (tmp); -} - -/****************************************************************************** -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. - *****************************************************************************/ -psSphere* psSphereTransformApply(psSphere* out, - const psSphereTransform* transform, - const psSphere* coord) -{ - 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->xo; - sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi; - cosY = sqrt(1.0 - sinY * sinY); - sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y); - cosX = cos(y) * cos(dx) / cos(y); - - out->r = atan2(sinX, cosX) + transform->Xo; - out->d = atan2(sinY, cosY); - - return (out); -} - -psSphereTransform* psSphereTransformICRStoEcliptic(psTime time) -{ - 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 Xo = 0.0; - double xo = 0.0; - - return (psSphereTransformAlloc(phi, Xo, xo)); -} - -psSphereTransform* psSphereTransformEcliptictoICRS(psTime time) -{ - 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 Xo = 0.0; - double xo = 0.0; - - return (psSphereTransformAlloc(phi, Xo, xo)); -} - -psSphereTransform* psSphereTransformICRStoGalatic(void) -{ - return (psSphereTransformAlloc(62.6, 282.25, 33.0)); -} - -psSphereTransform* psSphereTransformGalatictoICRS(void) -{ - return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); -} - -float cot(float x) +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +// None + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +static float p_psCot(float x); +static float p_psArg(float x, float y); + +/****************************************************************************** +XXX: Do this with a macro. + *****************************************************************************/ +float p_psCot(float x) { return (1.0 / atan(x)); @@ -170,6 +70,6 @@ XXX: Verify this arc tan function. *****************************************************************************/ -float arg(float x, - float y) +float p_psArg(float x, + float y) { if (x > 0) { @@ -185,6 +85,145 @@ } - psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y); + psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y); return (0.0); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ +psPlane* psPlaneTransformApply(psPlane* out, + const psPlaneTransform* transform, + const psPlane* coords) +{ + if (out == NULL) { + out = (psPlane* ) psAlloc(sizeof(psPlane)); + } + out->x = transform->x->coeff[0][0] + + (transform->x->coeff[1][0] * coords->x) + + (transform->x->coeff[0][1] * coords->y); + + out->y = transform->y->coeff[0][0] + + (transform->y->coeff[1][0] * coords->x) + + (transform->y->coeff[0][1] * coords->y); + + return (out); +} + +/****************************************************************************** +This transformation takes into account parameters beyond an objects spatial +coordinates: term3 and term4 (magnitude and color). + *****************************************************************************/ +psPlane* psPlaneDistortApply(psPlane* out, + const psPlaneDistort* transform, + const psPlane* coords, + float color, + float magnitude) +{ + if (out == NULL) { + out = (psPlane* ) psAlloc(sizeof(psPlane)); + } + + out->x = transform->x->coeff[0][0][0][0] + + (transform->x->coeff[1][0][0][0] * coords->x) + + (transform->x->coeff[0][1][0][0] * coords->y) + + (transform->x->coeff[0][0][1][0] * color) + + (transform->x->coeff[0][0][0][1] * magnitude); + + out->y = transform->y->coeff[0][0][0][0] + + (transform->y->coeff[1][0][0][0] * coords->x) + + (transform->y->coeff[0][1][0][0] * coords->y) + + (transform->y->coeff[0][0][1][0] * color) + + (transform->y->coeff[0][0][0][1] * magnitude); + + return (out); +} + +/****************************************************************************** +This function prototype has been modified since the SDRS. + *****************************************************************************/ +psSphereTransform* psSphereTransformAlloc(double NPlat, + double Xo, + double xo) +{ + psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform)); + + tmp->sinPhi = sin(NPlat); + tmp->cosPhi = cos(NPlat); + tmp->Xo = Xo; + tmp->xo = xo; + + return (tmp); +} + +/****************************************************************************** +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. + *****************************************************************************/ +psSphere* psSphereTransformApply(psSphere* out, + const psSphereTransform* transform, + const psSphere* coord) +{ + 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->xo; + sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi; + cosY = sqrt(1.0 - sinY * sinY); + sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y); + cosX = cos(y) * cos(dx) / cos(y); + + out->r = atan2(sinX, cosX) + transform->Xo; + out->d = atan2(sinY, cosY); + + return (out); +} + +psSphereTransform* psSphereTransformICRStoEcliptic(psTime time) +{ + 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 Xo = 0.0; + double xo = 0.0; + + return (psSphereTransformAlloc(phi, Xo, xo)); +} + +psSphereTransform* psSphereTransformEcliptictoICRS(psTime time) +{ + 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 Xo = 0.0; + double xo = 0.0; + + return (psSphereTransformAlloc(phi, Xo, xo)); +} + +psSphereTransform* psSphereTransformICRStoGalatic(void) +{ + return (psSphereTransformAlloc(62.6, 282.25, 33.0)); +} + +psSphereTransform* psSphereTransformGalatictoICRS(void) +{ + return (psSphereTransformAlloc(-62.6, 33.0, 282.25)); } @@ -201,5 +240,5 @@ if (projection->type == PS_PROJ_TAN) { - R = cot(coord->r) * (180.0 / M_PI); + R = p_psCot(coord->r) * (180.0 / M_PI); tmp->x = R * sin(coord->d); tmp->y = R * cos(coord->d); @@ -249,10 +288,10 @@ if (projection->type == PS_PROJ_TAN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = arg(-coord->y, coord->x); + tmp->d = p_psArg(-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->d = p_psArg(-coord->y, coord->x); tmp->r = acos((R * M_PI) / 180.0); @@ -271,5 +310,5 @@ 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->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); tmp->r = asin((coord->y * chu * M_PI) / 180.0); @@ -314,4 +353,5 @@ lin = psProject(position2, &proj); tmp = psDeproject(lin, &proj); + psFree(lin); // XXX: Do we need to convert units in tmp? @@ -349,8 +389,8 @@ /****************************************************************************** -XXX: Do I need to check for unacceptable transformation parameters? Maybe, +XXX: Do we need to check for unacceptable transformation parameters? Maybe, if the points are on the North/South Pole, etc? -XXX: Do I need to somehow scale this projection? +XXX: Do we need to somehow scale this projection? XXX: I copied the algorithm from the ADD exactly.