Index: trunk/psLib/src/astronomy/.cvsignore =================================================================== --- trunk/psLib/src/astronomy/.cvsignore (revision 4540) +++ (revision ) @@ -1,7 +1,0 @@ -Makefile.in -.deps -.libs -Makefile -*.lo -*.la - Index: trunk/psLib/src/astronomy/Makefile.am =================================================================== --- trunk/psLib/src/astronomy/Makefile.am (revision 4540) +++ (revision ) @@ -1,32 +1,0 @@ -#Makefile for astronomy functions of psLib -# -AM_CFLAGS=$(CFLAGS) -DPS_CONFIG_FILE_DEFAULT=\"$(sysconfdir)/pslib/psTime.config\" - -INCLUDES = \ - -I$(top_srcdir)/src/collections \ - -I$(top_srcdir)/src/dataManip \ - -I$(top_srcdir)/src/dataIO \ - -I$(top_srcdir)/src/image \ - -I$(top_srcdir)/src/sysUtils \ - $(all_includes) - -noinst_LTLIBRARIES = libpslibastronomy.la -libpslibastronomy_la_SOURCES = \ - psTime.c \ - psCoord.c \ - psAstrometry.c - -BUILT_SOURCES = psAstronomyErrors.h - -EXTRA_DIST = psAstronomyErrors.dat psAstronomyErrors.h astronomy.i - -psAstronomyErrors.h: psAstronomyErrors.dat - $(top_srcdir)/src/psParseErrorCodes --data=$? $@ - -pslibincludedir = $(includedir) -pslibinclude_HEADERS = \ - psTime.h \ - psCoord.h \ - psAstrometry.h \ - psPhotometry.h - Index: trunk/psLib/src/astronomy/astronomy.i =================================================================== --- trunk/psLib/src/astronomy/astronomy.i (revision 4540) +++ (revision ) @@ -1,23 +1,0 @@ -/* astronomy headers */ -%include "psAstrometry.h" -%include "psAstronomyErrors.h" -%include "psCoord.h" - -%include "psMetadata.h" -%extend psMetadataItem { - const char *get_STR(void) { - if (self->type != PS_META_STR) { - return NULL; - } else { - return self->data.V; - } - } -} - -%apply unsigned int *OUTPUT { unsigned int *nFail }; /* for psMetadataParseConfig */ -%include "psMetadataIO.h" -%clear psU32 *nFail; - -%include "psPhotometry.h" -%include "psTime.h" - Index: trunk/psLib/src/astronomy/psAstrometry.c =================================================================== --- trunk/psLib/src/astronomy/psAstrometry.c (revision 4540) +++ (revision ) @@ -1,826 +1,0 @@ -/** @file psAstrometry.c - * - * @brief This file defines the basic types for astronomical coordinate - * transformation - * - * @ingroup AstroImage - * - * @author GLG, MHPCC - * - * @version $Revision: 1.70 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-25 02:02:04 $ - * - * 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/psLib/src/astronomy/psAstrometry.h =================================================================== --- trunk/psLib/src/astronomy/psAstrometry.h (revision 4540) +++ (revision ) @@ -1,534 +1,0 @@ -/** @file psAstrometry.h -* -* @brief This file defines the basic types for astronomical coordinate -* transformation -* -* @ingroup AstroImage -* -* @author GLG, MHPCC -* -* @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-08 23:40:45 $ -* -* 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 Index: trunk/psLib/src/astronomy/psAstronomyErrors.dat =================================================================== --- trunk/psLib/src/astronomy/psAstronomyErrors.dat (revision 4540) +++ (revision ) @@ -1,69 +1,0 @@ -# -# This file is used to generate psAstronomyErrors.h content -# -# Format is: -# ERRORNAME(one word) ERROR_TEXT -# -# N.B. in code, the ERRORNAME appears as PS_ERRORTEXT_ERRORNAME -#################################################################### -psTime_FILE_NOT_FOUND Failed to open file %s. -psTime_FILE_TOO_MANY_ROWS Too many rows found in file %s. Max number of rows allowed is %d. -psTime_TIME_POSTDATES_TABLE Specified psTime postdates (%g) the table of %s information. -psTime_TIME_PREDATES_TABLE Specified psTime predates (%g) the table of %s information. -psTime_TIME_POSTDATES_TABLES Specified psTime postdates (%g) all tables of %s information. -psTime_TIME_PREDATES_TABLES Specified psTime predates (%g) all tables of %s information. -psTime_TABLE_DUPLICATE_ROWS The %s table was found to have two rows of the same time value. -psTime_TYPE_UNKNOWN Specified type, %d, is not supported. -psTime_TYPE_INCORRECT Specified type, %d, is incorrect. -psTime_TYPE_MISMATCH Specified psTime parameters must have same type. -psTime_GET_TOD_FAILED Failed to determine the current time from gettimeofday function. -psTime_CONVERT_TIME_TO_STRING_FAILED Failed to convert a time via strftime function. -psTime_APPEND_MSEC_FAILED Failed to append millisecond to time string with snprintf function. -psTime_USEC_INVALID The psTime usec attribute value, %u, is invalid. Must be less than 1e6. -psTime_ISOTIME_MALFORMED Specified ISO Time string, '%s', is malformed. Must be in 'YYYY-MM-DDThh:mm:ss.sss' format. -psTime_INTERPOLATION_FAILED Failed time table interpolation. -psTime_INTERPOLATION_FAILED_NAME Failed time table interpolation for '%s'. -psTime_LOOKUP_METADATA_FAILED Failed find '%s' in time metadata. -psTime_BAD_TABLE_COUNT Incorrect number of table files entered. Found: %d. Expected: %d. -psTime_BAD_VECTOR Incorrect vector size. Size: %d, Expected %d. -# -psCoord_PROJECTION_TYPE_UNDEFINED The projection type, %s, is undefined. -psCoord_PROJECTION_TYPE_UNKNOWN The projection type, %d, is unknown. -psCoord_UNITS_UNKNOWN Specified units, 0x%x, is not supported. -psCoord_OFFSET_MODE_UNKNOWN Specified offset mode, 0x%x, is not supported. -psCoord_INVALID_MJD Specified time is less than 1900. -# -psAstrometry_NONLINEAR_TRANSFORM The %s transfrom is not linear. Only linear transforms are supported. -# -psMetadata_METATYPE_INVALID Specified psMetadataType, %d, is not supported. -psMetadata_FORMAT_INVALID Specified print format, %%%c, is not supported. -psMetadata_METATYPE_MISMATCH Specified psMetadataType, %d, is incorrect. Expected %d. -psMetadata_ADD_LIST_FAILED Failed to add metadata item, %s, to items list. -psMetadata_ADD_TABLE_FAILED Failed to add metadata item, %s, to items table. -psMetadata_REMOVE_LIST_FAILED Failed to remove metadata item, %s, from metadata list. -psMetadata_REMOVE_LIST_INDEX_FAILED Failed to remove metadata item, at index %d, from metadata list. -psMetadata_REMOVE_TABLE_FAILED Failed to remove metadata item, %s, from metadata table. -psMetadata_ADD_COLLECTION_FAILED Failed to add metadata item, %s, to metadata collection list. -psMetadata_ADD_FAILED Failed to add metadata item to metadata collection list. -psMetadata_FIND_FAILED Could not find metadata item, %s. -psMetadata_FIND_INDEX_FAILED Could not find metadata item at index %d. -psMetadata_DUPLICATE_NOT_ALLOWED Duplicate metadata item name is not allowed. Use a psMetadataFlags option to allow such action. -psMetadata_REGEX_INVALID Specified regular expression is invalid. %s. -psMetadata_LOCATION_INVALID Specified location, %d, is invalid. -# -psMetadataIO_TYPE_INVALID Specified type, %d, is not supported. -psMetadataIO_EXTNUM_NOTPOSITIVE Specified extension number, %d, is invalid. Value must be positive if no extension name is given. -psMetadataIO_FITS_METATYPE_INVALID Specified FITS metadata type, %c, is not supported. -psMetadataIO_ADD_FAILED Failed to add metadata item, %s. -psMetadataIO_FILE_OPEN_FAILED Failed to open file '%s'. Check if it exists and it has the proper permissions. -psMetadataIO_FILE_MULTIPLE_CHAR More than one '%c' character not allowed. Found on line %u of %s. -psMetadataIO_FILE_ELEMENT_NULL Failed to read a metadata %s on line %u of %s. -psMetadataIO_FILE_TYPE_INVALID Metadata type '%s', found on line %u of %s, is invalid. -psMetadataIO_OVERWRITE_ITEM Duplicate Metadata item, %s, found on line %u of %s. Overwrite not allowed. -psMetadataIO_PARSE_FAILED Failed to parse the value '%s' of metadata item %s, type %s, on line %u of %s. -psMetadataIO_NO_NAME Failed to find key 'name' in table on line %u of %s. -psMetadataIO_TYPE_INVALID_LINE_FILE Specified type, %s, is not supported on line %u of %s. -psMetadataIO_TAG_MISMATCH Start tag, %s and end tag, %s do not agree. -psMetadataIO_TAG_UNKNOWN Invalid end tag name, %s. -psMetadataIO_TYPE_DUPLICATE Specified type, %s, on line %u of %s is already defined. -psMetadataIO_DUPLICATE_MULTI Duplicate MULTI specifier on line %u of %s. Index: trunk/psLib/src/astronomy/psAstronomyErrors.h =================================================================== --- trunk/psLib/src/astronomy/psAstronomyErrors.h (revision 4540) +++ (revision ) @@ -1,89 +1,0 @@ -/** @file psAstronomyErrors.h - * - * @brief Contains the error text for the astronomy functions - * - * @ingroup ErrorHandling - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-08 23:40:45 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PS_ASTRONOMY_ERRORS_H -#define PS_ASTRONOMY_ERRORS_H - -/* N.B., lines between '//~Start' and '//~End' are automatic generated from - * the template following the '//~Start'. The template is used to generate - * the other lines by, for each error text in psAstronomyErrors.dat, the following - * substitutions are made: - * $1 The error text macro name (first word in the psAstronomyErrors.dat lines) - * $2 The error text (rest of the line in psAstronomyErrors.dat) - * $n The order of the source line in psAstronomyErrors.dat (comments excluded) - * - * DO NOT EDIT THE LINES BETWEEN //~Start and //~End! ANY CHANGES WILL BE OVERWRITTEN. - */ - -//~Start #define PS_ERRORTEXT_$1 "$2" -#define PS_ERRORTEXT_psTime_FILE_NOT_FOUND "Failed to open file %s." -#define PS_ERRORTEXT_psTime_FILE_TOO_MANY_ROWS "Too many rows found in file %s. Max number of rows allowed is %d." -#define PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLE "Specified psTime postdates (%g) the table of %s information." -#define PS_ERRORTEXT_psTime_TIME_PREDATES_TABLE "Specified psTime predates (%g) the table of %s information." -#define PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES "Specified psTime postdates (%g) all tables of %s information." -#define PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES "Specified psTime predates (%g) all tables of %s information." -#define PS_ERRORTEXT_psTime_TABLE_DUPLICATE_ROWS "The %s table was found to have two rows of the same time value." -#define PS_ERRORTEXT_psTime_TYPE_UNKNOWN "Specified type, %d, is not supported." -#define PS_ERRORTEXT_psTime_TYPE_INCORRECT "Specified type, %d, is incorrect." -#define PS_ERRORTEXT_psTime_TYPE_MISMATCH "Specified psTime parameters must have same type." -#define PS_ERRORTEXT_psTime_GET_TOD_FAILED "Failed to determine the current time from gettimeofday function." -#define PS_ERRORTEXT_psTime_CONVERT_TIME_TO_STRING_FAILED "Failed to convert a time via strftime function." -#define PS_ERRORTEXT_psTime_APPEND_MSEC_FAILED "Failed to append millisecond to time string with snprintf function." -#define PS_ERRORTEXT_psTime_USEC_INVALID "The psTime usec attribute value, %u, is invalid. Must be less than 1e6." -#define PS_ERRORTEXT_psTime_ISOTIME_MALFORMED "Specified ISO Time string, '%s', is malformed. Must be in 'YYYY-MM-DDThh:mm:ss.sss' format." -#define PS_ERRORTEXT_psTime_INTERPOLATION_FAILED "Failed time table interpolation." -#define PS_ERRORTEXT_psTime_INTERPOLATION_FAILED_NAME "Failed time table interpolation for '%s'." -#define PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED "Failed find '%s' in time metadata." -#define PS_ERRORTEXT_psTime_BAD_TABLE_COUNT "Incorrect number of table files entered. Found: %d. Expected: %d." -#define PS_ERRORTEXT_psTime_BAD_VECTOR "Incorrect vector size. Size: %d, Expected %d." -#define PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED "The projection type, %s, is undefined." -#define PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN "The projection type, %d, is unknown." -#define PS_ERRORTEXT_psCoord_UNITS_UNKNOWN "Specified units, 0x%x, is not supported." -#define PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN "Specified offset mode, 0x%x, is not supported." -#define PS_ERRORTEXT_psCoord_INVALID_MJD "Specified time is less than 1900." -#define PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM "The %s transfrom is not linear. Only linear transforms are supported." -#define PS_ERRORTEXT_psMetadata_METATYPE_INVALID "Specified psMetadataType, %d, is not supported." -#define PS_ERRORTEXT_psMetadata_FORMAT_INVALID "Specified print format, %%%c, is not supported." -#define PS_ERRORTEXT_psMetadata_METATYPE_MISMATCH "Specified psMetadataType, %d, is incorrect. Expected %d." -#define PS_ERRORTEXT_psMetadata_ADD_LIST_FAILED "Failed to add metadata item, %s, to items list." -#define PS_ERRORTEXT_psMetadata_ADD_TABLE_FAILED "Failed to add metadata item, %s, to items table." -#define PS_ERRORTEXT_psMetadata_REMOVE_LIST_FAILED "Failed to remove metadata item, %s, from metadata list." -#define PS_ERRORTEXT_psMetadata_REMOVE_LIST_INDEX_FAILED "Failed to remove metadata item, at index %d, from metadata list." -#define PS_ERRORTEXT_psMetadata_REMOVE_TABLE_FAILED "Failed to remove metadata item, %s, from metadata table." -#define PS_ERRORTEXT_psMetadata_ADD_COLLECTION_FAILED "Failed to add metadata item, %s, to metadata collection list." -#define PS_ERRORTEXT_psMetadata_ADD_FAILED "Failed to add metadata item to metadata collection list." -#define PS_ERRORTEXT_psMetadata_FIND_FAILED "Could not find metadata item, %s." -#define PS_ERRORTEXT_psMetadata_FIND_INDEX_FAILED "Could not find metadata item at index %d." -#define PS_ERRORTEXT_psMetadata_DUPLICATE_NOT_ALLOWED "Duplicate metadata item name is not allowed. Use a psMetadataFlags option to allow such action." -#define PS_ERRORTEXT_psMetadata_REGEX_INVALID "Specified regular expression is invalid. %s." -#define PS_ERRORTEXT_psMetadata_LOCATION_INVALID "Specified location, %d, is invalid." -#define PS_ERRORTEXT_psMetadataIO_TYPE_INVALID "Specified type, %d, is not supported." -#define PS_ERRORTEXT_psMetadataIO_EXTNUM_NOTPOSITIVE "Specified extension number, %d, is invalid. Value must be positive if no extension name is given." -#define PS_ERRORTEXT_psMetadataIO_FITS_METATYPE_INVALID "Specified FITS metadata type, %c, is not supported." -#define PS_ERRORTEXT_psMetadataIO_ADD_FAILED "Failed to add metadata item, %s." -#define PS_ERRORTEXT_psMetadataIO_FILE_OPEN_FAILED "Failed to open file '%s'. Check if it exists and it has the proper permissions." -#define PS_ERRORTEXT_psMetadataIO_FILE_MULTIPLE_CHAR "More than one '%c' character not allowed. Found on line %u of %s." -#define PS_ERRORTEXT_psMetadataIO_FILE_ELEMENT_NULL "Failed to read a metadata %s on line %u of %s." -#define PS_ERRORTEXT_psMetadataIO_FILE_TYPE_INVALID "Metadata type '%s', found on line %u of %s, is invalid." -#define PS_ERRORTEXT_psMetadataIO_OVERWRITE_ITEM "Duplicate Metadata item, %s, found on line %u of %s. Overwrite not allowed." -#define PS_ERRORTEXT_psMetadataIO_PARSE_FAILED "Failed to parse the value '%s' of metadata item %s, type %s, on line %u of %s." -#define PS_ERRORTEXT_psMetadataIO_NO_NAME "Failed to find key 'name' in table on line %u of %s." -#define PS_ERRORTEXT_psMetadataIO_TYPE_INVALID_LINE_FILE "Specified type, %s, is not supported on line %u of %s." -#define PS_ERRORTEXT_psMetadataIO_TAG_MISMATCH "Start tag, %s and end tag, %s do not agree." -#define PS_ERRORTEXT_psMetadataIO_TAG_UNKNOWN "Invalid end tag name, %s." -#define PS_ERRORTEXT_psMetadataIO_TYPE_DUPLICATE "Specified type, %s, on line %u of %s is already defined." -#define PS_ERRORTEXT_psMetadataIO_DUPLICATE_MULTI "Duplicate MULTI specifier on line %u of %s." -//~End - -#endif // #ifndef PS_ASTRONOMY_ERRORS_H Index: trunk/psLib/src/astronomy/psCoord.c =================================================================== --- trunk/psLib/src/astronomy/psCoord.c (revision 4540) +++ (revision ) @@ -1,1265 +1,0 @@ -/** @file psCoord.c -* -* @brief Contains basic coordinate transformation definitions and operations -* -* This file defines the basic types for astronomical coordinate -* transformation -* -* @ingroup CoordinateTransform -* -* @author GLG, MHPCC -* -* @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-12 19:12:00 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ -/******************************************************************************/ -/* INCLUDE FILES */ -/******************************************************************************/ -#include "psType.h" -#include "psCoord.h" -#include "psMemory.h" -#include "psTime.h" -#include "psConstants.h" -#include "psError.h" -#include "psLogMsg.h" -#include "psAstronomyErrors.h" -#include "psAstrometry.h" -#include "psMatrix.h" -#include -#include - -// Modified Julian Day 01/01/1900 00:00:00 -#define MJD_1900 15021.0 - -// Days in Julian century -#define JULIAN_CENTURY 36525.0 - -static void planeFree(psPlane *p) -{ - // There are non dynamic allocated items -} - -static void sphereFree(psSphere *s) -{ - // There are non dynamic allocated items -} - -static void planeTransformFree(psPlaneTransform *pt) -{ - psFree(pt->x); - psFree(pt->y); -} - -static void planeDistortFree(psPlaneDistort *pt) -{ - psFree(pt->x); - psFree(pt->y); -} - -static void projectionFree(psProjection *p) -{ - // There are no dynamically allocated items -} - -/***************************************************************************** -p_psPlaneTransformLinearInvert(transform): : this is a private function which -simply inverts the supplied psPlaneTransform transform. It assumes that -"transform" is linear. - -This program assumes that the inverse of the following linear equations: - X2 = A + (B * X1) + (C * Y1); - Y2 = D + (E * X1) + (F * Y1); -is - Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B)); - X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C)); -or - X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) + - (X2 * -((F/C) / (E - ((F*B)/C)))) + - (Y2 * (1.0 / (E - ((F*B)/C)))); - Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) + - (X2 * -((E/B) / (F - ((C*E)/B)))) + - (Y2 * (1.0 / (F - ((C*E)/B)))); - -XXX: Since there is now a general psPlaneTransformInvert() function, we -should rename this. - *****************************************************************************/ -psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) -{ - PS_ASSERT_PTR_NON_NULL(transform, 0); - PS_ASSERT_PTR_NON_NULL(transform->x, 0); - PS_ASSERT_PTR_NON_NULL(transform->y, 0); - - psF64 A = 0.0; - psF64 B = 0.0; - psF64 C = 0.0; - psF64 D = 0.0; - psF64 E = 0.0; - psF64 F = 0.0; - - A = transform->x->coeff[0][0]; - if (transform->x->nX >= 2) { - B = transform->x->coeff[1][0]; - } - if (transform->x->nY >= 2) { - C = transform->x->coeff[0][1]; - } - D = transform->y->coeff[0][0]; - if (transform->y->nX >= 2) { - E = transform->y->coeff[1][0]; - } - if (transform->y->nY >= 2) { - F = transform->y->coeff[0][1]; - } - - psPlaneTransform *out = psPlaneTransformAlloc(2, 2); - - /* This is sample code from IfA. It didn't work initially, and I did not - spend any time debugging it. - - psF64 a = transform->x->coeff[1][0]; - psF64 b = transform->x->coeff[0][1]; - psF64 c = transform->y->coeff[1][0]; - psF64 d = transform->y->coeff[0][1]; - psF64 e = transform->x->coeff[0][0]; - psF64 f = transform->y->coeff[0][0]; - - psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant - - // Not entirely sure why this works, but it appears to do so....................................! - out->x->coeff[1][0] = invDet * a; - out->x->coeff[0][1] = - invDet * b; - out->y->coeff[1][0] = - invDet * c; - out->y->coeff[0][1] = invDet * d; - - out->x->coeff[0][0] = - invDet * (d * e + c * f); - out->y->coeff[0][0] = - invDet * (b * e + a * f); - */ - out->x->coeff[0][0] = (-D + ((F*A)/C)) / (E - ((F*B)/C)); - out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C)); - out->x->coeff[0][1] = 1.0 / (E - ((F*B)/C)); - out->y->coeff[0][0] = (-D + ((E*A)/B)) / (F - ((C*E)/B)); - out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B)); - out->y->coeff[0][1] = 1.0 / (F - ((C*E)/B)); - - return(out); -} - -/***************************************************************************** -p_psIsProjectionLinear(): this is a private function which simply determines -if the supplied psPlaneTransform transform is linear: if any of the -cooefficients of order 2 are higher are non-zero, then it is not linear. - *****************************************************************************/ -psS32 p_psIsProjectionLinear(psPlaneTransform *transform) -{ - PS_ASSERT_PTR_NON_NULL(transform, 0); - PS_ASSERT_PTR_NON_NULL(transform->x, 0); - PS_ASSERT_PTR_NON_NULL(transform->y, 0); - - for (psS32 i=0;i<(transform->x->nX);i++) { - for (psS32 j=0;j<(transform->x->nY);j++) { - if (transform->x->coeff[i][j] != 0.0) { - if (!(((i == 0) && (j == 0)) || - ((i == 0) && (j == 1)) || - ((i == 1) && (j == 0)))) { - return(0); - } - } - } - } - - for (psS32 i=0;i<(transform->y->nX);i++) { - for (psS32 j=0;j<(transform->y->nY);j++) { - if (transform->y->coeff[i][j] != 0.0) { - if (!(((i == 0) && (j == 0)) || - ((i == 0) && (j == 1)) || - ((i == 1) && (j == 0)))) { - return(0); - } - } - } - } - - return(1); -} - -// XXX: Must test psPlaneAlloc() and planeFree(). -// XXX: Must rewrite code and tests to use these functions. -psPlane* psPlaneAlloc(void) -{ - psPlane *p = psAlloc(sizeof(psPlane)); - - psMemSetDeallocator(p, (psFreeFunc) planeFree); - return(p); -} - -psSphere* psSphereAlloc(void) -{ - psSphere *s = psAlloc(sizeof(psSphere)); - - psMemSetDeallocator(s, (psFreeFunc) sphereFree); - return(s); -} - -psSphereRot* psSphereRotAlloc(double alphaP, - double deltaP, - double phiP) -{ - psSphereRot* rot = psAlloc(sizeof(psSphereRot)); - - double cosDelta = cos(deltaP); - double halfPhi = phiP / 2.0; - double sinHalfPhi = sin(halfPhi); - - // equations are directly from ADD - double vx = cosDelta*cos(alphaP); - double vy = cosDelta*sin(alphaP); - double vz = sin(deltaP); - - rot->q0 = vx*sinHalfPhi; - rot->q1 = vy*sinHalfPhi; - rot->q2 = vz*sinHalfPhi; - rot->q3 = cos(halfPhi); - - return rot; -} - -psSphereRot* psSphereRotQuat(double q0, - double q1, - double q2, - double q3) -{ - psSphereRot* rot = psAlloc(sizeof(psSphereRot)); - - double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); - rot->q0 = q0 / len; - rot->q1 = q1 / len; - rot->q2 = q2 / len; - rot->q3 = q3 / len; - - return rot; -} - -psPlaneTransform* psPlaneTransformAlloc(int n1, int n2) -{ - PS_ASSERT_INT_NONNEGATIVE(n1, NULL); - PS_ASSERT_INT_NONNEGATIVE(n2, NULL); - - psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform)); - pt->x = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD); - pt->y = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD); - - psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree); - return(pt); -} - -psPlane* psPlaneTransformApply(psPlane* out, - const psPlaneTransform* transform, - const psPlane* coords) -{ - PS_ASSERT_PTR_NON_NULL(transform, NULL); - PS_ASSERT_PTR_NON_NULL(transform->x, NULL); - PS_ASSERT_PTR_NON_NULL(transform->y, NULL); - PS_ASSERT_PTR_NON_NULL(coords, NULL); - - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - - out->x = psDPolynomial2DEval( - transform->x, - coords->x, - coords->y - ); - - out->y = psDPolynomial2DEval( - transform->y, - coords->x, - coords->y - ); - - return (out); -} - -psPlaneDistort* psPlaneDistortAlloc(int n1, int n2, int n3, int n4) -{ - PS_ASSERT_INT_NONNEGATIVE(n1, NULL); - PS_ASSERT_INT_NONNEGATIVE(n2, NULL); - PS_ASSERT_INT_NONNEGATIVE(n3, NULL); - PS_ASSERT_INT_NONNEGATIVE(n4, NULL); - - psPlaneDistort *pt = psAlloc(sizeof(psPlaneDistort)); - pt->x = psDPolynomial4DAlloc(n1, n2, n3, n4, PS_POLYNOMIAL_ORD); - pt->y = psDPolynomial4DAlloc(n1, n2, n3, n4, PS_POLYNOMIAL_ORD); - - psMemSetDeallocator(pt, (psFreeFunc) planeDistortFree); - return(pt); -} - -/****************************************************************************** -This transformation takes into account parameters beyond an objects spatial -coordinates: term3 and term4 (magnitude and color). - *****************************************************************************/ -psPlane* psPlaneDistortApply(psPlane* out, - const psPlaneDistort* distort, - const psPlane* coords, - float mag, - float color) -{ - PS_ASSERT_PTR_NON_NULL(distort, NULL); - PS_ASSERT_PTR_NON_NULL(distort->x, NULL); - PS_ASSERT_PTR_NON_NULL(distort->y, NULL); - PS_ASSERT_PTR_NON_NULL(coords, NULL); - - if (out == NULL) { - out = (psPlane* ) psAlloc(sizeof(psPlane)); - } - out->x = psDPolynomial4DEval( - distort->x, - coords->x, - coords->y, - mag, - color - ); - out->y = psDPolynomial4DEval( - distort->y, - coords->x, - coords->y, - mag, - color - ); - return (out); -} - -/****************************************************************************** -XXX: We convert Right Ascension angles to the range 0:PI. Is that acceptable? -XXX: Should we do something for Declination as well? - *****************************************************************************/ -psSphere* psSphereRotApply(psSphere* out, - const psSphereRot* transform, - const psSphere* coord) -{ - PS_ASSERT_PTR_NON_NULL(transform, NULL); - PS_ASSERT_PTR_NON_NULL(coord, NULL); - - if (out == NULL) { - out = psSphereAlloc(); - } - - - // apply the transform by creating a new psSphereRot from the input coord - // and combining it with the input transform (see ADD) - psSphereRot* coordRot = psSphereRotAlloc(coord->r, coord->d, 0); - coordRot->q3 = 0.0; - coordRot = psSphereRotCombine(coordRot, transform, coordRot); - // N.B., we can recycle coordRot right away due to the implementation of - // psSphereRotCombine puts the values of coordRot in a local variable first - - out->r = atan2(coordRot->q1,coordRot->q0); - out->d = atan2(coordRot->q2,sqrt(coordRot->q1*coordRot->q1+coordRot->q0*coordRot->q0)); - - return out; -} - -psSphereRot* psSphereRotCombine(psSphereRot* out, - const psSphereRot* rot1, - const psSphereRot* rot2) -{ - PS_ASSERT_PTR_NON_NULL(rot1, NULL); - PS_ASSERT_PTR_NON_NULL(rot2, NULL); - - if (out == NULL) { - out = (psSphereRot* ) psAlloc(sizeof(psSphereRot)); - } - - double a0 = rot1->q0; - double a1 = rot1->q1; - double a2 = rot1->q2; - double a3 = rot1->q3; - double b0 = rot2->q0; - double b1 = rot2->q1; - double b2 = rot2->q2; - double b3 = rot2->q3; - - // following came from ADD - out->q0 = b3*a0 + b2*a1 - b1*a2 + b0*a3; - out->q1 = b3*a1 - b2*a0 + b1*a3 + b0*a2; - out->q2 = b3*a2 + b2*a3 + b1*a0 - b0*a1; - out->q3 = b3*a3 - b3*a2 - b1*a1 - b0*a0; - - return out; -} - -psSphereRot *psSphereRotInvert(psSphereRot *rot) -{ - PS_ASSERT_PTR_NON_NULL(rot, NULL); - - double norm = sqrt(rot->q0*rot->q0 + rot->q1*rot->q1 + rot->q2*rot->q2 + rot->q3*rot->q3); - rot->q1 = -rot->q1 / norm; - rot->q2 = -rot->q2 / norm; - rot->q3 = -rot->q3 / norm; - - return rot; -} - -psSphereRot* psSphereRotICRSToEcliptic(const psTime *time) -{ - psF64 T; - - // Check for null parameter - PS_ASSERT_PTR_NON_NULL(time, NULL); - - // Convert psTime to MJD - psF64 MJD = psTimeToMJD(time); - - // Check the specified MJD is greater than 1900 - if ( MJD < MJD_1900 ) { - psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD); - return NULL; - } - - // Calculate number of Julian centuries since 1900 - T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; - - psF64 alphaP = 0.0; - psF64 deltaP = DEG_TO_RAD(23.0) + - MIN_TO_RAD(27.0) + - SEC_TO_RAD(8.26) - - (SEC_TO_RAD(46.845) * T) - - (SEC_TO_RAD(0.0059) * T * T) + - (SEC_TO_RAD(0.00181) * T * T * T); - psF64 phiP = 0.0; - - // Don't neglect the minus sign on deltaP (bug 244): - return (psSphereRotAlloc(alphaP, deltaP, phiP)); -} - - -psSphereRot* psSphereRotEclipticToICRS(const psTime *time) -{ - psF64 T; - - // Check for null parameter - PS_ASSERT_PTR_NON_NULL(time, NULL); - - // Convert psTime to MJD - psF64 MJD = psTimeToMJD(time); - - // Check the specified MJD is greater than 1900 - if ( MJD < MJD_1900 ) { - psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD); - return NULL; - } - - // Calculate number of Julian centuries since 1900 - T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; - - psF64 alphaP = 0.0; - psF64 deltaP = DEG_TO_RAD(23.0) + - MIN_TO_RAD(27.0) + - SEC_TO_RAD(8.26) - - (SEC_TO_RAD(46.845) * T) - - (SEC_TO_RAD(0.0059) * T * T) + - (SEC_TO_RAD(0.00181) * T * T * T); - psF64 phiP = 0.0; - - return (psSphereRotAlloc(alphaP, -deltaP, phiP)); -} - -// XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS() -psSphereRot* psSphereRotGalacticToICRS(void) -{ - psF64 alphaP = DEG_TO_RAD(32.93192); - psF64 deltaP = DEG_TO_RAD(-62.87175); - psF64 phiP = DEG_TO_RAD(282.85948); - - return (psSphereRotAlloc(alphaP, deltaP, phiP)); -} - -psSphereRot* psSphereRotICRSToGalactic(void) -{ - psF64 alphaP = DEG_TO_RAD(282.85948); - psF64 deltaP = DEG_TO_RAD(62.87175); - psF64 phiP = DEG_TO_RAD(32.93192); - - return (psSphereRotAlloc(alphaP, deltaP, phiP)); -} - -psProjection* psProjectionAlloc( - double R, - double D, - double Xs, - double Ys, - psProjectionType type) -{ - psProjection *p = psAlloc(sizeof(psProjection)); - p->D = D; - p->R = R; - p->Xs = Xs; - p->Ys = Ys; - p->type = type; - - psMemSetDeallocator(p, (psFreeFunc) projectionFree); - return(p); -} - -psPlane* psProject(const psSphere* coord, - const psProjection* projection) -{ - PS_ASSERT_PTR_NON_NULL(coord, NULL); - PS_ASSERT_PTR_NON_NULL(projection, NULL); - - psF64 theta = 0.0; - psF64 phi = 0.0; - - // Allocate return value - psPlane* out = psPlaneAlloc(); - - // Convert to projection spherical coordinate system - theta = asin( sin(coord->d)*sin(projection->D) + - cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R)); - phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R), - sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) ); - - // Perform the specified projection - // Gnomonic projection - if (projection->type == PS_PROJ_TAN) { - out->x = (cos(theta)*sin(phi))/sin(theta); - out->y = (-1.0*cos(theta)*cos(phi))/sin(theta); - // Othrographic projection - } else if (projection->type == PS_PROJ_SIN) { - out->x = cos(theta)*sin(phi); - out->y = -1.0*cos(theta)*cos(phi); - // Hammer-Aitoff projection - } else if ( projection->type == PS_PROJ_AIT) { - psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); - out->x = 2.0*zeta*cos(theta)*sin(phi/2.0); - out->y = zeta*sin(theta); - // Parabolic projection - } else if ( projection->type == PS_PROJ_PAR) { - out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); - out->y = M_PI*sin(theta/3.0); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, - projection->type); - psFree(out); - return NULL; - } - - // Apply plate scales - out->x *= projection->Xs; - out->y *= projection->Ys; - - // Return output - return out; -} - -psSphere* psDeproject(const psPlane* coord, - const psProjection* projection) -{ - PS_ASSERT_PTR_NON_NULL(coord, NULL); - PS_ASSERT_PTR_NON_NULL(projection, NULL); - - psF64 theta = 0.0; - psF64 phi = 0.0; - - // Allocate return sphere structure - psSphere* out = psSphereAlloc(); - - // Remove plate scales - psF64 x = coord->x/projection->Xs; - psF64 y = coord->y/projection->Ys; - - // Perform inverse projection - // Gnonomic deprojection - if ( projection->type == PS_PROJ_TAN) { - phi = atan(-1.0*x/y); - theta = atan(1.0/sqrt(x*x+y*y)); - // Orhtographic deprojection - } else if ( projection->type == PS_PROJ_SIN) { - phi = atan((-1.0*x)/y); - theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y)); - // Hammer-Aitoff deprojection - } else if ( projection->type == PS_PROJ_AIT) { - psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0))); - phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) ); - theta = asin(y*z); - // Parabolic deprojection - } else if ( projection->type == PS_PROJ_PAR) { - psF64 rho = y/M_PI; - phi = x/(1.0 - 4.0*rho*rho); - theta = 3.0*asin(rho); - // Invalid deprojection type - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, - projection->type); - psFree(out); - return NULL; - } - - // Convert from projection spherical coordinates - out->d = asin( sin(theta)*sin(projection->D) + - cos(theta)*cos(projection->D)*cos(phi) ); - out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi), - sin(theta)*cos(projection->D) - - cos(theta)*sin(projection->D)*cos(phi) ); - - // Return sphere coordinate - return out; -} - -/****************************************************************************** -The basic idea is to project both positions onto the linear plane, with -position1 at the center, then calculate the linear offset between those -projections. - -XXX: Do I 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: Does PS_LINEAR mode make sense? The result must be returned in psSphere - regardless of the mode. - -XXX: How to compound errors? - *****************************************************************************/ -psSphere* psSphereGetOffset(const psSphere* position1, - const psSphere* position2, - psSphereOffsetMode mode, - psSphereOffsetUnit unit) -{ - PS_ASSERT_PTR_NON_NULL(position1, NULL); - PS_ASSERT_PTR_NON_NULL(position2, NULL); - - // Check positions near 90 degree and issue warnings if necessary - if (position1->d >= DEG_TO_RAD(90.0)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psDeproject(): position1->d is larger than 90 degrees. Returning NULL."); - return NULL; - } - if (position2->d >= DEG_TO_RAD(90.0)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psDeproject(): position2->d is larger than 90 degrees. Returning NULL."); - return NULL; - } - - // Allocate return structure - psSphere* tmp = psSphereAlloc(); - - // Mode is LINEAR - Use first position as projection center and project second point - // onto tangent plane, set point projected into psSphere structure x->r y->d - if (mode == PS_LINEAR) { - psProjection* proj = psProjectionAlloc(position1->r, - position1->d, - 1.0, - 1.0, - PS_PROJ_TAN); - - // Perform projection onto tangent plane - psPlane* lin = psProject(position2, proj); - - // Set return values - tmp->r = lin->x; - tmp->d = lin->y; - - // Free data structures allocated - psFree(proj); - psFree(lin); - - // Mode is SPHERICAL - Get difference between positiion 1 and position 2 and convert - // offset value from radians to desired units and return - } else if (mode == PS_SPHERICAL) { - tmp->r = position2->r - position1->r; - tmp->d = position2->d - position1->d; - - // Wrap these to an acceptable range. This assumes that all - // angles are in radians. - tmp->r = fmod(tmp->r, 2*M_PI); - tmp->d = fmod(tmp->d, 2*M_PI); - tmp->rErr = 0.0; - tmp->dErr = 0.0; - - // Convert to desired units - if (unit == PS_ARCSEC) { - tmp->r = RAD_TO_SEC(tmp->r); - tmp->d = RAD_TO_SEC(tmp->d); - } else if (unit == PS_ARCMIN) { - tmp->r = RAD_TO_MIN(tmp->r); - tmp->d = RAD_TO_MIN(tmp->d); - } else if (unit == PS_DEGREE) { - tmp->r = RAD_TO_DEG(tmp->r); - tmp->d = RAD_TO_DEG(tmp->d); - } else if (unit == PS_RADIAN) {} - else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psCoord_UNITS_UNKNOWN, - unit); - psFree(tmp); - return NULL; - } - // Invalid mode - } else { - - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN, - mode); - psFree(tmp); - return NULL; - } - - // Return value - return tmp; -} - -/****************************************************************************** -XXX: Do we need to check for unacceptable transformation parameters? Maybe, - if the points are on the North/South Pole, etc? - -XXX: Do we need to somehow scale this projection? - -XXX: I copied the algorithm from the ADD exactly. - -XXX: Should we compound errors? - *****************************************************************************/ - -psSphere* psSphereSetOffset(const psSphere* position, - const psSphere* offset, - psSphereOffsetMode mode, - psSphereOffsetUnit unit) -{ - PS_ASSERT_PTR_NON_NULL(position, NULL); - PS_ASSERT_PTR_NON_NULL(offset, NULL); - - psSphere* tmp; - psF64 tmpR = 0.0; - psF64 tmpD = 0.0; - - // If mode is linear then set position to projection center - // and offset to linear coordinate then deproject to obtain - // new sphere coordinate - if (mode == PS_LINEAR) { - - // Allocate plane coordinate and set coordinate - psPlane* lin = psPlaneAlloc(); - lin->x = offset->r; - lin->y = offset->d; - - // Allocate and set projection structure - psProjection* proj = psProjectionAlloc(position->r, - position->d, - 1.0, - 1.0, - PS_PROJ_TAN); - - // Project tangent plane coord to spherical coord - tmp = psDeproject(lin, proj); - - // Free data structures used - psFree(proj); - psFree(lin); - - // If mode is spherical then convert offset to radians, add the offset - // to the position and wrap to 0 to 2pi - } else if (mode == PS_SPHERICAL) { - - // Convert offset unit to radians - if (unit == PS_ARCSEC) { - tmpR = SEC_TO_RAD(offset->r); - tmpD = SEC_TO_RAD(offset->d); - } else if (unit == PS_ARCMIN) { - tmpR = MIN_TO_RAD(offset->r); - tmpD = MIN_TO_RAD(offset->d); - } else if (unit == PS_DEGREE) { - tmpR = DEG_TO_RAD(offset->r); - tmpD = DEG_TO_RAD(offset->d); - } else if (unit == PS_RADIAN) { - tmpR = offset->r; - tmpD = offset->d; - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psCoord_UNITS_UNKNOWN, - unit); - return NULL; - } - - // Allocate sphere structure to return - tmp = psSphereAlloc(); - - // Add offset and wrap to 0 to 2PI if necessary - tmp->r = position->r + tmpR; - tmp->r = fmod(tmp->r, 2.0*M_PI); - tmp->d = position->d + tmpD; - tmp->d = fmod(tmp->d, 2.0*M_PI); - tmp->rErr = 0.0; - tmp->dErr = 0.0; - - // Invalid mode report error - } else { - - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN, - mode); - return NULL; - } - - return tmp; -} - - - -/****************************************************************************** -psSpherePrecess(coords, fromTime, toTime): - -XXX: Use static memory for tmpST. - *****************************************************************************/ -psSphere *psSpherePrecess(psSphere *coords, - const psTime *fromTime, - const psTime *toTime) -{ - // Check input for NULL pointers - PS_ASSERT_PTR_NON_NULL(coords, NULL); - PS_ASSERT_PTR_NON_NULL(fromTime, NULL); - PS_ASSERT_PTR_NON_NULL(toTime, NULL); - - // Calculate Julian centuries - psF64 fromMJD = psTimeToMJD(fromTime); - psF64 toMJD = psTimeToMJD(toTime); - psF64 T = (toMJD - fromMJD) / JULIAN_CENTURY; - - // Calculate conversion constants - psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + - (DEG_TO_RAD(0.0000839) * T * T) + - (DEG_TO_RAD(0.000005) * T * T * T)); - - psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) - - (DEG_TO_RAD(0.0001185) * T * T) - - (DEG_TO_RAD(0.0000116) * T * T * T); - - psF64 phiP = DEG_TO_RAD(90.0) + ((DEG_TO_RAD(0.6406161) * T) + - (DEG_TO_RAD(0.0003041) * T * T) + - (DEG_TO_RAD(0.0000051) * T * T * T)); - - // Create transform with proper constants - psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP); - - // Apply transform to coordinates - psSphere *out = psSphereRotApply(NULL, tmpST, coords); - - psFree(tmpST); - - return(out); -} - -/***************************************************************************** -multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and -multiplies them. Basically, for each non-zero coeff in the trans1 coeff[][] -array, you must multiply by all non-zero coeffs in trans2. - -XXX: Inefficient in that the out polynomial is allocated every time. - *****************************************************************************/ - -static psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1, - psDPolynomial2D *trans2) -{ - //TRACE: printf("multiplyDPoly2D(%d %d: %d %d)\n", trans1->nX, trans1->nY, trans2->nX, trans2->nY); - psS32 orderX = (trans1->nX + trans2->nX) - 1; - psS32 orderY = (trans1->nY + trans2->nY) - 1; - - psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); - //TRACE: printf("Creating poly (%d, %d)\n", orderX, orderY); - for (psS32 i = 0 ; i < out->nX; i++) { - for (psS32 j = 0 ; j < out->nY; j++) { - out->coeff[i][j] = 0.0; - out->mask[i][j] = 0; - } - } - - for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) { - for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) { - if (0.0 != trans1->coeff[t1x][t1y]) { - for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) { - for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) { - /* Possible debug-only macro which checks these coords? - if ((t1x+t2x) >= orderX) - printf("BAD 1\n"); - if ((t1y+t2y) >= orderY) - printf("BAD 2\n"); - */ - out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]); - } - } - } - } - } - return(out); -} - - -/***************************************************************************** -psPlaneTransformCombine(out, trans1, trans2) - -XXX: Much room for optimization. Currently, we call the polyMultiply -routine far too many times. - *****************************************************************************/ -psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out, - const psPlaneTransform *trans1, - const psPlaneTransform *trans2, - psRegion region, - int nSamples) -{ - - // XXX: This does not yet use region and nSamples: need to modify -rdd - - PS_ASSERT_PTR_NON_NULL(trans1, NULL); - PS_ASSERT_PTR_NON_NULL(trans2, NULL); - //TRACE: printf("psPlaneTransformCombine(%d, %d, %d, %d: %d, %d, %d, %d)\n", trans1->x->nX, trans1->x->nY, trans1->y->nX, trans1->y->nY, trans2->x->nX, trans2->x->nY, trans2->y->nX, trans2->y->nY); - // - // Determine the size of the new psPlaneTransform. - // - // PS_MAX( Number of x terms in T2->x * number of x terms in T1->x, - // Number of y terms in T2->x * number of x terms in T1->y, - psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX), - (trans2->x->nY * trans1->y->nX)); - psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY), - (trans2->x->nY * trans1->y->nY)); - - psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX), - (trans2->y->nY * trans1->y->nX)); - psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY), - (trans2->y->nY * trans1->y->nY)); - psS32 orderX = PS_MAX(orderXnX, orderYnX); - psS32 orderY = PS_MAX(orderXnY, orderYnY); - - // - // Allocate the new psPlaneTransform, if necessary. - // - psPlaneTransform *myPT = NULL; - if (out == NULL) { - myPT = psPlaneTransformAlloc(orderX, orderY); - } else { - if ((out->x->nX == orderX) && (out->x->nY == orderY) && - (out->y->nX == orderX) && (out->y->nY == orderY)) { - myPT = out; - } else { - psFree(out); - myPT = psPlaneTransformAlloc(orderX, orderY); - } - } - - // - // Initialize the new psPlaneTransform, if necessary. - // - for (psS32 i = 0 ; i < orderX ; i++) { - for (psS32 j = 0 ; j < orderY ; j++) { - myPT->x->coeff[i][j] = 0.0; - myPT->x->mask[i][j] = 0; - myPT->y->coeff[i][j] = 0.0; - myPT->y->mask[i][j] = 0; - } - } - - // - // For each term (a * x^i * y^j) in trans2, we substitute the appropriate - // equation from trans1, and raise it to the appropriate power. This is - // done via the multiplyDPoly2D(). The result is a polynomial (currPoly) - // and its coefficients are added into the myPT coeff matrix. - // - // XXX: This is horribly inefficient in that the trans1 polys are repeatedly - // multiplied against themselves. This can easily be improved. - // - - for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) { - for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) { - psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); - - currPoly->coeff[0][0] = 1.0; - currPoly->mask[0][0] = 0; - psDPolynomial2D *newPoly = NULL; - - if (trans2->x->mask[t2x][t2y] == 0) { - // Must raise trans1->y to the t2y-power. - for (psS32 c = 0 ; c < t2y; c++) { - newPoly = multiplyDPoly2D(currPoly, trans1->y); - psFree(currPoly); - currPoly = newPoly; - } - - // Must raise trans1->x to the t2x-power. - for (psS32 c = 0 ; c < t2x; c++) { - newPoly = multiplyDPoly2D(currPoly, trans1->x); - psFree(currPoly); - currPoly = newPoly; - } - - // Set the appropriate coeffs in myPT->x - for (psS32 i = 0 ; i < currPoly->nX ; i++) { - for (psS32 j = 0 ; j < currPoly->nY ; j++) { - myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y]; - } - } - } - psFree(currPoly); - } - } - - - for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) { - for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) { - psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); - currPoly->coeff[0][0] = 1.0; - currPoly->mask[0][0] = 0; - psDPolynomial2D *newPoly = NULL; - - if (trans2->y->mask[t2x][t2y] == 0) { - - // Must raise trans1->y to the t2y-power. - for (psS32 c = 0 ; c < t2y; c++) { - newPoly = multiplyDPoly2D(currPoly, trans1->y); - psFree(currPoly); - currPoly = newPoly; - } - - // Must raise trans1->x to the t2x-power. - for (psS32 c = 0 ; c < t2x; c++) { - newPoly = multiplyDPoly2D(currPoly, trans1->x); - psFree(currPoly); - currPoly = newPoly; - } - - // Set the appropriate coeffs in myPT->x - for (psS32 i = 0 ; i < currPoly->nX ; i++) { - for (psS32 j = 0 ; j < currPoly->nY ; j++) { - myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y]; - } - } - } - psFree(currPoly); - } - } - - //TRACE: printf("Exiting combine()\n"); - return(myPT); -} - -/***************************************************************************** -psPlaneTransformFit(trans, source, dest, nRejIter, sigmaClip) - -XXX: What about nRejIter? Iterations? -XXX: Use static vectors for internal data. -XXX: This code has problems with data that corresponds to a non-linear fit. - *****************************************************************************/ -bool psPlaneTransformFit(psPlaneTransform *trans, - const psArray *source, - const psArray *dest, - int nRejIter, - float sigmaClip) -{ - PS_ASSERT_PTR_NON_NULL(trans, NULL); - PS_ASSERT_PTR_NON_NULL(source, NULL); - PS_ASSERT_PTR_NON_NULL(dest, NULL); - - psS32 numCoords = PS_MIN(source->n, dest->n); - psS32 order = PS_MAX(trans->x->nX, trans->x->nY); - order = PS_MAX(order, trans->y->nX); - order = PS_MAX(order, trans->y->nY); - - // - // Create fake polynomial to use in evaluation - // - psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD); - for (int i = 0; i < order; i++) { - for (int j = 0; j < order; j++) { - fakePoly->coeff[i][j] = 1.0; - fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate - } - } - - // - // Initialize the matrix and vectors - // - psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients - psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution - psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x - psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y - for (psS32 i = 0; i < nCoeff; i++) { - for (psS32 j = 0; j < nCoeff; j++) { - matrix->data.F64[i][j] = 0.0; - } - xVector->data.F64[i] = 0.0; - yVector->data.F64[i] = 0.0; - } - - // - // Iterate over the grid points - // - for (psS32 g = 0; g < numCoords; g++) { - // Iterate over the polynomial coefficients, accumulating the matrix and vectors - - for (psS32 i = 0, ijIndex = 0; i < order; i++) { - for (psS32 j = 0; j < order - i; j++, ijIndex++) { - fakePoly->mask[i][j] = 0; - psF64 xIn = ((psPlane *) source->data[g])->x; - psF64 yIn = ((psPlane *) source->data[g])->y; - psF64 xOut = ((psPlane *) dest->data[g])->x; - psF64 yOut = ((psPlane *) dest->data[g])->y; - psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn); - fakePoly->mask[i][j] = 1; - - for (psS32 m = 0, mnIndex = 0; m < order; m++) { - for (psS32 n = 0; n < order - m; n++, mnIndex++) { - fakePoly->mask[m][n] = 0; - psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn); - fakePoly->mask[m][n] = 1; - - matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly; - } - } - - xVector->data.F64[ijIndex] += ijPoly * xOut; - yVector->data.F64[ijIndex] += ijPoly * yOut; - } - } - } - - // - // Solution via LU Decomposition - // - psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition - psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix - psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x - psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y - - // - // XXX: Should check the output of the matrix routines and return false if bad. - // - - // - // Stuff coefficients into transformation - // - for (psS32 i = 0, ijIndex = 0; i < order; i++) { - for (psS32 j = 0; j < order - i; j++, ijIndex++) { - trans->x->coeff[i][j] = xSolution->data.F64[ijIndex]; - trans->y->coeff[i][j] = ySolution->data.F64[ijIndex]; - } - } - - psFree(fakePoly); - psFree(permutation); - psFree(luMatrix); - psFree(xSolution); - psFree(ySolution); - psFree(matrix); - psFree(xVector); - psFree(yVector); - - return(true); -} - - -/***************************************************************************** -psPlaneTransformInvert(out, in, region, nSamples) - -// XXX: Use static data structures. - *****************************************************************************/ -psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, - const psPlaneTransform *in, - psRegion region, - int nSamples) -{ - PS_ASSERT_PTR_NON_NULL(in, NULL); - // - // If the transform is linear, then invert it exactly and return. - // - if (p_psIsProjectionLinear((psPlaneTransform *) in)) { - return(p_psPlaneTransformLinearInvert((psPlaneTransform *) in)); - } - PS_INT_COMPARE(1, nSamples, NULL); - - // Ensure that the input transformation is symmetrical. - if ((in->x->nX != in->x->nY) || - (in->y->nX != in->y->nY) || - (in->x->nX != in->y->nX)) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY."); - } - psS32 order = in->x->nX; - - psPlaneTransform *myPT = NULL; - psPlane *inCoord = psPlaneAlloc(); - psPlane *outCoord = psPlaneAlloc(); - - // - // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size. - // - if (out == NULL) { - myPT = psPlaneTransformAlloc(order, order); - } else { - if ((out->x->nX == order) && (out->x->nY == order) && - (out->y->nX == order) && (out->y->nY == order)) { - myPT = out; - } else { - psFree(out); - myPT = psPlaneTransformAlloc(order, order); - } - } - - // - // Copy the input transform to myPT. - // - for (psS32 i = 0 ; i < in->x->nX ; i++) { - for (psS32 j = 0 ; j < in->x->nY ; j++) { - myPT->x->coeff[i][j] = in->x->coeff[i][j]; - } - } - for (psS32 i = 0 ; i < in->y->nX ; i++) { - for (psS32 j = 0 ; j < in->y->nY ; j++) { - myPT->y->coeff[i][j] = in->y->coeff[i][j]; - } - } - - // - // Create a grid of xin,yin --> xout,yout - // - psArray *inData = psArrayAlloc(nSamples * nSamples); - psArray *outData = psArrayAlloc(nSamples * nSamples); - for (psS32 i = 0 ; i < inData->n; i++) { - inData->data[i] = (psPtr *) psPlaneAlloc(); - outData->data[i] = (psPtr *) psPlaneAlloc(); - } - - // - // Initialize the grid. Since we want the inverse of the transformation, the - // inCoords are written to the outData vector, and the outCoords are written - // to the inData vector. - // - psS32 cnt = 0; - for (int yint = 0; yint < nSamples; yint++) { - inCoord->y = region.y0 + ((psF32) yint) * ((region.y1 - region.y0) / ((psF32) nSamples)); - for (int xint = 0; xint < nSamples; xint++) { - inCoord->x = region.x0 + ((psF32) xint) * ((region.x1 - region.x0) / ((psF32) nSamples)); - (void)psPlaneTransformApply(outCoord, in, inCoord); - ((psPlane *) outData->data[cnt])->x = inCoord->x; - ((psPlane *) outData->data[cnt])->y = inCoord->y; - ((psPlane *) inData->data[cnt])->x = outCoord->x; - ((psPlane *) inData->data[cnt])->y = outCoord->y; - - cnt++; - } - } - // XXX: what values should be used here? - bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0); - - psFree(inCoord); - psFree(outCoord); - psFree(inData); - psFree(outData); - - if (rc == true) { - return(myPT); - } - - // XXX: Generate an error message, or warning message. - return(NULL); -} Index: trunk/psLib/src/astronomy/psCoord.h =================================================================== --- trunk/psLib/src/astronomy/psCoord.h (revision 4540) +++ (revision ) @@ -1,488 +1,0 @@ -/** @file psCoord.h -* -* @brief Contains basic coordinate transformation definitions and operations -* -* This file defines the basic types for astronomical coordinate -* transformation -* -* @ingroup CoordinateTransform -* -* @author GLG, MHPCC -* -* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-12 19:12:00 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ - -#ifndef PS_COORD_H -#define PS_COORD_H - -#include "psType.h" -#include "psImage.h" -#include "psArray.h" -#include "psList.h" -#include "psFunctions.h" -// N.B. inclusion of psTime.h was done to after the typedefs to handle cross-dependency of typedefs - -/// @addtogroup CoordinateTransform -/// @{ - -/** Euclidiean Coordinate System. - * - * Both detector and sky positions will be used extensively in the IPP. One - * coordinate system to be used is linear coordinates which conform to - * Euclidean geometry. - * - */ -typedef struct -{ - double x; ///< x position - double y; ///< y position - double xErr; ///< Error in x position - double yErr; ///< Error in y position -} -psPlane; - -/** Angular Coordinate System - * - * Both detector and sky positions will be used extensively in the IPP. One - * coordinate system to be used is angular coordinates for which additional - * care must often be taken in comparison to a euclidiean coordinate system. - * - */ -typedef struct -{ - double r; ///< RA - double d; ///< Dec - double rErr; ///< Error in RA - double dErr; ///< Error in Dec -} -psSphere; - -/** Cubic Coordinate System - * - */ -typedef struct -{ - double x; ///< cos (DEC) cos (RA) - double y; ///< cos (DEC) sic (RA) - double z; ///< sin (DEC) - double xErr; ///< Error in x - double yErr; ///< Error in y - double zErr; ///< Error in z -} -psCube; - -/** Spherical rotations represent coordinate transformation in 3-D, as well as - * the effects of precession and nutation. The structure contains the - * elements of a quaternion to represent the spherical rotational. - * - */ -typedef struct -{ - double q0; ///< first element of the quaternion - double q1; ///< second element of the quaternion - double q2; ///< third element of the quaternion - double q3; ///< fourth element of the quaternion -} -psSphereRot; - -/** 2D Polynomial Transform - * - * A transform between coordinate systems that consists simply of two 2D - * polynomials to transform both components - the output coordinates depend - * only on the input coordinates and no other quantities of objects at those - * coordinates. - * - */ -typedef struct -{ - psDPolynomial2D* x; ///< 2D polynomial transform of X coordinates - psDPolynomial2D* y; ///< 2D polynomial transform of Y coordinates -} -psPlaneTransform; - -/** 4D Polynomial Transform - * - * A transform between coordinate systems that consists of two 4D polynomials - * in which the output coordinates are also specified to be a function of the - * magnitude and color of the object with the given coordinates. This type of - * coordinate transformation is necessary to represent the (color-dependent) - * optical distortions caused by the atmosphere and camera optics, and the - * possibly effects of charge transfer inefficiency. - * - * The lowest two terms are the x and y axis of the target system. The higher - * two terms may represent magnitude and color terms. - */ -typedef struct -{ - psDPolynomial4D* x; ///< 4D polynomial transform of X coordinates - psDPolynomial4D* y; ///< 4D polynomial transform of Y coordinates -} -psPlaneDistort; - -/** Projection type for projection/deprojection - * - * @see psProject, psDeproject - * - */ -typedef enum { - PS_PROJ_TAN, ///< Tangent projection - PS_PROJ_SIN, ///< Sine projection - PS_PROJ_AIT, ///< Aitoff projection - PS_PROJ_PAR, ///< Par projection - // PS_PROJ_GLS, ///< GLS projection - // PS_PROJ_CAR, ///< CAR projection - // PS_PROJ_MER, ///< MER projection - PS_PROJ_NTYPE ///< Number of types; must be last. -} psProjectionType; - -/** Parameter set for projection/deprojection - * - * @see psProject, psDeproject - * - */ -typedef struct -{ - double R; ///< Coordinates of projection center - double D; ///< Coordinates of projection center - double Xs; ///< plate-scale in X direction - double Ys; ///< plate-scale in Y direction - psProjectionType type; ///< Projection type -} -psProjection; - -/** Mode for Offset calculation between two sky positions - * - * @see psSphereGetOffset, psSphereSetOffset - * - */ -typedef enum { - PS_SPHERICAL, ///< offset corresponds to an angular offset - PS_LINEAR ///< offset corresponds to a linear offset -} psSphereOffsetMode; - -/** The units of the offset - * - * @see psSphereGetOffset, psSphereSetOffset - * - */ -typedef enum { - PS_ARCSEC, ///< Arcseconds - PS_ARCMIN, ///< Arcminutes - PS_DEGREE, ///< Degrees - PS_RADIAN ///< Radians -} psSphereOffsetUnit; - -#include "psTime.h" - -/** Allocates a psPlane - * - * @return psPlane* resulting plane structure. - */ - -psPlane* psPlaneAlloc(void); - -/** Allocates a psSphere - * - * @return psSphere* resulting sphere structure. - */ -psSphere* psSphereAlloc(void); - -/** psSphereRot allocator which defines the rotation in terms of the coordinate - * of the pole and the rotation about that pole. - * - * @return psSphereRot* Newly allocated psSphereRot object - */ -psSphereRot* psSphereRotAlloc( - double alphaP, - double deltaP, - double phiP -); - -/** psSphereRot allocator which defines the rotation from the elements of the - * quaternion. - * - * @return psSphereRot* Newly allocated psSphereRot object - */ -psSphereRot* psSphereRotQuat( - double q0, - double q1, - double q2, - double q3 -); - -/** Allocates a psPlaneTransform transform. - * - * @return psPlaneTransform* resulting plane transform - */ -psPlaneTransform* psPlaneTransformAlloc( - int n1, ///< The order of the x term in the transform. - int n2 ///< The order of the y term in the transform. -); - -/** Applies the psPlaneTransform transform to a specified coordinate - * - * @return psPlane* resulting coordinate based on transform - */ -psPlane* psPlaneTransformApply( - psPlane* out, ///< a psPlane to recycle. If NULL, a new one is generated. - const psPlaneTransform* transform, ///< the transform to apply - const psPlane* coords ///< the coordinate to apply the transform above. -); - -/** Allocates a psPlaneDistort transform. - * - * @return psPlaneDistort* resulting plane distort transform - */ - -psPlaneDistort* psPlaneDistortAlloc( - int n1, ///< The order of the w term in the transform. - int n2, ///< The order of the x term in the transform. - int n3, ///< The order of the y term in the transform. - int n4 ///< The order of the z term in the transform. -); - - -/** Applies the psPlaneDistort transform to a specified coordinate - * - * @return psPlane* resulting coordinate based on transform - */ -psPlane* psPlaneDistortApply( - psPlane* out, ///< a psPlane to recycle. If NULL, a new one is generated. - const psPlaneDistort* distort, ///< the transform to apply - const psPlane* coords, ///< the coordinate to apply the transform above. - float mag, ///< third term -- maybe magnitude - float color ///< forth term -- maybe color -); - - -/** Applies the psSphereRot transform for a specified coordinate - * - * @return psSphere* resulting coordinate based on transform - */ -psSphere* psSphereRotApply( - psSphere* out, ///< a psSphere to recycle. If NULL, a new one is generated. - const psSphereRot* transform, ///< the transform to apply - const psSphere* coord ///< the coordinate to apply the transform above.x -); - -/** Combines two rotations to produce a single rotation which is equivalent of - * applying the first rotation and then the second. - * - * @return psSphereRot* new psSphereRot transform - */ -psSphereRot* psSphereRotCombine( - psSphereRot* out, - const psSphereRot* rot1, - const psSphereRot* rot2 -); - -/** Inverts a psSphereRot's rotation. - * - * @return psSphereRot* The inverted psSphereRot - */ -psSphereRot* psSphereRotInvert( - psSphereRot* rot ///< the psSphereRot to invert -); - -/** Creates the appropriate transform for converting from ICRS to Ecliptic - * coordinate systems. - * - * @return psSphereRot* transform for ICRS->Ecliptic coordinate systems - */ -psSphereRot* psSphereRotICRSToEcliptic( - const psTime* time ///< the time for which the resulting transform will be valid -); - -/** Creates the appropriate transform for converting from Ecliptic to ICRS - * coordinate systems. - * - * @return psSphereRot* transform for Ecliptic->ICRS coordinate systems - */ -psSphereRot* psSphereRotEclipticToICRS( - const psTime* time ///< the time for which the resulting transform will be valid -); - -/** Creates the appropriate transform for converting from ICRS to Galactic - * coordinate systems. - * - */ -psSphereRot* psSphereRotICRSToGalactic(void); - -/** Creates the appropriate transform for converting from Galactic to ICRS - * coordinate systems. - * - */ -psSphereRot* psSphereRotGalacticToICRS(void); - -/** Allocates memory for a psProjection structure - * - * @return psProjection* psProjection structure - */ -psProjection* psProjectionAlloc( - double R, ///< Right-ascension of projection center. - double D, ///< Declination of projection center. - double Xs, ///< Scale in x-dimension - double Ys, ///< Scale in y-dimension - psProjectionType type -); - -/** Projects a spherical coordinate to a linear coordinate system - * - * @return psPlane* projected coordinate - */ -psPlane* psProject( - const psSphere* coord, ///< coordinate to project - const psProjection* projection ///< parameters of the projection -); - -/** Reverse projection of a linear coordinate to a spherical coordinate system - * - * @return psPlane* projected coordinate - */ -psSphere* psDeproject( - const psPlane* coord, ///< coordinate to project - const psProjection* projection ///< parameters of the projection -); - -/** Determines the offset (RA,Dec) on the sky between two positions. - * - * Both an offset mode and an offset unit may be defined. The mode may be - * either PS_SPHERICAL, in which case the specified offset corresponds to an - * offset in angles, or it may be PS_LINEAR, in which case the offset - * corresponds to a linear offset in a local projection. The offset unit may - * be in one of PS_ARCSEC, PS_ARCMIN, PS_DEGREE, and PS_RADIAN, which - * specifies the units of the offset only. - * - * @return psSphere* the offset between position1 and position2 - */ -psSphere* psSphereGetOffset( - const psSphere* position1, ///< first position for calculating offset - const psSphere* position2, ///< second position for calculating offset - psSphereOffsetMode mode, ///< type of offset can be PS_SPHERICAL or PS_LINEAR - psSphereOffsetUnit unit ///< specifies the units of offset only -); - -/** Applies the given offset to a coordinate. - * - * Both an offset mode and an offset unit may be defined. The mode may be - * either PS_SPHERICAL, in which case the specified offset corresponds to an - * offset in angles, or it may be PS_LINEAR, in which case the offset - * corresponds to a linear offset in a local projection. The offset unit may - * be in one of PS_ARCSEC, PS_ARCMIN, PS_DEGREE, and PS_RADIAN, which - * specifies the units of the offset only. - * - * @return psSphere* the original position with the given offset applied. - */ -psSphere* psSphereSetOffset( - const psSphere* position, ///< coordinate of origin - const psSphere* offset, ///< coordinate of offset to apply - psSphereOffsetMode mode, ///< corresponds to an offset in angles or local projection - psSphereOffsetUnit unit ///< specifies the units of offset only -); - -/** Generates the complete spherical rotation to account for precession - * between two times. The equinoxes shall be Julian equinoxes. - * - * @return psSphere* the resulting spherical rotation - */ -psSphere* psSpherePrecess( - psSphere *coords, ///< coordinates (modified in-place) - const psTime *fromTime, ///< equinox of coords input - const psTime *toTime ///< equinox of coords output -); - -/** Takes a given transform and inverts it linearly if possible. - * - * @return psPlaneTransform - * the linearly inverted transform -*/ -psPlaneTransform *p_psPlaneTransformLinearInvert( - psPlaneTransform *transform ///< transform to invert -); - - -/** Takes a transform and tests whether or not it is a linear projection. - * - * @return psS32 - * the order of the projection -*/ -psS32 p_psIsProjectionLinear( - psPlaneTransform *transform ///< transform to test for linearity -); - -/** inverts a given transformation. - * - * It may assume that the input transformation is one-to-one, and that the - * inverse transformation may be specified through using polynomials of the - * same type and order as the forward transformation. In the event that the - * input transformation is linear, an exact solution may be calculated; - * otherwise nSamples samples in each axis, covering the region specified by - * region shall be used as a grid to fit the best inverse transformation. The - * function shall return NULL if it was unable to generate the inverse - * transformation; otherwise it shall return the inverse transformation. In - * the event that out is NULL, a new psPlaneTransform shall be allocated and - * returned. - * - * @return psPlaneTransform* the resulting inverted transform - */ -psPlaneTransform* psPlaneTransformInvert( - psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. - const psPlaneTransform *in, ///< transform to invert - psRegion region, ///< region to fit for non-linear transform inversion - int nSamples ///< number of samples in each axis for fit -); - -/** Creates a single transformation that has the effect of performing trans1 - * followed by trans2. - * - * psPlaneTransformCombine takes two transformations (trans1 and trans2) and - * returns a single transformation that has the effect of performing trans1 - * followed by trans2. In the event that the input transformation is linear, - * an exact solution may be calculated; otherwise nSamples samples in each - * axis, covering the region specified by region shall be used as a grid to - * fit the best inverse transformation. The function shall return NULL if it - * was unable to generate the transformation; otherwise it shall return the - * transformation. - * - * @return psPlaneTransform* resulting transformation - */ -psPlaneTransform* psPlaneTransformCombine( - psPlaneTransform *out, ///< a transform to recycle, or NULL if one is to be created. - const psPlaneTransform *trans1, ///< first transform to combine - const psPlaneTransform *trans2, ///< first transform to combine - psRegion region, ///< region to cover (for non-linear transforms) - int nSamples ///< number of samples on each axis (for non-linear transforms) -); - - -/** takes two arrays containing matched coordinates and returns the - * best-fitting transformation. - * - * psPlaneTransformFit takes two arrays containing matched coordinates (i.e., - * coordinates in the source array correspond to the coordinates in the dest - * array) and returns the best-fitting transformation. The source and dest - * will contain psCoords. In the event that the number of coordinates in each - * is not identical, the function shall generate a warning, and extra - * coordinates in the longer of the two shall be ignored. The trans transform - * may not be NULL, since it specifies the desired order, polynomial type and - * any polynomial terms to mask. nRejIter rejection iterations shall be - * performed, wherein coordinates lying more than sigmaClip standard - * deviations from the fit shall be rejected. - * - * @return bool TRUE if successful, otherwise FALSE. - */ -bool psPlaneTransformFit( - psPlaneTransform *trans, - const psArray *source, - const psArray *dest, - int nRejIter, - float sigmaClip -); -//XXX: need to add doxygen comments on the parameters above. -rdd - -/// @} - -#endif // #ifndef PS_COORD_H Index: trunk/psLib/src/astronomy/psPhotometry.h =================================================================== --- trunk/psLib/src/astronomy/psPhotometry.h (revision 4540) +++ (revision ) @@ -1,74 +1,0 @@ -/** @file psPhotometry.h -* -* @brief Contains basic photometric structures. -* -* This file defines the basic photometric structures. -* -* @ingroup Photometry -* -* @author GLG, MHPCC -* -* XXX: There is no code associated with this header file. Perhaps we should -* incorporate it into psAstrometry.h -* -* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-08 23:40:45 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ - -#ifndef PS_PHOTOMETRIC_H -#define PS_PHOTOMETRIC_H - -#include "psType.h" -#include "psFunctions.h" - -/// @addtogroup Photometry -/// @{ - -/** The photometric system description - * - * The photometric system is defined by the psPhotSystem structure. A - * photometric system is identified by a human-readable name (ie, SDSS.g, - * Landolt92.B, GPC1.OTA32.r). Each photometric system is given a unique - * identifier ID. Observations taken with a specific camera, detector, and - * filter represent their own photometric system, and it may be necessary to - * perform transformations between these systems. Photometric systems - * associated with observations from a specific camera/ detector/filter - * combination can be associated with those components. - * - */ - -typedef struct -{ - const psS32 ID; ///< ID number for this photometric system - const char *name; ///< Name of photometric system - const char *camera; ///< Camera for photometric system - const char *filter; ///< Filter used for photometric system - const char *detector; ///< Detector used for photometric system -} -psPhotSystem; - -/** Photometric system transformation - * - * This structure defines the transformation between two photometric systems. - * - */ - -typedef struct -{ - const psPhotSystem src; ///< Source photometric system - const psPhotSystem dst; ///< Destination photometric system - const psPhotSystem pP; ///< Primary color reference - const psPhotSystem pM; ///< Primary color reference - const psPhotSystem sP; ///< Secondary color reference - const psPhotSystem sM; ///< Secondary color reference - psF32 pA; ///< Color offset for references - psF32 sA; ///< Color offset for references - psPolynomial3D transform; ///< Transformation from source to destination -} -psPhotTransform; - -/// @} - -#endif // #ifndef PS_PHOTOMETRIC_H Index: trunk/psLib/src/astronomy/psTime.c =================================================================== --- trunk/psLib/src/astronomy/psTime.c (revision 4540) +++ (revision ) @@ -1,1598 +1,0 @@ -/** @file psTime.c - * - * @brief Definitions for time, time utilities, and conversion functions for use with psLib astronomy - * functions. - * - * A collection of functions are required by psLib to manipulate time data. These functions primarily consist - * of conversions between specific time formats. They use the UNIX timeval time system as the - * base upon which International Atomic Time (TAI) and Universal Time Coordinated (UTC) are calculated. - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.66 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-12 19:12:00 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#include -#include -#include -#include -#include - -#include "psTime.h" -#include "psError.h" -#include "psLogMsg.h" -#include "psMemory.h" -#include "psAbort.h" -#include "psImage.h" -#include "psCoord.h" -#include "psString.h" -#include "psMetadata.h" -#include "psMetadataIO.h" -#include "psLookupTable.h" -#include "psConstants.h" -#include "psAstronomyErrors.h" - -#include "config.h" - -#define MAX_STRING_LENGTH 256 - -/** Sidereal angular conversion from seconds to radians for GMST in seconds (i.e. pi/(180*240)) */ -#define S2R (7.272205216643039903848711535369e-5) - -/** Two times pi with double precision accuracy */ -#define TWOPI (2.0*M_PI) - -/** Conversion from radians to degrees */ -#define R2DEG = (180.0/M_PI) - - /** Maximum length of time string */ - #define MAX_TIME_STRING_LENGTH 256 - - /** Seconds per minute */ - #define SEC_PER_MINUTE 60.0 - - /** Seconds per hour */ - #define SEC_PER_HOUR (60.0*SEC_PER_MINUTE) - - /** Seconds per day */ - #define SEC_PER_DAY (24.0*SEC_PER_HOUR) - - /** Seconds per year */ - #define SEC_PER_YEAR (365.0*SEC_PER_DAY) - - /** Microseconds per day */ - #define NSEC_PER_DAY 86400000000000.0 - - /** Time metadata read from config file */ - static psMetadata *timeMetadata = NULL; - -// Offset to convert terrestrial time(TT) to international atomic time(TAI) -#define TAI_TT_OFFSET_SECONDS 32 -#define TAI_TT_OFFSET_NANOSECONDS 184000000 - -// Offset from converting to MJD -#define MJD_EPOCH_OFFSET 40587.0 - -// Offset for converting to JD -#define JD_EPOCH_OFFSET 2440587.5 - -// Offset of year 0000 from epoch -#define YEAR_0000_SEC -62125920000.0 - -// Offset of year 9999 from epoch -#define YEAR_9999_SEC 253202544000.0 - -/** Static function prototypes */ -static char *cleanString(char *inString, int sLen); -static char* getToken(char **inString, char *delimiter, psParseErrorType *status); -static psF64 searchTables(psF64 index, psU64 column, char *metadataTableNames[], - psU32 nTables, psLookupStatusType* status); -static psTime* convertTimeTAIUTC(psTime* time); -static psTime* convertTimeUTCTAI(psTime* time); -static psTime* convertTimeTAITT(psTime* time); -static psTime* convertTimeTTTAI(psTime* time); -static psTime* convertTimeUTCUT1(psTime* time); - -/** Removes leading and trailing whitespace and # characters from a string. The cleaned string is a new null - * terminated copy of the original input string. */ -static char *cleanString(char *inString, int sLen) -{ - char *ptrB = NULL; - char *ptrE = NULL; - char *cleaned = NULL; - - ptrB = inString; - - // Skip over leading # or whitespace - while (isspace(*ptrB) || *ptrB=='#') { - ptrB++; - } - - // Skip over trailing whitespace, null terminators, and # characters - ptrE = inString + sLen - 1; - while(isspace(*ptrE) || *ptrE=='\0' || *ptrE=='#') { - ptrE--; - } - - // Length, sLen, does not include '\0' - sLen = ptrE - ptrB + 1; - - // Adds '\0' to end of string and +1 to sLen - cleaned = psStringNCopy(ptrB, sLen); - - return cleaned; -} - -/** Returns cleaned token based on delimiter, but not including delimiter. Also changes the pointer location - * the beginning of the string. Tokens are newly allocated null terminated strings. */ -static char* getToken(char **inString, char *delimiter, psParseErrorType *status) -{ - char *cleanToken = NULL; - int sLen = 0; - - // Skip over leading whitespace - while(isspace(**inString)) { - (*inString)++; - } - - // Length of token, not including delimiter - sLen = strcspn(*inString, delimiter); - - if(sLen) { - - // Create new, cleaned, and null terminated token - cleanToken = cleanString(*inString, sLen); - - // Move to end of token - (*inString) += (sLen+1); - - } else if(**inString!='\0' && sLen==0) { - *status = PS_PARSE_ERROR_GENERAL; - } - - return cleanToken; -} - -// get the psTime.config filename by checking environment variable first, then original installation area. -char* p_psGetConfigFileName() -{ - char* filename = getenv("PS_CONFIG_FILE"); - - if (filename == NULL) { // environment variable not found - filename = PS_CONFIG_FILE_DEFAULT; // this should come from configure.ac - } - - return filename; -} - - -// Searches time tables in priority order and performs interpolation if input index value is within a table. -// If the index value is out of range, the status is set accordingly. -static psF64 searchTables(psF64 index, psU64 column, char *metadataTableNames[], - psU32 nTables, psLookupStatusType* status) -{ - char* tableName = NULL; - psF64 result = NAN; - psLookupTable* table = NULL; - psMetadataItem* tableMetadataItem = NULL; - - // Check time metadata. Function call reports errors. - if(timeMetadata == NULL) { - if(!p_psTimeInit(p_psGetConfigFileName())) - *status = PS_LOOKUP_ERROR; - return NAN; - } - - // Search each table in priority order: daily, eopc,finals - for(psS32 i = 0; i < nTables; i++) { - - // Get table name from list of tables to search - tableName = metadataTableNames[i]; - - // Lookup table name in time metadata - tableMetadataItem = psMetadataLookup(timeMetadata, tableName); - - // Check if table not a metadata item - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - tableName); - *status = PS_LOOKUP_ERROR; - return NAN; - } - - // Get table from metadata - table = (psLookupTable*)tableMetadataItem->data.V; - - // Check that table is not NULL - PS_ASSERT_PTR_NON_NULL(table,NAN); - - // Check if index within to/from range - if(index >= table->validFrom ) { - if(index <= table->validTo) { - // Attempt to interpolate table - result = psLookupTableInterpolate(table, index, column); - *status = PS_LOOKUP_SUCCESS; - if(!isnan(result)) { - break; - } - } else { - *status = PS_LOOKUP_PAST_BOTTOM; - } - } else { - *status = PS_LOOKUP_PAST_TOP; - } - } - - return result; -} - -bool p_psTimeInit(const char *fileName) -{ - psS32 numLines = 0; - bool foundTable = false; - char *tableDir = NULL; - char *tableNames = NULL; - char *tableFormats = NULL; - char *namesPtr = NULL; - char *formatPtr = NULL; - char *metadataNamesPtr = NULL; - char *tableName = NULL; - char *tableFormat = NULL; - char *fullTableName = NULL; - psS32 i = 0; - psS32 j = 0; - psS32 numTables = 0; - psU32 nFail = 0; - psVector *tablesFrom = NULL; - psVector *tablesTo = NULL; - psVector *tablesIndex = NULL; - psMetadataItem *metadataItem = NULL; - psLookupTable *table = NULL; - psParseErrorType status = PS_PARSE_SUCCESS; - char metadataTableNames[4][MAX_STRING_LENGTH] = {"daily", "eopc", "finals", "tai"}; - - // Read time config file - timeMetadata = psMetadataConfigParse(timeMetadata, &nFail, fileName, true); - if(timeMetadata == NULL) { - return false; - } else if(nFail != 0) { - return false; - } - - // Get number of tables - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.n"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.n"); - return false; - } - numTables = (psS32)metadataItem->data.S32; - - // Get lower range of tables - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.from"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.from"); - return false; - } - tablesFrom = psVectorCopy(tablesFrom, metadataItem->data.V, PS_TYPE_F64); - if(tablesFrom->n != numTables) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesFrom->n, numTables); - psFree(tablesFrom); - return false; - } - - // Get upper range of tables - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.to"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.to"); - psFree(tablesFrom); - return false; - } - tablesTo = psVectorCopy(tablesTo, metadataItem->data.V, PS_TYPE_F64); - if(tablesTo->n != numTables) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_VECTOR, tablesTo->n, numTables); - psFree(tablesFrom); - psFree(tablesTo); - return false; - } - - // Get index columns for the tables - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.index"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.index"); - psFree(tablesFrom); - psFree(tablesTo); - return false; - } - tablesIndex = psVectorCopy(tablesIndex, metadataItem->data.V, PS_TYPE_S32); - if(tablesIndex->n != numTables) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_BAD_VECTOR,tablesIndex->n,numTables); - psFree(tablesFrom); - psFree(tablesTo); - psFree(tablesIndex); - return false; - } - - // Get path to time data files - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.dir"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.dir"); - psFree(tablesFrom); - psFree(tablesTo); - psFree(tablesIndex); - return false; - } - tableDir = psStringCopy(metadataItem->data.V); - - // Table file names - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.files"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.files"); - - psFree(tablesFrom); - psFree(tablesTo); - psFree(tablesIndex); - psFree(tableDir); - return false; - } - tableNames = psStringCopy(metadataItem->data.V); - - // Get table format strings - metadataItem = psMetadataLookup(timeMetadata, "psLib.time.tables.format"); - if(metadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.format"); - psFree(tablesFrom); - psFree(tablesTo); - psFree(tablesIndex); - psFree(tableDir); - psFree(tableNames); - return false; - } - tableFormats = psStringCopy(metadataItem->data.V); - formatPtr = tableFormats; - - // Read time tables - namesPtr = tableNames; - while((tableName=getToken(&namesPtr, " ", &status)) != NULL) { - - // Form path with table name, adding one to length for last '/' that may not occur - // in string in cong file - fullTableName = (char*)psAlloc(strlen(tableDir)+strlen(tableName)+1+1); - - // Old strings may come back from psAlloc(), so set initial position to EOL - fullTableName[0]='\0'; - strcat(fullTableName, tableDir); - strcat(fullTableName, "/"); - strcat(fullTableName, tableName); - - // Get table format - tableFormat = getToken(&formatPtr,",",&status); - if(tableFormat == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.tables.format"); - } - - // Create and read table - if(i < numTables) { - table = psLookupTableAlloc(fullTableName, (const char*)tableFormat, tablesIndex->data.S32[i]); - numLines = psLookupTableRead(table); - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_BAD_TABLE_COUNT, i+1, numTables); - } - - // Place tables into metadata slightly altered names as keys to create consistent naming conventions - foundTable = false; - for(j=0; jsec = 0; - outTime->nsec = 0; - outTime->type = type; - outTime->leapsecond = false; - - return outTime; -} - -psTime* psTimeGetNow(psTimeType type) -{ - struct timeval now; - psTime *time = NULL; - - // Allocate psTime struct - time = psTimeAlloc(type); - - // Verify time structure allocated - if(time == NULL) { - return NULL; - } - - // Get the system time - if (gettimeofday(&now, (struct timezone *)0) == -1) { - psError(PS_ERR_OS_CALL_FAILED, true, - PS_ERRORTEXT_psTime_GET_TOD_FAILED); - return NULL; - } - - // Convert timeval time to psTime - time->sec = now.tv_sec; - time->nsec = now.tv_usec*1000; - - // Add most leapseconds to UTC time to get TAI time if necessary - if(type == PS_TIME_TAI) { - time->sec += p_psTimeGetTAIDelta(time); - } - - return time; -} - -static psTime* convertTimeTAIUTC(psTime* time) -{ - psF64 deltaTAI = 0.0; - psS64 deltaSec = 0; - psU32 deltaNsec = 0; - psF64 deltaUTC = 0.0; - - // Determine delta to convert between UTC and TAI - deltaTAI = p_psTimeGetTAIDelta(time); - deltaSec = (psS64)(deltaTAI); - deltaNsec = (psU32)((deltaTAI - (psF64)deltaSec) * 1e9); - - // Determine seconds - time->sec -= deltaSec; - - // Check for underflow in nsec - if(deltaNsec > time->nsec) { - // Borrow second - time->nsec += 1e9; - time->sec--; - } - - // Determine nsec - time->nsec -= deltaNsec; - - // Check for overflow in nsec - if(time->nsec >= 1e9) { - time->nsec -= 1e9; - time->sec++; - } - - // Set new type - time->type = PS_TIME_UTC; - - // Check if leapsecond present in delta - deltaUTC = p_psTimeGetTAIDelta(time); - if(fabs(deltaTAI-deltaUTC) >= 1.0) { - time->sec++; - } - - return time; -} - -static psTime* convertTimeUTCTAI(psTime* time) -{ - psF64 delta = 0.0; - psS64 deltaSec = 0; - psU32 deltaNsec = 0; - - // Determine delta to convert between UTC and TAI - delta = p_psTimeGetTAIDelta(time); - - deltaSec = (psS64)(delta); - deltaNsec = (psU32)((delta - (psF64)deltaSec) * 1e9); - - // Determine seconds - time->sec += deltaSec; - - // Determine nsec - time->nsec += deltaNsec; - - // Check for overflow in nsec - if(time->nsec >= 1e9) { - time->nsec -= 1e9; - time->sec++; - } - - // Set new type - time->type = PS_TIME_TAI; - - return time; -} - -static psTime* convertTimeTAITT(psTime* time) -{ - // Add TT offset - time->sec += TAI_TT_OFFSET_SECONDS; - time->nsec += TAI_TT_OFFSET_NANOSECONDS; - - // Check for overflow in nsec - if(time->nsec >= 1e9) { - time->nsec -= 1e9; - time->sec++; - } - - // Set new type - time->type = PS_TIME_TT; - - return time; -} - -static psTime* convertTimeTTTAI(psTime* time) -{ - // Subtract TT offset - time->sec -= TAI_TT_OFFSET_SECONDS; - - // Check for nsec underflow - if(TAI_TT_OFFSET_NANOSECONDS > time->nsec) { - // Borrow second - time->sec--; - time->nsec += 1e9; - } - time->nsec -= TAI_TT_OFFSET_NANOSECONDS; - - // Check for overflow in nsec - if(time->nsec >= 1e9) { - time->nsec -= 1e9; - time->sec++; - } - - // Set new type - time->type = PS_TIME_TAI; - - return time; -} - -static psTime* convertTimeUTCUT1(psTime* time) -{ - psS64 ut1utc = 0; - - // Get UT1-UTC value - ut1utc = (psS64)(psTimeGetUT1Delta(time,PS_IERS_A) * 1e9); - - // Since UTC is within 0.9 sec of UT1 then nsec member is the member affected - if((ut1utc < 0) && (abs(ut1utc) > time->nsec)) { - // Borrow from sec - time->sec--; - if(time->leapsecond) { - time->leapsecond = false; - } else { - time->leapsecond = psTimeIsLeapSecond(time); - } - // Add to nsec - time->nsec += 1e9; - } - time->nsec += ut1utc; - - // Check for overflow in nsec - if(time->nsec >= 1e9) { - time->nsec -= 1e9; - time->sec++; - if(time->leapsecond) { - time->leapsecond = false; - time->sec--; - } else { - time->leapsecond = psTimeIsLeapSecond(time); - } - } - - // Set new type - time->type = PS_TIME_UT1; - - return time; -} - -psTime* psTimeConvert(psTime *time, psTimeType type) -{ - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),time); - - // If the input type is UT1 then return time and generate error message - if(time->type == PS_TIME_UT1) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Cannot convert from UT1 time type"); - return time; - } - - // If the time to convert to is the same as psTime the return time - if (time->type == type) { - return time; - } - - // Convert from TAI to UTC, TT, UT1 - if(time->type == PS_TIME_TAI) { - // Convert from TAI to UTC - if(type == PS_TIME_UTC) { - time = convertTimeTAIUTC(time); - // Convert from TAI to TT - } else if(type == PS_TIME_TT) { - time = convertTimeTAITT(time); - // Convert from TAI to UT1 - } else if(type == PS_TIME_UT1) { - // Convert to UTC first - time = convertTimeTAIUTC(time); - // Convert UTC to UT1 - time = convertTimeUTCUT1(time); - // Convert from TAI to unknown time type - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); - } - // Convert from TT to TAI, UTC, UT1 - } else if(time->type == PS_TIME_TT) { - // Convert from TT to UTC - if(type == PS_TIME_UTC) { - // Convert to TAI time first - time = convertTimeTTTAI(time); - // Convert from TAI to UTC - time = convertTimeTAIUTC(time); - // Convert from TT to TAI - } else if(type == PS_TIME_TAI) { - time = convertTimeTTTAI(time); - // Convert from TT to UT1 - } else if(type == PS_TIME_UT1) { - // Convert to UTC first - // Convert to TAI time first - time = convertTimeTTTAI(time); - // Convert from TAI to UTC - time = convertTimeTAIUTC(time); - // Convert from UTC to UT1 - time = convertTimeUTCUT1(time); - // Convert from TT to unknown time type - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); - } - // Convert from UTC to TAI, TT, UT1 - } else if(time->type == PS_TIME_UTC) { - // Convert UTC to TAI - if(type == PS_TIME_TAI) { - time = convertTimeUTCTAI(time); - // Convert UTC to TT - } else if(type == PS_TIME_TT) { - // Convert to TAI time first - time = convertTimeUTCTAI(time); - // Convert TAI to TT - time = convertTimeTAITT(time); - // Convert UTC to UT1 - } else if(type == PS_TIME_UT1) { - time = convertTimeUTCUT1(time); - // Convert UTC to unknown time type - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, type); - } - // Convert unknown time type - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_UNKNOWN, time->type); - } - - return time; -} - -double psTimeToLMST(psTime *time, double longitude) -{ - psF64 jdTdtDays = 0.0; - psF64 jdUt1Days = 0.0; - psF64 mjdUt1Days = 0.0; - psF64 lmstRad = 0.0; - psF64 fracDays = 0.0; - psF64 gmstRad = 0.0; - psF64 t = 0.0; - psF64 tu = 0.0; - psF64 const1 = 24110.5493771; - psF64 const2 = 8639877.3173760; - psF64 const3 = 307.4771600; - psF64 const4 = 0.0931118; - psF64 const5 = -0.0000062; - psF64 const6 = 0.0000013; - psTime *tdtTime = NULL; - psTime *ut1Time = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); - - // Verify input time is not in UT1 seconds - if(time->type == PS_TIME_UT1) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time->type); - return NAN; - } - - // Determine time reference to TT - tdtTime = psTimeAlloc(time->type); - tdtTime->sec = time->sec; - tdtTime->nsec = time->nsec; - tdtTime->leapsecond = time->leapsecond; - tdtTime = psTimeConvert(tdtTime,PS_TIME_TT); - - // Determine time reference to UT1 - ut1Time = psTimeAlloc(time->type); - ut1Time->sec = time->sec; - ut1Time->nsec = time->nsec; - ut1Time->leapsecond = time->leapsecond; - ut1Time = psTimeConvert(ut1Time,PS_TIME_UT1); - - // Calculate UT1 as Julian Centuries since J2000.0 - jdUt1Days = psTimeToJD(ut1Time); - mjdUt1Days = psTimeToMJD(ut1Time); - t = (jdUt1Days - 2451545.0)/36525.0; - - // Calculate TDT as Julian centuries since J2000.0 - jdTdtDays = psTimeToJD(tdtTime); - tu = (jdTdtDays - 2451545.0)/36525.0; - - // Calculate fractional part of MJD - fracDays = fmod(mjdUt1Days, 1.0); - - // Calculate Greenwich Mean Sidereal Time (GMST) in radians. - // Equation set up to minimize multiplications. - gmstRad = fracDays*TWOPI - + (const1+const2*tu+t*(const3+t*(const4+t*(const5+const6*t))))*S2R; - - // Place GMST between 0 and 2*pi - gmstRad = fmod(gmstRad, TWOPI); - - // Calculate Local Mean Sidereal Time (LMST) in radians - lmstRad = gmstRad + longitude; - - // Free temporary structs - psFree(ut1Time); - psFree(tdtTime); - - return lmstRad; -} - -double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin) -{ - psU32 nTables = 2; - psF64 mjd = 0.0; - psF64 result = 0.0; - psU64 tableColumn = 0; - psF64 dut2ut1 = 0.0; - psF64 t = 0.0; - psVector* dut = NULL; - psMetadataItem* tableMetadataItem = NULL; - psLookupStatusType status = PS_LOOKUP_SUCCESS; - char* metadataTableNames[2] = {"dailyTable", "finalsTable"}; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); - - // Check for invalid bulletin specified - if((bulletin != PS_IERS_A) && (bulletin != PS_IERS_B)) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,"Invalid bulletin specified %d",bulletin); - return NAN; - } - - // Set lookup table column based on Bullentin - if(bulletin == PS_IERS_A) { - tableColumn = 3; - } else { - tableColumn = 6; - } - - // Attempt to find value through table lookup and interpolation - mjd = psTimeToMJD(time); - result = searchTables(mjd,tableColumn,metadataTableNames,nTables,&status); - - // Value could not be found through table lookup and interpolation - if(status == PS_LOOKUP_PAST_TOP) { - - // Date too early for tables. Get default time delta value from metadata, and issue warning. - psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"UT1-UTC"); - - // Lookup value from time metadata loaded from psTime.config - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.dut"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.before.dut"); - return NAN; - } - result = tableMetadataItem->data.F64; - - } else if(status == PS_LOOKUP_PAST_BOTTOM) { - /* Date too late for tables. Issue warning and use following formulae for predicting - ahead of the most recent available table entry. - ut1-utc = [0] + [1]*(MJD - [2]) - (ut2-ut1) - [0, 1, 2] = @psLib.time.predict.dut - ut2-ut1 = 0.022 sin(2*pi*t) - 0.012 cos(2*pi*t) - 0.006 sin(4*pi*t) + 0.007 cos(4*pi*t) - t = 2000.0 + (MJD - 51544.03)/365.2422 - */ - // Generate warning of postdate information - psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "UT1-UTC"); - - // Lookup values to calculate prediction - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.dut"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.predict.dut"); - return NAN; - } - dut = (psVector*)tableMetadataItem->data.V; - PS_ASSERT_PTR_NON_NULL(dut,NAN); - - // Calculate predication of future UT1-UTC - t = 2000.0 + (mjd - 51544.03)/365.2422; - dut2ut1 = 0.022*sin(TWOPI*t) - 0.012*cos(TWOPI*t) - 0.006*sin(4.0*M_PI*t) + 0.007*cos(4.0*M_PI*t); - result = dut->data.F64[0] + dut->data.F64[1]*(mjd - dut->data.F64[2]) - dut2ut1; - - } else if(status != PS_LOOKUP_SUCCESS) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); - return NAN; - } - - return result; -} - -psSphere* p_psTimeGetPoleCoords(const psTime* time) -{ - psU32 nTables = 3; - psF64 x = 0.0; - psF64 y = 0.0; - psF64 mjd = 0.0; - psF64 a = 0.0; - psF64 c = 0.0; - psF64 mjdPred = 0.0; - psSphere* output = NULL; - psLookupStatusType xStatus = PS_LOOKUP_SUCCESS; - psLookupStatusType yStatus = PS_LOOKUP_SUCCESS; - psMetadataItem *tableMetadataItem = NULL; - char *metadataTableNames[3] = {"dailyTable", "eopcTable", "finalsTable"}; - psVector *xp = NULL; - psVector *yp = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - - if(time->type != PS_TIME_TAI) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type); - return NULL; - } - - // Attempt to find value through table lookup and interpolation - mjd = psTimeToMJD(time); - // x = searchTables(mjd, 0, &xStatus, metadataTableNames, nTables); - x = searchTables(mjd, 0, metadataTableNames, nTables,&xStatus); - // y = searchTables(mjd, 0, &yStatus, metadataTableNames, nTables); - y = searchTables(mjd, 0, metadataTableNames, nTables,&yStatus); - - // Value could not be found through table lookup and interpolation - if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP) { - - // Date too earlier for tables. Get default polar coodinate values from metadata, and issue warning. - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES, mjd, "polar motion"); - - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.xp"); - return NULL; - } - x = tableMetadataItem->data.F64; - - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.yp"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.before.yp"); - return NULL; - } - y = tableMetadataItem->data.F64; - - } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM) { - - /* Date too late for tables. Issue warning and use following formulae for predicting - ahead of the most recent available table entry. - x = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c - [0], [1], [2], [3] = @psLib.time.predict.xp - y = [0] + [1]*cos a + [2]*sin a + [3]*cos c + [4]*sin c - [0], [1], [2], [3] = @psLib.time.predict.yp - a = 2*pi*(mjd - pslib.time.predict.mjd)/365.25 - c = 2*pi*(mjd - pslib.time.predict.mjd)/435.0 - */ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TIME_POSTDATES_TABLES, mjd, "polar motion"); - - // Get predicted MJD - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, - "psLib.time.predict.mjd"); - return NULL; - } - mjdPred = tableMetadataItem->data.F64; - - // Get xp - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp"); - return NULL; - } - xp = (psVector*)tableMetadataItem->data.V; - PS_ASSERT_PTR_NON_NULL(xp,NULL); - - // Get yp - tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp"); - return NULL; - } - yp = (psVector*)tableMetadataItem->data.V; - PS_ASSERT_PTR_NON_NULL(yp,NULL); - - // Calculate "a" and "c" constants - a = TWOPI*(mjd - mjdPred)/365.25; - c = TWOPI*(mjd - mjdPred)/435.0; - - // Calculate x and y polar coordinates - x = xp->data.F64[0] + - xp->data.F64[1]*cos(a) + - xp->data.F64[2]*sin(a) + - xp->data.F64[3]*cos(c) + - xp->data.F64[4]*sin(c); - - y = yp->data.F64[0] + - yp->data.F64[1]*cos(a) + - yp->data.F64[2]*sin(a) + - yp->data.F64[3]*cos(c) + - yp->data.F64[4]*sin(c); - - } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); - return NULL; - } - - // Create output sphere and convert arcsec to radians (i.e. x/60/60*PS_PI/180) - output = psAlloc(sizeof(psSphere)); - output->r = x * M_PI / 648000.0; - output->d = y * M_PI / 648000.0; - - return output; -} - -psF64 p_psTimeGetTAIDelta(const psTime *time) -{ - psF64 jd = 0.0; - psF64 mjd = 0.0; - psF64 out = 0.0; - psF64 const1 = 0.0; - psF64 const2 = 0.0; - psF64 const3 = 0.0; - psLookupTable* table = NULL; - psMetadataItem *tableMetadataItem = NULL; - psVector *results = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); - - // Check time metadata - if(timeMetadata == NULL) { - if(!p_psTimeInit(p_psGetConfigFileName())) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_FILE_NOT_FOUND, "psTime.config"); - return 0.0; - } - } - - // Get table from metadata - tableMetadataItem = psMetadataLookup(timeMetadata, "taiTable"); - if(tableMetadataItem == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "taiTable"); - return 0.0; - } - table = (psLookupTable*)tableMetadataItem->data.V; - PS_ASSERT_PTR_NON_NULL(table,0); - - // Determine Julian and modified Julian dates used in table lookup and time delta calculation - jd = psTimeToJD(time); - mjd = psTimeToMJD(time); - - // Set ceiling of the julian date to the last entry in the lookup table - if(table->validTo < jd) { - jd = table->validTo; - } - - // Interpolation of look up table - results = psLookupTableInterpolateAll(table, jd); - - // Check for successful interpolation - if(results == NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); - } - - // Set constants from table - const1 = results->data.F64[1]; - const2 = results->data.F64[2]; - const3 = results->data.F64[3]; - - // If const3 not equal to zero solve for difference else floor of const1 - if(fabs(const3-0.0) > FLT_EPSILON) { - out = const1 + (mjd - const2) * const3; - } else { - out = floor(const1); - } - - psFree(results); - - return out; -} - -long psTimeLeapSecondDelta(const psTime *time1, const psTime *time2) -{ - psS64 diff = 0; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time1,0); - PS_ASSERT_PTR_NON_NULL(time2,0); - PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0); - PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0); - diff = abs((psS64)p_psTimeGetTAIDelta((psTime*)time1)-(psS64)p_psTimeGetTAIDelta((psTime*)time2)); - - return diff; -} - -bool psTimeIsLeapSecond(const psTime* utc) -{ - psTime* prevUtc = NULL; - psBool returnValue = false; - - // Check for valid time - PS_ASSERT_PTR_NON_NULL(utc,false); - - // Verify time is UTC type - if(utc->type != PS_TIME_UTC) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,utc->type); - return false; - } - - // Allocate time to hold utc - 1 second - prevUtc = psTimeAlloc(PS_TIME_UTC); - prevUtc->sec = utc->sec - 1; - - // Check the absolute difference between the two times for leapsecond - if(psTimeLeapSecondDelta(utc,prevUtc) >= 1.0) { - returnValue = true; - } else { - returnValue = false; - } - - // Free prevUtc - psFree(prevUtc); - - return returnValue; -} - -double psTimeToJD(const psTime *time) -{ - psF64 jd = NAN; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); - - // Julian date conversion - if(time->sec < 0) { - // psTime earlier than epoch - jd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET; - } else { - // psTime greater than epoch - jd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + JD_EPOCH_OFFSET; - } - - return jd; -} - -double psTimeToMJD(const psTime *time) -{ - psF64 mjd = NAN; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NAN); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NAN); - - // Modified Julian date conversion - if(time->sec < 0) { - // psTime earlier than epoch - mjd = time->sec / SEC_PER_DAY - time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET; - } else { - // psTime greater than epoch - mjd = time->sec / SEC_PER_DAY + time->nsec / NSEC_PER_DAY + MJD_EPOCH_OFFSET; - } - - return mjd; -} - -psString psTimeToISO(const psTime *time) -{ - psS32 ds = 0; - char *timeString = NULL; - char *tempString = NULL; - struct tm *tmTime = NULL; - time_t sec; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - - // Check valid year range - PS_ASSERT_LONG_WITHIN_RANGE(time->sec,YEAR_0000_SEC,YEAR_9999_SEC,NULL) - - // Allocate temp strings - tempString = psAlloc(MAX_TIME_STRING_LENGTH); - timeString = psAlloc(MAX_TIME_STRING_LENGTH); - - // Convert nanoseconds to decaseconds - ds = time->nsec / 100000000; - sec = time->sec; - - // If leapsecond use previous day - if(time->leapsecond) { - sec--; - } - - // tmTime variable is statically allocated, no need to free - tmTime = gmtime(&sec); - - // Converts psTime to YYYY-MM-DDThh:mm:ss.sss in string form - if (!strftime(tempString, MAX_TIME_STRING_LENGTH, "%Y-%m-%dT%H:%M:%S", tmTime)) { - psError(PS_ERR_OS_CALL_FAILED, true, PS_ERRORTEXT_psTime_CONVERT_TIME_TO_STRING_FAILED); - } - - // Check if time is UTC and leapsecond - if(((time->type==PS_TIME_UTC)||(time->type==PS_TIME_UT1)) && (time->leapsecond)) { - // Modify second to be 60 - tempString[17] = '6'; - tempString[18] = '0'; - } - - // Create string with milliseconds - if (snprintf(timeString, MAX_TIME_STRING_LENGTH, "%s,%1dZ", tempString, ds) < 0) { - psError(PS_ERR_OS_CALL_FAILED, true, PS_ERRORTEXT_psTime_APPEND_MSEC_FAILED); - } - psFree(tempString); - - return timeString; -} - -struct timeval* psTimeToTimeval(const psTime *time) -{ - struct timeval *timevalTime = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->sec,0,INT32_MAX,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - - // Allocate structure timeval - timevalTime = (struct timeval*)psAlloc(sizeof(struct timeval)); - - // Set structure members - timevalTime->tv_sec = time->sec; - timevalTime->tv_usec = time->nsec / 1000; - - return timevalTime; -} - -/* -struct tm* p_psTimeToTM(const psTime *time) -{ - psS64 cent = 0; - psS64 year = 0; - psS64 month = 0; - psS64 day = 0; - psS64 hour = 0; - psS64 minute = 0; - psS64 seconds = 0; - psS64 temp = 0; - struct tm* tmTime = NULL; - - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - - seconds = time->sec%60; - minute = time->sec/60%60; - hour = time->sec/3600%24; - day = (time->sec+62135596800)/86400; - - // Add 306 days to make relative to Mar 1, 0; also adjust day to be within a range (1..2**28-1) where our - // calculations will work with 32bit ints - if(day > (pow(2, 28)-307)) - { - temp = (day - 146097+306)/146097+1; // Avoid overflow if day close to maxint - day -= temp * 146097-306; - } else if((day += 306) <= 0) - { - temp = -( -day / 146097 + 1); // Avoid ambiguity in C division of negatives - day -= temp * 146097; - } - - cent = (day*4-1)/146097; // Calc number of centuries day is after 29 Feb of yr 0 - day -= cent*146097/4; // 4 centuries = 146097 days - year = (day*4-1)/1461; // Calc number of years into the century - day -= year*1461/4; // Again March-based (4 yrs =\u02dc 146[01] days) - month = (day*12+1093)/367; // Get the month (3..14 represent March through - day -= (month*367-1094)/12; // February of following year) - year += cent*100+temp*400; // Get the real year, which is off by - - // One if month is January or February - if(month > 12) - { - year++; - month -= 12; - } - - // Allocate output - tmTime = (struct tm*)psAlloc(sizeof(struct tm)); - - tmTime->tm_year = year - 1900; - tmTime->tm_mon = month - 1; - tmTime->tm_mday = day + 1; - tmTime->tm_hour = hour; - tmTime->tm_min = minute; - tmTime->tm_sec = seconds; - tmTime->tm_isdst = -1; - - return tmTime; -} -*/ - -psTime* psTimeFromJD(double jd) -{ - psF64 days = 0.0; - psF64 seconds = 0.0; - psTime *outTime = NULL; - - // Allocate psTime struct - outTime = psTimeAlloc(PS_TIME_TAI); - - // Julian date conversion courtesy of Eugene Magnier - days = jd - 2440587.5; - seconds = days * SEC_PER_DAY; - if(seconds < 0.0) { - outTime->nsec = (seconds - (psS64)seconds) * -1000000000.0; // psTime earlier than epoch - } else { - outTime->nsec = (seconds - (psS64)seconds) * 1000000000.0; // psTime greater than epoch - } - outTime->sec = seconds; - - // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); - - return outTime; -} - -psTime* psTimeFromMJD(double mjd) -{ - psF64 days = 0.0; - psF64 seconds = 0.0; - psTime *outTime = NULL; - - // Allocate psTime struct - outTime = psTimeAlloc(PS_TIME_TAI); - - // Modified Julian date conversion courtesy of Eugene Magnier - days = mjd - 40587.0; - seconds = days * SEC_PER_DAY; - - if(seconds < 0.0) { - outTime->nsec = (seconds - (psS64)seconds) * -1000000000.0; // psTime earlier than epoch - } else { - outTime->nsec = (seconds - (psS64)seconds) * 1000000000.0; // psTime greater than epoch - } - outTime->sec = seconds; - - // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),NULL); - - return outTime; -} - -psTime* psTimeFromISO(const char *input) -{ - psS32 millisecond; - struct tm tmTime; - psTime *outTime = NULL; - - // Check for NULL string - PS_ASSERT_PTR_NON_NULL(input,NULL); - - // Convert YYYY-MM-DDThh:mm:ss.sss in string form to tm time - if (sscanf(input, "%d-%d-%dT%d:%d:%d,%d", &tmTime.tm_year, &tmTime.tm_mon, &tmTime.tm_mday, - &tmTime.tm_hour, &tmTime.tm_min, &tmTime.tm_sec,&millisecond) < 7) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_ISOTIME_MALFORMED, input); - return NULL; - } - - PS_ASSERT_INT_NONNEGATIVE(tmTime.tm_year, outTime); - PS_ASSERT_INT_WITHIN_RANGE(tmTime.tm_mon,1,12,outTime); - PS_ASSERT_INT_WITHIN_RANGE(tmTime.tm_mday,1,31,outTime); - PS_ASSERT_INT_WITHIN_RANGE(tmTime.tm_hour,0,23,outTime); - PS_ASSERT_INT_WITHIN_RANGE(tmTime.tm_min,0,59,outTime); - PS_ASSERT_INT_WITHIN_RANGE(tmTime.tm_sec,0,59,outTime); - PS_ASSERT_INT_WITHIN_RANGE(millisecond,0,999,outTime); - - tmTime.tm_year -= 1900; - tmTime.tm_mon--; - tmTime.tm_isdst = -1; - - // Convert tm time to psTime - outTime = p_psTimeFromTM(&tmTime); - outTime->nsec = millisecond * 1000000; - - return outTime; -} - -psTime* psTimeFromTT(psS64 sec, psU32 nsec) -{ - psTime* outTime = NULL; - - // Verify nsec within range - PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL); - - // Allocate psTime data - outTime = psTimeAlloc(PS_TIME_TT); - - // Set data members - outTime->sec = sec; - outTime->nsec = nsec; - - // Return data structure - return outTime; -} - -psTime* psTimeFromUTC(psS64 sec, psU32 nsec, bool leapsecond) -{ - psTime* outTime = NULL; - - // Verify nsec within range - PS_ASSERT_INT_WITHIN_RANGE(nsec,0,(psU32)((1e9)-1),NULL); - - // Allocate psTime data - outTime = psTimeAlloc(PS_TIME_UTC); - - // Set data members - outTime->sec = sec; - outTime->nsec = nsec; - - // Set leapsecond flag if necessary - outTime->leapsecond = psTimeIsLeapSecond(outTime); - - return outTime; -} - -psTime* psTimeFromTimeval(const struct timeval *input) -{ - psTime *outTime = NULL; - - - // Error check - PS_ASSERT_PTR_NON_NULL(input,NULL); - - // Allocate psTime struct - outTime = psTimeAlloc(PS_TIME_TAI); - - // Convert to psTime - outTime->sec = input->tv_sec; - outTime->nsec = input->tv_usec * 1000; - - // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); - - return outTime; -} - -psTime* p_psTimeFromTM(const struct tm* time) -{ - psS64 year; - psS64 month; - psS64 day; - psS64 hour; - psS64 minute; - psS64 seconds; - psS64 temp; - psTime *outTime = NULL; - - // Error check - PS_ASSERT_PTR_NON_NULL(time,NULL); - - // Allocate psTime struct - outTime = psTimeAlloc(PS_TIME_TAI); - - // Extract data from TM struct - year = time->tm_year + 1900; - month = time->tm_mon + 1; - day = time->tm_mday; - hour = time->tm_hour; - minute = time->tm_min; - seconds = time->tm_sec; - - // Make month in range 3..14 (treat Jan & Feb as months 13..14 of prev year) - if( month <= 2 ) - { - temp = (14 - month) / 12; - // year -= (temp = (14 - month) / 12); - year -= temp; - month += 12 * temp; - } else if(month > 14) - { - temp = (month - 3) / 12; - // year += (temp = (month - 3) / 12); - year += temp; - month -= 12 * temp; - } - - // Make year positive - if (year < 0 ) - { - day -= 146097 * (temp = (399 - year) / 400); - year += 400 * temp; - } - - // Add day of month, days of previous 0-11 month period that began w/March, days of previous 0-399 year - // period that began w/March of a 400-multiple year), days of any 400-year periods before that, and 306 - // days to adjust from Mar 1, year 0-relative to Jan 1, year 1-relative. Add hours, minutes, and seconds. - day += (month * 367 - 1094) / 12 + year % 100 * 1461 / 4 + (year/100 * 36524 + year/400) - 306; - outTime->sec = (((day - 1) * SEC_PER_DAY) - 62135596800) + hour*SEC_PER_HOUR - + minute*SEC_PER_MINUTE + seconds; - - // C's TM does not define a microsecond field. Microseconds must be manipulated by calling function. - outTime->nsec = 0; - - // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); - - return outTime; -} - -psTime* psTimeMath(const psTime *time, double delta) -{ - psF64 sec = 0.0; - psTime *outTime = NULL; - psTime *tempTime = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time,NULL); - PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - - // Convert time to TAI if necessary, but without changing input arguments - if(time->type == PS_TIME_UTC) { - tempTime = psTimeAlloc(PS_TIME_UTC); - tempTime->sec = time->sec; - tempTime->nsec = time->nsec; - tempTime = psTimeConvert(tempTime, PS_TIME_TAI); - outTime = psTimeAlloc(PS_TIME_TAI); - } else { - tempTime = psMemIncrRefCounter((psTime*)time); - outTime = psTimeAlloc(time->type); - } - - // Create output time - sec = delta + (psF64)tempTime->sec + (psF64)tempTime->nsec/1e9; - PS_ASSERT_LONG_WITHIN_RANGE((psS64)sec,0,PS_MAX_S64,outTime); - outTime->sec = (psS64)sec; - outTime->nsec = (psU32)((sec - (psF64)outTime->sec)*1e9); - - // Error check - PS_ASSERT_INT_WITHIN_RANGE(outTime->nsec,0,(psU32)((1e9)-1),outTime); - - // Convert result to same time type as input - if(time->type == PS_TIME_UTC) { - outTime = psTimeConvert(outTime, PS_TIME_UTC); - } - - psFree(tempTime); - - return outTime; -} - -double psTimeDelta(const psTime *time1, const psTime *time2) -{ - psF64 out = 0.0; - psF64 uSec1 = 0.0; - psF64 uSec2 = 0.0; - psTime *tempTime1 = NULL; - psTime *tempTime2 = NULL; - - // Error checks - PS_ASSERT_PTR_NON_NULL(time1,0.0); - PS_ASSERT_INT_WITHIN_RANGE(time1->nsec,0,(psU32)((1e9)-1),0.0); - PS_ASSERT_PTR_NON_NULL(time2,0.0); - PS_ASSERT_INT_WITHIN_RANGE(time2->nsec,0,(psU32)((1e9)-1),0.0); - - // Verify both times of the same type - if(time1->type != time2->type) { - psError(PS_ERR_BAD_PARAMETER_VALUE,true,PS_ERRORTEXT_psTime_TYPE_INCORRECT,time1->type); - return out; - } - - // Convert time to TAI if necessary, but without changing input arguments - if(time1->type == PS_TIME_UTC) { - tempTime1 = psTimeAlloc(PS_TIME_UTC); - tempTime1->sec = time1->sec; - tempTime1->nsec = time1->nsec; - tempTime1 = psTimeConvert(tempTime1, PS_TIME_TAI); - } else { - tempTime1 = psMemIncrRefCounter((psTime*)time1); - } - if(time2->type == PS_TIME_UTC) { - tempTime2 = psTimeAlloc(PS_TIME_UTC); - tempTime2->sec = time2->sec; - tempTime2->nsec = time2->nsec; - tempTime2 = psTimeConvert(tempTime2, PS_TIME_TAI); - } else { - tempTime2 = psMemIncrRefCounter((psTime*)time2); - } - - uSec1 = tempTime1->sec >= 0 ? 1.0 : -1.0; - uSec1 = uSec1*tempTime1->nsec/1e9; - uSec2 = tempTime2->sec >= 0 ? 1.0 : -1.0; - uSec2 = uSec2*tempTime2->nsec/1e9; - out = (tempTime1->sec-tempTime2->sec) + (uSec1-uSec2); - - psFree(tempTime1); - psFree(tempTime2); - - return out; -} - Index: trunk/psLib/src/astronomy/psTime.h =================================================================== --- trunk/psLib/src/astronomy/psTime.h (revision 4540) +++ (revision ) @@ -1,349 +1,0 @@ -/** @file psTime.h - * - * @brief Definitions for time, time utilities, and conversion functions for use - * with psLib astronomy functions. - * - * A collection of functions are required by psLib to manipulate time data. These - * functions primarily consist of conversions between specific time formats. They - * use the UNIX timeval time system as the base upon which International Atomic - * Time (TAI) and Universal Time Coordinated (UTC) are calculated. - * - * @author Ross Harman, MHPCC - * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-12 19:12:00 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PSTIME_H -#define PSTIME_H - -#include -#include -#include - -#include "psType.h" -// N.B. inclusion of psCoord.h was done to after the typedefs to handle cross-dependency of typedefs - -/// @addtogroup Time -/// @{ - -/** Time type. - * - * Enumeration for psTime types, TAI or UTC time. - */ -typedef enum { - PS_TIME_TAI, ///< Temps Atomique International (TAI) time (time with leapseconds) - PS_TIME_UTC, ///< Universal Time Coordinated (UTC) time (time without leapseconds) - PS_TIME_UT1, ///< Universal Time corrected for polar motion - PS_TIME_TT, ///< Terrestrial Time -} psTimeType; - -/** Time Bulletin type - * - * Enumeration for psTimeBulletin type, A or B. - */ -typedef enum { - PS_IERS_A, ///< IERS Bulletin A - PS_IERS_B, ///< IERS Bulletin B -} psTimeBulletin; - -/** Definition of psTime. - * - * The psTime struct is used by psLib to represent time values critical to - * astronomical calculations. This structure represents a time which is - * equivalent to TAI (International Atomic Time) and is measured in both - * seconds and microseconds. - */ -typedef struct psTime -{ - psS64 sec; ///< Seconds since epoch, Jan 1, 1970. - psU32 nsec; ///< Nanoseconds since last second. - psBool leapsecond; ///< if time falls on UTC leapsecond - psTimeType type; ///< Type of time. -} -psTime; - -#include "psCoord.h" -#include "psImage.h" - -/** Initialize time data. - * - * Reads config and data files associated with various time conversions. - * - * @return bool: True for success, false for failure. - */ -psBool p_psTimeInit( - const char *fileName ///< File name containing config/data info -); - -/** Free memory persistant time data. - * - * Frees time data to be held in memory until the end of successful program execution. - * - * @return void: void. - */ -psBool p_psTimeFinalize(void); - -/** Allocate time struct. - * - * Allocates an empty time struct. User must specify the psTimeType - * (PS_TIME_TAI or PS_TIME_UTC) in the argument. The seconds and microseconds members - * of the struct are set to zero. - * - * @return psTime*: Struct with empty time. - */ -psTime* psTimeAlloc( - psTimeType type ///< Type of time to create (UTC or TAI). -); - -/** Get current time. - * - * Gets current time from the system clock. User must specify the psTimeType - * (PS_TIME_TAI or PS_TIME_UTC) in the argument. - * - * @return psTime*: Struct with current time. - */ -psTime* psTimeGetNow( - psTimeType type ///< Type of time to get (UTC or TAI). -); - -/** Convert psTime to UTC or TAI time. - * - * Converts psTime to UTC or TAI time based on the psTimeType argument. - * - * @return psTime*: Pointer to psTime. - */ -psTime* psTimeConvert( - psTime *time, ///< Time to be converted. - psTimeType type ///< Type to be converted to. -); - -/** Convert psTime to Local Mean Sidereal Time (LMST). - * - * Converts psTime at the given longitude to LMST time. If the input time is not - * in UTC format, then it is converted. - * - * @return double: LST Time. - */ -double psTimeToLMST( - psTime *time, ///< psTime to be converted. - double longitude ///< Longitude. -); - -/** Determine UT1 - UTC from table lookup. - * - * This function is necessary to for various SLALIB functions. - * - * @return double: Time difference. - */ -double psTimeGetUT1Delta( - const psTime *time, ///< psTime to be looked up. - psTimeBulletin bulletin ///< IERS bulletin to use -); - -/** Determine TAI - UTC from table lookup. - * - * This function is necessary to for various psTime functions. - * - * @return psF64: Time difference. - */ -psF64 p_psTimeGetTAIDelta( - const psTime *time ///< psTime to be looked up. -); - -/** Determine polar coordinates at a given time. - * - * Determines the orientation of the polar axis at the given time. - * - * @return psSphere*: Spherical coordinates of Earth's polar axias. - */ -psSphere* p_psTimeGetPoleCoords( - const psTime *time ///< psTime determine polar orientation. -); - -/** Calculate the number of leapseconds between two times. - * - * Calculates the number of leapseconds between two times. - * - * @return long: leapseconds added between given times - */ -long psTimeLeapSecondDelta( - const psTime* time1, ///< First input time. - const psTime* time2 ///< Second input time. -); - -/** Determine if UTC time is a leapsecond. - * - * Determines if the specified UTC time is a valid leapsecond. - * - * @return bool: valid leap second - */ -bool psTimeIsLeapSecond( - const psTime* utc ///< UTC to verify if leap second -); - -/** Convert psTime to Julian date time. - * - * Converts psTime to Julian date (JD) time. This function does not add or - * subtract leapseconds. - * - * @return double: Julian Date (JD) time. - */ -double psTimeToJD( - const psTime* time ///< Input time to be converted. -); -/** Convert psTime to modified Julian date time. - * - * Converts psTime to modified Julian date (MJD) time. This function does not - * add or subtract leapseconds. - * - * @return double: Modified Julian Days (MJD) time. - */ -double psTimeToMJD( - const psTime* time ///< Input time to be converted. -); - -/** Convert psTime to ISO8601 formatted string. - * - * Converts psTime to a null terminated string in the form of YYYY-MM-DDThh:mm:ss.sss. - * This function does not add or subtract leapseconds. - * - * @return psString: Pointer null terminated array of chars in ISO time. - */ -psString psTimeToISO( - const psTime* time ///< Input time to be converted. -); - -/** Convert psTime to timeval time. - * - * Converts psTime to timeval time. This function does not add or subtract leapseconds. - * - * @return timeval*: timeval struct time. - */ -struct timeval* psTimeToTimeval( - const psTime* time ///< Input time to be converted. - ); - -/* - * Convert psTime to tm time. - * - * Converts psTime to tm time. This function is based on a Perl algorithm availble - * in the Pan-STARRS Image processing Algorithm Design Description (ADD). This function - * does not add or subtract leapseconds. - * - * @return tm: tm struct time. - * -struct tm* p_psTimeToTM( - const psTime *time ///< Input time to be converted. - ); -*/ -/** Convert JD to psTime. - * - * Converts JD time to psTime. This function does not add or subtract leapseconds. - * - * @return psTime: time. - */ -psTime* psTimeFromJD( - double jd ///< Input time to be converted. -); - -/** Convert MJD to psTime. - * - * Converts MJD time to psTime. This function does not add or subtract leapseconds. - * - * @return psTime: time. - */ -psTime* psTimeFromMJD( - double mjd ///< Input time to be converted. -); - -/** Convert ISO to psTime. - * - * Converts ISO time to psTime. This function does not add or subtract leapseconds. - * - * @return psTime*: time - */ -psTime* psTimeFromISO( - const char* input ///< Input time to be converted. -); - -/** Convert timeval to psTime. - * - * Converts timeval time to psTime. This function does not add or subtract leapseconds. - * - * @return psTime*: time. - */ -psTime* psTimeFromTimeval( - const struct timeval *input ///< Input time to be converted. -); - -/** Convert Terrestrial Time to psTime - * - * Converts Terrestial Time to psTime. This function assumes resultant time is of type TT. - * - * @return psTime*: time (TT) - */ -psTime* psTimeFromTT( - psS64 sec, ///< Input terrestrial time in seconds - psU32 nsec ///< Input terrestrial time fraction of seconds (nanoseconds) -); - -/** Convert UTC time to psTime - * - * Converts UTC time to psTime. It will verify if time specified is a leapsecond. - * - * @return psTime*: time (UTC)time - */ -psTime* psTimeFromUTC( - psS64 sec, ///< Input time in seconds - psU32 nsec, ///< Input time fraction of seconds (nanoseconds) - bool leapsecond ///< Input time is a leapsecond -); - -/** Convert tm time to psTime. - * - * Converts tm time to psTime. This function is based on a Perl algorithm availble - * in the Pan-STARRS Image processing Algorithm Design Description (ADD). This function - * does not add or subtract leapseconds. - * - * @return psTime*: time. - */ -psTime* p_psTimeFromTM( - const struct tm *time ///< Input time to be converted. -); - -/** Adds delta to time. Result is in TAI time. - * - * Adds delta to time. Input time is converted to TAI format if necessary. - * - * @return psTime*: time. - */ -psTime* psTimeMath( - const psTime *time, ///< Time. - double delta ///< Time delta. -); - -/** Determine difference between two times. Result is in TAI time. - * - * Determine difference between two times. Input times are converted to TAI format if necessary. - * - * @return double: Time difference. - */ -double psTimeDelta( - const psTime *time1, ///< First time. - const psTime *time2 ///< Second time. -); - -/** Get the filename of the psLib configuration file. - * - * @return char* If a PS_CONFIG_FILE environment variable exists, - * that is returned, otherwise the default location - * dependent on the installation location. - */ -char* p_psGetConfigFileName(); - -/// @} - -#endif // #ifndef PSTIME_H