Index: /trunk/psModules/src/psAstrometry.c =================================================================== --- /trunk/psModules/src/psAstrometry.c (revision 4577) +++ /trunk/psModules/src/psAstrometry.c (revision 4577) @@ -0,0 +1,826 @@ +/** @file psAstrometry.c + * + * @brief This file defines the basic types for astronomical coordinate + * transformation + * + * @ingroup AstroImage + * + * @author GLG, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-18 18:48:29 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include + +#include "psFunctions.h" +#include "psAstrometry.h" +#include "psMemory.h" +#include "psError.h" +#include "psConstants.h" +#include "psAstronomyErrors.h" +#include "psMatrix.h" +#include "psTrace.h" +#include "psLogMsg.h" + +/***************************************************************************** +checkValidImageCoords(): this is a private function which simply +determines if the supplied x,y coordinates are in the range for the supplied +psImage. + *****************************************************************************/ +static psS32 checkValidImageCoords(double x, + double y, + psImage* tmpImage) +{ + PS_ASSERT_IMAGE_NON_NULL(tmpImage, 0); + + if ((x < 0.0) || (x > (double)tmpImage->numCols) || + (y < 0.0) || (y > (double)tmpImage->numRows)) { + return (0); + } + + return (1); +} + + +static void FPAFree(psFPA* fpa) +{ + if (fpa != NULL) { + psFree(fpa->chips); + psFree(fpa->grommit); + psFree(fpa->exposure); + psFree(fpa->metadata); + psFree(fpa->fromTangentPlane); + psFree(fpa->toTangentPlane); + psFree(fpa->pattern); + psFree(fpa->colorPlus); + psFree(fpa->colorMinus); + psFree(fpa->projection); + } +} + +static void chipFree(psChip* chip) +{ + if (chip != NULL) { + psFree(chip->cells); + psFree(chip->metadata); + psFree(chip->toFPA); + psFree(chip->fromFPA); + } +} + +static void cellFree(psCell* cell) +{ + if (cell != NULL) { + psFree(cell->readouts); + psFree(cell->metadata); + psFree(cell->toChip); + psFree(cell->fromChip); + psFree(cell->toFPA); + psFree(cell->toTP); + psFree(cell->toSky); + } +} + +static void readoutFree(psReadout* readout) +{ + if (readout != NULL) { + psFree(readout->image); + psFree(readout->mask); + psFree(readout->objects); + psFree(readout->metadata); + } +} + +static void observatoryFree(psObservatory* obs) +{ + if (obs != NULL) { + psFree(obs->name); + } +} + +static void exposureFree(psExposure* exp) +{ + if (exp != NULL) { + psFree(exp->time); + psFree(exp->observatory); + psFree(exp->cameraName); + psFree(exp->telescopeName); + } +} + +static void fixedPatternFree(psFixedPattern* fp) +{ + if (fp != NULL) { + for (psS32 i = 0; i < fp->p_ps_xRows; i++) { + psFree(fp->x[i]); + } + + for (psS32 j = 0; j < fp->p_ps_yRows; j++) { + psFree(fp->y[j]); + } + + psFree(fp->x); + psFree(fp->y); + } +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +/* + * XXX: Verify that you interpreted the SDR correctly. + * + * XXX: This assumes that x,y must be of type F64 + */ +psFixedPattern* psFixedPatternAlloc(double x0, + double y0, + double xScale, + double yScale, + const psImage *x, + const psImage *y) +{ + psFixedPattern *tmp; + psS32 i; + psS32 j; + + PS_ASSERT_IMAGE_NON_NULL(x, NULL); + PS_ASSERT_IMAGE_NON_NULL(y, NULL); + PS_ASSERT_IMAGE_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_IMAGE_TYPE(y, PS_TYPE_F64, NULL); + + tmp = (psFixedPattern *) psAlloc(sizeof(psFixedPattern)); + // XXX: Is this correct? + tmp->nX = (x->numCols * x->numRows); + tmp->nY = (y->numCols * y->numRows); + tmp->x0 = x0; + tmp->y0 = y0; + tmp->xScale = xScale; + tmp->yScale = yScale; + tmp->p_ps_xRows = x->numRows; + tmp->p_ps_xCols = x->numCols; + tmp->p_ps_yRows = y->numRows; + tmp->p_ps_yCols = y->numCols; + tmp->x = (double **) psAlloc(x->numRows * sizeof(double *)); + for (i=0;inumRows;i++) { + (tmp->x)[i] = (double *) psAlloc(x->numCols * sizeof(double)); + } + for (i=0;inumRows;i++) { + for (j=0;jnumCols;j++) { + (tmp->x)[i][j] = x->data.F64[i][j]; + } + } + + tmp->y = (double **) psAlloc(y->numRows * sizeof(double *)); + for (i=0;inumRows;i++) { + (tmp->y)[i] = (double *) psAlloc(y->numCols * sizeof(double)); + } + for (i=0;inumRows;i++) { + for (j=0;jnumCols;j++) { + (tmp->y)[i][j] = y->data.F64[i][j]; + } + } + + psMemSetDeallocator(tmp,(psFreeFunc)fixedPatternFree); + + return(tmp); +} + + +psExposure* psExposureAlloc(double ra, + double dec, + double hourAngle, + double zenithDistance, + double azimuth, + const psTime* time, + float rotAngle, + float temperature, + float pressure, + float humidity, + float exposureTime, + float wavelength, + const psObservatory* observatory) +{ + PS_ASSERT_PTR_NON_NULL(observatory, NULL); + + psExposure* exp = psAlloc(sizeof(psExposure)); + *(double *)&exp->ra = ra; + *(double *)&exp->dec = dec; + *(double *)&exp->hourAngle = hourAngle; + *(double *)&exp->zenithDistance = zenithDistance; + *(double *)&exp->azimuth = azimuth; + *(float *)&exp->rotAngle = rotAngle; + *(float *)&exp->temperature = temperature; + *(float *)&exp->pressure = pressure; + *(float *)&exp->humidity = humidity; + *(float *)&exp->exposureTime = exposureTime; + *(float *)&exp->wavelength = wavelength; + + exp->time = psMemIncrRefCounter((psPtr)time); + exp->observatory = psMemIncrRefCounter((psPtr)observatory); + + // XXX: how is this value derived? + *(double *)&exp->lst = psTimeToLMST((psTime*)time,observatory->longitude); + *(float *)&exp->positionAngle = 0.0f; // XXX: need input, see Bug #207 + *(float *)&exp->parallacticAngle = 0.0f; // XXX: need input, see Bug #207 + *(float *)&exp->airmass = 0.0f; // XXX: needs calculation! = slaAirmas(zenithDistance); + *(float *)&exp->parallacticFactor = 0.0f; + exp->cameraName = NULL; + exp->telescopeName = NULL; + + psMemSetDeallocator(exp,(psFreeFunc)exposureFree); + + return exp; +} + +psObservatory* psObservatoryAlloc(const char* name, + double latitude, + double longitude, + double height, + double tlr) +{ + psObservatory* obs = psAlloc(sizeof(psObservatory)); + + if (name == NULL) { + obs->name = NULL; + } else { + obs->name = psAlloc(strlen(name)+1); + strcpy((char*)obs->name, name); + } + + *(double *)&obs->latitude = latitude; + *(double *)&obs->longitude = longitude; + *(double *)&obs->height = height; + *(double *)&obs->tlr = tlr; + + psMemSetDeallocator(obs,(psFreeFunc)observatoryFree); + + return obs; +} + +psFPA* psFPAAlloc(psS32 nChips, + const psExposure* exp) +{ + PS_ASSERT_INT_NONNEGATIVE(nChips, NULL); + + psFPA* newFPA = psAlloc(sizeof(psFPA)); + + // create array of NULL chips of the size nChips + newFPA->chips = psArrayAlloc(nChips); + psPtr* chips = newFPA->chips->data; + for (psS32 i=0;ichips->n = 0; // per requirement + + newFPA->metadata = NULL; + newFPA->fromTangentPlane = NULL; + newFPA->toTangentPlane = NULL; + newFPA->pattern = NULL; + + if (exp != NULL) { + newFPA->exposure = psMemIncrRefCounter((psExposure*)exp); + newFPA->grommit = psGrommitAlloc(exp); + } else { + newFPA->exposure = NULL; + newFPA->grommit = NULL; + } + + newFPA->colorPlus = NULL; + newFPA->colorMinus = NULL; + newFPA->projection = NULL; + + newFPA->rmsX = 0.0f; + newFPA->rmsY = 0.0f; + newFPA->chi2 = 0.0f; + + psMemSetDeallocator(newFPA,(psFreeFunc)FPAFree); + + return newFPA; +} + +/* + * psChip constructor + */ +psChip* psChipAlloc(psS32 nCells, + psFPA *parentFPA) +{ + PS_ASSERT_INT_NONNEGATIVE(nCells, NULL); + + psChip* chip = psAlloc(sizeof(psChip)); + + // create array of NULL psCells + int n = (nCells > 0) ? nCells : 1; + chip->cells = psArrayAlloc(n); + psPtr* cells = chip->cells->data; + for (psS32 i=0;icells->n = 0; // per requirement + + *(int*)&chip->row0 = 0; + *(int*)&chip->col0 = 0; + + chip->metadata = NULL; + + chip->toFPA = NULL; + chip->fromFPA = NULL; + + chip->parent = parentFPA; + + psMemSetDeallocator(chip,(psFreeFunc)chipFree); + + return chip; + +} + +/* + * psCell constructor + */ +psCell* psCellAlloc(psS32 nReadouts, + psChip* parentChip) +{ + PS_ASSERT_INT_NONNEGATIVE(nReadouts, NULL); + + psCell* cell = psAlloc(sizeof(psCell)); + + // create array of NULL psReadouts + int n = (nReadouts > 0) ? nReadouts : 1; + cell->readouts = psArrayAlloc(n); + psPtr* readouts = cell->readouts->data; + for (psS32 i=0;ireadouts->n = 0; // per requirement + + *(int*)&cell->row0 = 0; + *(int*)&cell->col0 = 0; + + cell->metadata = NULL; + + cell->toChip = NULL; + cell->fromChip = NULL; + cell->toFPA = NULL; + cell->toTP = NULL; + cell->toSky = NULL; + + cell->parent = parentChip; + + psMemSetDeallocator(cell,(psFreeFunc)cellFree); + + return cell; + + +} + +psReadout* psReadoutAlloc() +{ + psReadout* readout = psAlloc(sizeof(psReadout)); + + *(psU32*)&readout->colBins = 1; + *(psU32*)&readout->rowBins = 1; + *(psU32*)&readout->rowParity = 0; + *(psU32*)&readout->colParity = 0; + *(psS32*)&readout->col0 = 0; + *(psS32*)&readout->row0 = 0; + + readout->image = NULL; + readout->mask = NULL; + readout->objects = NULL; + readout->metadata = NULL; + + psMemSetDeallocator(readout,(psFreeFunc)readoutFree); + + return readout; +} + +psGrommit* psGrommitAlloc(const psExposure* exp) +{ + PS_ASSERT_PTR_NON_NULL(exp, NULL); + + psSphere* polarMotion = p_psTimeGetPoleCoords(exp->time); + + psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit)); + + *(double*)&grommit->latitude = exp->observatory->latitude; + *(double*)&grommit->longitude = exp->observatory->longitude; + *(double*)&grommit->height = exp->observatory->height; + *(double*)&grommit->abberationMag = 0.0; // XXX: need to figure out what to set here. + *(double*)&grommit->temperature = exp->temperature; + *(double*)&grommit->pressure = exp->pressure; + *(double*)&grommit->humidity = exp->humidity; + *(double*)&grommit->wavelength = exp->wavelength; + *(double*)&grommit->lapseRate = exp->observatory->tlr; + *(double*)&grommit->refractA = polarMotion->r; // XXX: need to figure out what to set here too. + *(double*)&grommit->refractB = polarMotion->d; // XXX: need to figure out what to set here too. + *(double*)&grommit->siderealTime = psTimeToMJD(exp->time); // XXX: this is probably not correct + + psFree(polarMotion); + + return (grommit); +} + +psCell* psCellInFPA(const psPlane* fpaCoord, + const psFPA* FPA) +{ + PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL); + PS_ASSERT_PTR_NON_NULL(FPA, NULL); + + psChip* tmpChip = NULL; + psPlane chipCoord; + psCell* outCell = NULL; + + // Determine which chip contains the fpaCoords. + tmpChip = psChipInFPA(fpaCoord, FPA); + if (tmpChip == NULL) { + return(NULL); + } + + // Convert to those chip coordinates. + psCoordFPAToChip(&chipCoord, fpaCoord, tmpChip); + + // Determine which cell contains those chip coordinates. + outCell = psCellInChip(&chipCoord, tmpChip); + + return (outCell); +} + +psChip* psChipInFPA(const psPlane* fpaCoord, + const psFPA* FPA) +{ + PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL); + PS_ASSERT_PTR_NON_NULL(FPA, NULL); + PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL); + + psArray* chips = FPA->chips; + psS32 nChips = chips->n; + psPlane chipCoord; + psCell *tmpCell = 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 (psS32 i = 0; i < nChips; i++) { + psChip* tmpChip = chips->data[i]; + PS_ASSERT_PTR_NON_NULL(tmpChip, NULL); + PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL); + + psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord); + + tmpCell = psCellInChip(&chipCoord, tmpChip); + if (tmpCell != NULL) { + return(tmpChip); + } + } + + // XXX: Print warning here? + return (NULL); +} + +psCell* psCellInChip(const psPlane* chipCoord, + const psChip* chip) +{ + PS_ASSERT_PTR_NON_NULL(chipCoord, NULL); + PS_ASSERT_PTR_NON_NULL(chip, NULL); + + psPlane cellCoord; + psArray* cells; + + 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 (psS32 i = 0; i < cells->n; i++) { + psCell* tmpCell = (psCell* ) cells->data[i]; + PS_ASSERT_PTR_NON_NULL(tmpCell, NULL); + PS_ASSERT_PTR_NON_NULL(tmpCell->fromChip, NULL); + psArray* readouts = tmpCell->readouts; + + if (readouts != NULL) { + for (psS32 j = 0; j < readouts->n; j++) { + psReadout* tmpReadout = readouts->data[j]; + PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL); + + psPlaneTransformApply(&cellCoord, + tmpCell->fromChip, + chipCoord); + + if (checkValidImageCoords(cellCoord.x, + cellCoord.y, + tmpReadout->image)) { + return (tmpCell); + } + } + } + } + + return (NULL); +} + +psPlane* psCoordCellToChip(psPlane* outCoord, + const psPlane* inCoord, + const psCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(inCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + + return (psPlaneTransformApply(outCoord, cell->toChip, inCoord)); +} + +psPlane* psCoordChipToFPA(psPlane* outCoord, + const psPlane* inCoord, + const psChip* chip) +{ + PS_ASSERT_PTR_NON_NULL(inCoord, NULL); + PS_ASSERT_PTR_NON_NULL(chip, NULL); + + return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord)); +} + +psPlane* psCoordFPAToTP(psPlane* outCoord, + const psPlane* inCoord, + double color, + double magnitude, + const psFPA* fpa) +{ + PS_ASSERT_PTR_NON_NULL(inCoord, NULL); + PS_ASSERT_PTR_NON_NULL(fpa, 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) +{ + PS_ASSERT_PTR_NON_NULL(tpCoord, NULL); + PS_ASSERT_PTR_NON_NULL(grommit, NULL); + + if (outSphere == NULL) { + outSphere = (psSphere* ) psAlloc(sizeof(psSphere)); + } + + // XXX: this was done by a SLALIB call -- needs to be reimplemented + psWarning("Warning! psCoordTPToSky functionality is no longer implemented"); + /* 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) +{ + PS_ASSERT_PTR_NON_NULL(cellCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + + return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord)); +} + +psSphere* psCoordCellToSky(psSphere* skyCoord, + const psPlane* cellCoord, + double color, + double magnitude, + const psCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(cellCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent->exposure, 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) +{ + PS_ASSERT_PTR_NON_NULL(cellCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL); + if (cell->toSky) { + // XXX: Should we use toTP or toSky? + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psCoordCellToSkyQuick(): The cell->toSky transform is ignored. The cell->toTP transform is being used."); + } + + 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(NULL, 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) +{ + PS_ASSERT_PTR_NON_NULL(in, NULL); + PS_ASSERT_PTR_NON_NULL(grommit, NULL); + + // char* type = "RA"; + + if (tpCoord == NULL) { + tpCoord = (psPlane* ) psAlloc(sizeof(psPlane)); + } + + // XXX: this was done by a SLALIB call -- needs to be reimplemented + psWarning("Warning! psCoordSkyToTP functionality is no longer implemented"); + /* slaOapqk(type, 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) +{ + PS_ASSERT_PTR_NON_NULL(tpCoord, NULL); + PS_ASSERT_PTR_NON_NULL(fpa, NULL); + PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL); + + return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, + tpCoord, color, magnitude)); +} + +psPlane* psCoordFPAToChip(psPlane* chipCoord, + const psPlane* fpaCoord, + const psChip* chip) +{ + PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL); + PS_ASSERT_PTR_NON_NULL(chip, NULL); + PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL); + + chipCoord = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord); + return(chipCoord); +} + +psPlane* psCoordChipToCell(psPlane* cellCoord, + const psPlane* chipCoord, + const psCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(chipCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->fromChip, NULL); + + cellCoord = psPlaneTransformApply(cellCoord, cell->fromChip, chipCoord); + return(cellCoord); +} + +psPlane* psCoordSkyToCell(psPlane* cellCoord, + const psSphere* skyCoord, + double color, + double magnitude, + const psCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(skyCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent->grommit, 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); +} + +psPlane* psCoordSkyToCellQuick(psPlane* cellCoord, + const psSphere* skyCoord, + const psCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(skyCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL); + PS_ASSERT_PTR_NON_NULL(cell->toTP, NULL); + if (cell->toSky) { + // XXX: Should we use toTP or toSky? + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psCoordSkyToCellQuick(): The cell->toSky transform is ignored. The cell->toTP transform is being used."); + } + + psPlane *tpCoord = NULL; + psChip *whichChip = cell->parent; + psFPA *whichFPA = whichChip->parent; + psProjectionType oldProjectionType; + psPlaneTransform *TPtoCell = NULL; + + // Save the old projection type and set the new projection type to TAN. + oldProjectionType = whichFPA->projection->type; + whichFPA->projection->type = PS_PROJ_TAN; + + if (cellCoord == NULL) { + cellCoord = (psPlane* ) psAlloc(sizeof(psPlane)); + } + + tpCoord = psProject(skyCoord, whichFPA->projection); + + // generate an error if cell->toTP is not linear. + if (0 == p_psIsProjectionLinear(cell->toTP)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM, + "cell to tangent plane"); + } + + TPtoCell = p_psPlaneTransformLinearInvert(cell->toTP); + cellCoord = psPlaneTransformApply(cellCoord, TPtoCell, tpCoord); + + // Restore old projection type. Free memory. + whichFPA->projection->type = oldProjectionType; + psFree(tpCoord); + return (cellCoord); +} + + + Index: /trunk/psModules/src/psAstrometry.h =================================================================== --- /trunk/psModules/src/psAstrometry.h (revision 4577) +++ /trunk/psModules/src/psAstrometry.h (revision 4577) @@ -0,0 +1,534 @@ +/** @file psAstrometry.h +* +* @brief This file defines the basic types for astronomical coordinate +* transformation +* +* @ingroup AstroImage +* +* @author GLG, MHPCC +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-18 18:48:29 $ +* +* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifndef PS_ASTROMETRY_H +#define PS_ASTROMETRY_H + +#include "psType.h" +#include "psImage.h" +#include "psArray.h" +#include "psList.h" +#include "psFunctions.h" +#include "psMetadata.h" +#include "psCoord.h" +#include "psPhotometry.h" + +struct psCell; +struct psChip; +struct psFPA; +struct psExposure; + +/// @addtogroup AstroImage +/// @{ + +/** Wallace's Grommit + * + * SLALib requires several elements to perform the transformations between + * the tangent plane and the sky. Pre-computing these quantities for each + * exposure means that subsequent transformations are faster. For historical + * reasons, this structure is known colloquially as "Wallace's Grommit". + * + */ +typedef struct +{ + const double latitude; ///< geodetic latitude (radians) + const double longitude; ///< longitude + ... (radians) + const double height; ///< height (HM) + const double abberationMag; ///< magnitude of diurnal aberration vector + const double temperature; ///< ambient temperature (TDK) + const double pressure; ///< pressure (PMB) + const double humidity; ///< relative humidity (RH) + const double wavelength; ///< wavelength (WL) + const double lapseRate; ///< lapse rate (TLR) + const double refractA, refractB; ///< refraction constants A and B (radians) + const double siderealTime; ///< local apparent sidereal time (radians) +} +psGrommit; + +/** Fixed Pattern Corrections + * + * The fixed pattern is a correction to the general astrometric solution + * formed by summing the residuals from many observations. The intent is to + * correct for higher-order distortions in the camera system on a coarse + * grid (larger than individual pixels, but smaller than a single cell). + * Hence, in addition to the offsets, we need to specify the size and scale + * of the grid in x and y as well as the origin of the grid. + */ +typedef struct +{ + psS32 nX; ///< Number of elements in x direction + psS32 nY; ///< Number of elements in y direction + double x0; ///< X Position of 0,0 corner on focal plane + double y0; ///< Y Position of 0,0 corner on focal plane + double xScale; ///< Scale of the grid in x direction + double yScale; ///< Scale of the grid in x direction + /// XXX: I added the following memvers to facilitate the psFreeing of the x,y data structures. + psS32 p_ps_xRows; ///< Number of rows in the x member + psS32 p_ps_xCols; ///< Number of cols in the x member + psS32 p_ps_yRows; ///< Number of rows in the y member + psS32 p_ps_yCols; ///< Number of cols in the y member + double **x; ///< The grid of offsets in x + double **y; ///< The grid of offsets in y +} +psFixedPattern; + +/** Readout data structure. + * + * A readout is the result of a single read of a cell (or a portion thereof). + * It contains a pointer to the pixel data, and additional pointers to the + * objects found in the readout, and the readout metadata. It also contains + * the offset from the lower-left corner of the chip, in the case that the + * CCD was windowed. + * + */ +typedef struct +{ + const psS32 col0; ///< Offset from the left of chip. + const psS32 row0; ///< Offset from the bottom of chip. + + const int colParity; ///< Column Readout Direction + const int rowParity; ///< Row Readout Direction + + const psU32 colBins; ///< Amount of binning in x-dimension + const psU32 rowBins; ///< Amount of binning in y-dimension + + psImage* image; ///< Imaging area of readout + psImage* mask; ///< Mask area for readout + psList* objects; ///< Objects derived from Readout + psMetadata* metadata; ///< Readout-level metadata +} +psReadout; + +/** Cell data structure + * + * A cell consists of one or more readouts. It also contains a pointer to the + * cell's metadata, and its parent chip. On the astrometry side, it also + * contains coordinate transforms from the cell to chip, from the cell to + * focal-plane, as well as a "quick and dirty" tranform from the cell to + * sky coordinates. + * + */ +typedef struct +{ + const psS32 col0; ///< Offset from the left of chip + const psS32 row0; ///< Offset from the bottom of chip + + psArray* readouts; ///< readouts from the cell + + psMetadata* metadata; ///< cell-level metadata + + psPlaneTransform* toChip; ///< transformations from cell to chip coordinates + psPlaneTransform* fromChip; ///< transformations from cell to chip coordinates + psPlaneTransform* toFPA; ///< transformations from cell to FPA coordinates + psPlaneTransform* toTP; ///< transformations from cell to FPA coordinates + psPlaneTransform* toSky; ///< transformations from cell to tangent plane coordinates + + struct psChip* parent; ///< chip in which contains this cell +} +psCell; + +/** Chip data structure + * + * A chip consists of one or more cells (according to the number of amplifiers + * on the CCD). It contains a pointer to the chip's metadata, and a pointer + * to the parent focal plane. For astrometry, it contains a coordinate + * transform from the chip to the focal plane, and vis-versa. + * + */ +typedef struct psChip +{ + const psS32 col0; ///< Offset from the left of FPA + const psS32 row0; ///< Offset from the bottom of FPA + + psArray* cells; ///< cells in the chip + + psMetadata* metadata; ///< chip-level metadata + + psPlaneTransform* toFPA; ///< transformation from chip to FPA coordinates + psPlaneTransform* fromFPA; ///< transformation from FPA to chip coordinates + + struct psFPA* parent; ///< FPA which contains this chip +} +psChip; + +/** A Focal-Plane + * + * A focal plane consists of one or more chips (according to the number of + * contiguous silicon). It contains pointers to the focal-plane's metadata + * and the exposure information. For astrometry, it contains a transformation + * from the focal plane to the tangent plane and the fixed pattern residuals. + * Since colors are involved in the transformation, it is necessary to specify + * the color the transformation is defined. We also include some values to + * characterize the quality of the transformation: the root square deviation + * for the x and y transformation fits, and the chi-squared for the + * transformation fit. + * + */ +typedef struct psFPA +{ + psArray* chips; ///< chips in the focal plane array + psMetadata* metadata; ///< focal-plane's metadata + + psPlaneDistort* fromTangentPlane; ///< transformation from tangent plane to focal plane + psPlaneDistort* toTangentPlane; ///< transformation from focal plane to tangent plane + psFixedPattern* pattern; ///< fixed pattern residual offsets + + const struct psExposure* exposure; ///< information about this exposure + psGrommit *grommit; ///< Wallace's grommit + + psPhotSystem* colorPlus; ///< Color reference + psPhotSystem* colorMinus; ///< Color reference + psProjection *projection; ///< projection + + float rmsX; ///< RMS for x transformation fits + float rmsY; ///< RMS for y transformation fits + float chi2; ///< chi^2 of astrometric solution +} +psFPA; + +/** Observatory Information + * + * A container for the observatory data that doesn't change per exposure. + * + */ +typedef struct +{ + const char* name; ///< Name of observatory + const double latitude; ///< Latitude of observatory, east positive (degrees?) + const double longitude; ///< Longitude of observatory (degrees?) + const double height; ///< Height of observatory in meters + const double tlr; ///< Tropospheric Lapse Rate +} +psObservatory; + +/** Exposure Information + * + * Several quantities from the telescope in order to make a first guess at + * the astrometric solution. From these quantities, further quantities can + * be derivedand stored for later use. + * + */ +typedef struct psExposure +{ + const double ra; ///< Telescope boresight, right ascention + const double dec; ///< Telescope boresight, declination + const double hourAngle; ///< Hour angle + const double zenithDistance; ///< Zenith distance + const double azimuth; ///< Azimuth + const psTime* time; ///< Time of observation + const float rotAngle; ///< Rotator position angle in degrees? XXX: see bug#209 + const float temperature; ///< Air temperature in Kelvin + const float pressure; ///< Air pressure in mB + const float humidity; ///< Relative humidity, for refraction + const float exposureTime; ///< Exposure time + const float wavelength; ///< Wavelength in microns + const psObservatory* observatory; ///< Observatory data + + /* Derived quantities */ + const double lst; ///< Local Sidereal Time + const float positionAngle; ///< Position angle + const float parallacticAngle; ///< Parallactic angle + const float airmass; ///< Airmass, calculated from zenith distance + const float parallacticFactor; ///< Parallactic factor + const char* cameraName; ///< name of camera which provided exposure + const char* telescopeName; ///< name of telescope which provided exposure +} +psExposure; + +/** Allocator for psFixedPattern struct + * + * Allocates a new psFixedPattern struct with the attributes coorsponding + * to the parameters set to the said input values. + * + * @return psFixedPattern* New psFixedPattern struct. + */ +psFixedPattern* psFixedPatternAlloc( + double x0, ///< X Position of 0,0 corner on focal plane + double y0, ///< Y Position of 0,0 corner on focal plane + double xScale, ///< Scale of the grid in x direction + double yScale, ///< Scale of the grid in x direction + const psImage *x, ///< The grid of offsets in x + const psImage *y ///< The grid of offsets in y +); + + +/** Allocator for psExposure + * + * We need several quantities from the telescope in order to make a first + * guess at the astrometric solution. From these quantities, further + * quantities can be derived and stored for later use. + * + * @return psExposure* New psExposure struct + */ +psExposure* psExposureAlloc( + double ra, ///< Telescope boresight, right ascention + double dec, ///< Telescope boresight, declination + double hourAngle, ///< Hour angle + double zenithDistance, ///< Zenith distance + double azimuth, ///< Azimuth + const psTime* time, ///< time of observation + float rotAngle, ///< Rotator position angle + float temperature, ///< Temperature + float pressure, ///< Pressure + float humidity, ///< Relative humidity + float exposureTime, ///< Exposure time + float wavelength, ///< wavelength + const psObservatory* observatory ///< Observatory data +); + +/** Allocator for psObservatory + * + * This function shall construct a new psObservatory with attributes + * cooresponding to the function parameters. + * + * @return psObservatory* new psObservatory struct + */ +psObservatory* psObservatoryAlloc( + const char* name, ///< Name of observatory + double latitude, ///< Latitude of observatory, east positive + double longitude, ///< Longitude of observatory + double height, ///< Height of observatory + double tlr ///< Tropospheric Lapse Rate +); + +/** Allocator for psFPA + * + * This function shall make an empty psFPA, with the nChips allocated + * pointers to psChips being set to NULL; all other pointers in the structure + * shall be initialized to NULL, apart from the grommit, which shall be + * constructed on the basis of the exp parameter. + * + * @return psFPA* a newly allocated psFPA + */ +psFPA* psFPAAlloc( + psS32 nChips, ///< number of chips in the FPA + const psExposure* exp ///< the exposure information +); + +/** Allocates a psChip + * + * This allocator shall make an empty psChip, with the nCells allocated + * pointers to psCells being set to NULL; all other pointers in the structure + * shall be initialized to NULL. + * + * @return psChip* newly allocated psChip + */ +psChip* psChipAlloc( + psS32 nCells, ///< number of cells in Chip + psFPA* parentFPA ///< parent FPA +); + +/** Allocates a psCell + * + * The constructor shall make an empty psCell, with the nReadouts allocated + * pointers to psReadouts being set to NULL; all other pointers in the + * structure shall be initialized to NULL. + * + * @return psCell* newly allocated psCell + */ +psCell* psCellAlloc( + psS32 nReadouts, ///< number of readouts in cell + psChip* parentChip ///< parent Chip +); + +/** Allocates a psReadout + * + * All pointers in the structure other than the image shall be initialized + * to NULL. + * + * @return psReadout* newly allocated psReadout with all internal pointers set to NULL + */ +psReadout* psReadoutAlloc(); + +/** Allocates a Wallace's Grommit structure. + * + * The psGrommit is calculated from telescope information for the particular + * exposure. + * + * @return psGrommit* New grommit structure. + */ +psGrommit* psGrommitAlloc( + const psExposure* exp ///< the cooresponding exposure structure. +); + +/** Find cooresponding cell for given FPA coordinate + * + * @return psCell* the cell cooresponding to the coord in FPA + */ +psCell* psCellInFPA( + const psPlane* coord, ///< the coordinate in FPA plane + const psFPA* FPA ///< the FPA to search for the cell +); + +/** Find cooresponding chip for given FPA coordinate + * + * @return psChip* the chip cooresponding to coord + */ +psChip* psChipInFPA( + const psPlane* coord, ///< the coordinate in FPA plane + const psFPA* FPA ///< the FPA to search for the cell +); + +/** Find cooresponding cell for given Chip coordinate + * + * @return psCell* the cell cooresponding to coord + */ +psCell* psCellInChip( + const psPlane* coord, ///< the coordinate in Chip plane + const psChip* chip ///< the chip to search for the cell +); + +/** Translate a cell coordinate into a chip coordinate + * + * @return psPlane* the resulting chip coordinate + */ +psPlane* psCoordCellToChip( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within Cell + const psCell* cell ///< the Cell in interest +); + +/** Translate a chip coordinate into a FPA coordinate + * + * @return psPlane* the resulting FPA coordinate + */ +psPlane* psCoordChipToFPA( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within Chip + const psChip* chip ///< the chip in interest +); + +/** Translate a FPA coordinate into a Tangent Plane coordinate + * + * @return psPlane* the resulting Tangent Plane coordinate + */ +psPlane* psCoordFPAToTP( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within FPA + double color, ///< Color of source + double magnitude, ///< Magnitude of source + const psFPA* fpa ///< the FPA in interest +); + +/** Translate a Tangent Plane coordinate into a Sky coordinate + * + * @return psSphere* the resulting Sky coordinate + */ +psSphere* psCoordTPToSky( + psSphere* out, ///< a sphere struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within Tangent Plane + const psGrommit* grommit ///< the grommit of the tangent plane +); + +/** Translate a cell coordinate into a FPA coordinate + * + * @return psPlane* the resulting FPA coordinate + */ +psPlane* psCoordCellToFPA( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within cell + const psCell* cell ///< the cell in interest +); + +/** Translate a cell coordinate into a Sky coordinate + * + * @return psSphere* the resulting Sky coordinate + */ +psSphere* psCoordCellToSky( + psSphere* out, ///< a sphere struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within cell + double color, ///< Color of source + double magnitude, ///< Magnitude of source + const psCell* cell ///< the cell in interest +); + +/** Translate a cell coordinate into a Sky coordinate using a 'quick and + * dirty' method + * + * @return psSphere* the resulting Sky coordinate + */ +psSphere* psCoordCellToSkyQuick( + psSphere* out, ///< a sphere struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within cell + const psCell* cell ///< the cell in interest +); + +/** Translate a Sky coordinate into a Tangent Plane coordinate + * + * @return psPlane* the resulting Tangent Plane coordinate + */ +psPlane* psCoordSkyToTP( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psSphere* in, ///< the sky coordinate + const psGrommit* grommit ///< the grommit +); + +/** Translate a Tangent Plane coordinate into a FPA coordinate + * + * @return psPlane* the resulting FPA coordinate + */ +psPlane* psCoordTPToFPA( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the coordinate within tangent plane + double color, ///< Color of source + double magnitude, ///< Magnitude of source + const psFPA* fpa ///< the FPA of interest +); + +/** Translate a FPA coordinate into a chip coordinate + * + * @return psPlane* the resulting chip coordinate + */ +psPlane* psCoordFPAToChip( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the FPA coordinate + const psChip* chip ///< the chip of interest +); + +/** Translate a chip coordinate into a cell coordinate + * + * @return psPlane* the resulting cell coordinate + */ +psPlane* psCoordChipToCell( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the Chip coordinate + const psCell* cell ///< the cell of interest +); + +/** Translate a sky coordinate into a cell coordinate + * + * @return psPlane* the resulting cell coordinate + */ +psPlane* psCoordSkyToCell( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psSphere* in, ///< the Sky coordinate + double color, ///< Color of source + double magnitude, ///< Magnitude of source + const psCell* cell ///< the cell of interest +); + +/** Translate a sky coordinate into a cell coordinate using a 'quick and + * dirty' method + * + * @return psPlane* the resulting cell coordinate + */ +psPlane* psCoordSkyToCellQuick( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psSphere* in, ///< the Sky coordinate + const psCell* cell ///< the cell of interest +); + + +#endif // #ifndef PS_ASTROMETRY_H