Index: /trunk/psastro/Makefile =================================================================== --- /trunk/psastro/Makefile (revision 5564) +++ /trunk/psastro/Makefile (revision 5565) @@ -27,4 +27,5 @@ $(SRC)/psastroArguments.$(ARCH).o \ $(SRC)/psastroUtils.$(ARCH).o \ +$(SRC)/psModUtils.$(ARCH).o \ $(SRC)/psLibUtils.$(ARCH).o \ $(SRC)/pmAstrom.$(ARCH).o \ @@ -34,8 +35,8 @@ $(SRC)/psastro-mktest.$(ARCH).o \ $(SRC)/psastroArguments.$(ARCH).o \ +$(SRC)/psastroIO.$(ARCH).o \ +$(SRC)/psModUtils.$(ARCH).o \ $(SRC)/psLibUtils.$(ARCH).o \ $(SRC)/pmAstrom.$(ARCH).o - -# $(SRC)/psModUtils.$(ARCH).o \ psastro: $(BIN)/psastro.$(ARCH) @@ -48,4 +49,6 @@ INSTALL = psastro psastro-mktest + +all: psastro.install psastro-mktest.install # dependancy rules for binary code ######################### Index: /trunk/psastro/doc/mktest.txt =================================================================== --- /trunk/psastro/doc/mktest.txt (revision 5565) +++ /trunk/psastro/doc/mktest.txt (revision 5565) @@ -0,0 +1,23 @@ + +the program, psastro-mktest, will generate a set of fake test data for +psastro. the config file 'psastro-mktest.conf' gives an example setup +for this program. It generates 4 output files: + +- cmpfile : this is a basic elixir-style cmp file (stars plus header) + which can be taken as input to psastro + +- catfile : this is a basic text reference catalog, which the current + psastro loads as a reference + +- rawfile : this is a text listing of all astrometry stages for the + generated stars: readout, cell, chip, fpa, tpa, sky, with the last + two columns being the magnitudes. this file is generated by + applying the transformations starting at the readout and going + step-by-step to the sky. + +- reffile : this is the inverse of the preceeding file, with the + resulting sources transformed from the sky coordinates down + step-by-step to the readout pixels. If the transformations are + functioning, these two files (rawfile & reffile) should be identical + within floating point errors. + Index: /trunk/psastro/src/pmAstrom.c =================================================================== --- /trunk/psastro/src/pmAstrom.c (revision 5564) +++ /trunk/psastro/src/pmAstrom.c (revision 5565) @@ -13,102 +13,16 @@ } -psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout) { - - if (cellpix == NULL) { - cellpix = psPlaneAlloc (); - } - - cellpix->x = readpix->x*readout->colBins + readout->col0; - cellpix->y = readpix->y*readout->rowBins + readout->row0; - - return (cellpix); -} - -psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout) { - - if (readpix == NULL) { - readpix = psPlaneAlloc (); - } - - readpix->x = (cellpix->x - readout->col0) / readout->colBins; - readpix->y = (cellpix->y - readout->row0) / readout->rowBins; - - return (readpix); -} - -psPlane* psCoordChipToCell_EAM(psPlane* cellCoord, - const psPlane* chipCoord, - const pmCell* cell) +// sort by mag (descending) +int pmAstromObjSortByMag (const void **a, const void **b) { - PS_ASSERT_PTR_NON_NULL(chipCoord, NULL); - PS_ASSERT_PTR_NON_NULL(cell, NULL); - PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); - - // XXX EAM : why was this being done? - // pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent); - PS_ASSERT_PTR_NON_NULL(cell->toChip, NULL); - psPlaneTransform *tmpChipToCell = p_psPlaneTransformLinearInvert(cell->toChip); - PS_ASSERT_PTR_NON_NULL(tmpChipToCell, NULL); - cellCoord = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord); - psFree(tmpChipToCell); - return(cellCoord); -} - -bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky) { - - bool status; - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; - for (int k = 0; k < cell->readouts->n; k++) { - pmReadout *readout = cell->readouts->data[k]; - psArray *stars = psMetadataLookupPtr (&status, readout->analysis, starlist); - if (toSky) { - psastroProjectRawstars (stars, readout); - } else { - psastroProjectRefstars (stars, readout); - } - } - } - } - return true; -} - -bool psastroProjectRawstars (psArray *stars, pmReadout *readout) { - - pmCell *cell = readout->parent; - pmChip *chip = cell->parent; - pmFPA *fpa = chip->parent; - - for (int i = 0; i < stars->n; i++) { - pmAstromObj *star = stars->data[i]; - psCoordReadoutToCell (&star->cell, &star->pix, readout); - psCoordCellToChip (&star->chip, &star->cell, cell); - psCoordChipToFPA (&star->FP, &star->chip, chip); - psCoordFPAToTP (&star->TP, &star->FP, 0.0, 0.0, fpa); - psCoordTPToSky (&star->sky, &star->TP, fpa->projection); - } - return true; -} - -bool psastroProjectRefstars (psArray *stars, pmReadout *readout) { - - pmCell *cell = readout->parent; - pmChip *chip = cell->parent; - pmFPA *fpa = chip->parent; - - for (int i = 0; i < stars->n; i++) { - pmAstromObj *star = stars->data[i]; - psCoordSkyToTP (&star->TP, &star->sky, fpa->projection); - psCoordTPToFPA (&star->FP, &star->TP, 0.0, 0.0, fpa); - psCoordFPAToChip (&star->chip, &star->FP, chip); - psCoordChipToCell_EAM (&star->cell, &star->chip, cell); - psCoordCellToReadout (&star->pix, &star->cell, readout); - } - return true; -} - + pmAstromObj *A = *(pmAstromObj **)a; + pmAstromObj *B = *(pmAstromObj **)b; + + psF32 diff = A->Mag - B->Mag; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) { Index: /trunk/psastro/src/pmAstrom.h =================================================================== --- /trunk/psastro/src/pmAstrom.h (revision 5564) +++ /trunk/psastro/src/pmAstrom.h (revision 5565) @@ -29,4 +29,5 @@ } pmAstromStats; +/* in pmAstrom.c */ pmAstromObj *pmAstromObjAlloc (void); pmAstromObj *pmAstromObjCopy (pmAstromObj *old); @@ -35,16 +36,9 @@ psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle); psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config); +int pmAstromObjSortByFPX (const void **a, const void **b); +int pmAstromObjSortByMag (const void **a, const void **b); -int pmAstromObjSortByFPX (const void **a, const void **b); - -bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky); -bool psastroProjectRawstars (psArray *stars, pmReadout *readout); -bool psastroProjectRefstars (psArray *stars, pmReadout *readout); - -psPlane* psCoordChipToCell_EAM( - psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created - const psPlane* in, ///< the Chip coordinate - const pmCell* cell ///< the cell of interest -); - -psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform); +/* in pmAstromGrid.c */ +pmAstromStats pmAstromGridAngle (psArray *st1, psArray *st2, psMetadata *config); +psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat); +pmAstromStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config); Index: /trunk/psastro/src/psLibUtils.c =================================================================== --- /trunk/psastro/src/psLibUtils.c (revision 5564) +++ /trunk/psastro/src/psLibUtils.c (revision 5565) @@ -171,2 +171,34 @@ return(out); } + +// returns the rotation term, forcing positive parity +double psPlaneTransformGetRotation (psPlaneTransform *map) { + + if (map->x->nX < 1) return 0; + if (map->x->nY < 1) return 0; + + if (map->y->nX < 1) return 0; + if (map->y->nY < 1) return 0; + + double pc1_1 = map->x->coeff[1][0]; + double pc1_2 = map->x->coeff[0][1]; + double pc2_1 = map->y->coeff[1][0]; + double pc2_2 = map->y->coeff[0][1]; + + double px = SIGN (pc1_1); + double py = SIGN (pc2_2); + + // both x and y terms imply an angle. take the average + double t1 = -atan2 (px*pc1_2, px*pc1_1); + double t2 = +atan2 (py*pc2_1, py*pc2_2); + + // careful near -pi,+pi boundary... + if (t1 - t2 > M_PI/2) t2 += 2*M_PI; + if (t2 - t1 > M_PI/2) t1 += 2*M_PI; + + double theta = 0.5*(t1 + t2); + while (theta < M_PI) theta += 2*M_PI; + while (theta > M_PI) theta -= 2*M_PI; + + return (theta); +} Index: /trunk/psastro/src/psLibUtils.h =================================================================== --- /trunk/psastro/src/psLibUtils.h (revision 5564) +++ /trunk/psastro/src/psLibUtils.h (revision 5565) @@ -1,11 +1,4 @@ # ifndef PS_LIB_UTILS # define PS_LIB_UTILS - -// structure to carry a dynamic string -typedef struct { - int NLINE; - int Nline; - char *line; -} psLine; # define psMemCopy(A)(psMemIncrRefCounter((A))) @@ -13,12 +6,17 @@ # define RAD_DEG 0.017453292519943 +# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) + // psLib extra utilities psS32 psLogArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) psS32 psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -// psLine functions -- keep out for now? -psLine *psLineAlloc (int Nline); -bool psLineInit (psLine *line); -bool psLineAdd (psLine *line, char *format, ...); +psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input); +psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input); +psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input); +psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input); +psProjection *psProjectionCopy (psProjection *input); +psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort); +double psPlaneTransformGetRotation (psPlaneTransform *map); # endif Index: /trunk/psastro/src/psModUtils.c =================================================================== --- /trunk/psastro/src/psModUtils.c (revision 5564) +++ /trunk/psastro/src/psModUtils.c (revision 5565) @@ -1,116 +1,146 @@ # include "psastro.h" -// these functions temporarily replace existing psModules code +psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout) { -// template sample alloc/free pair: -# if (0) -static void psFooFree (psFoo *foo) { + if (cellpix == NULL) { + cellpix = psPlaneAlloc (); + } - if (foo == NULL) return; - return; -} -psFoo *psFooAlloc (int i1, int i2) { + cellpix->x = readpix->x*readout->colBins + readout->col0; + cellpix->y = readpix->y*readout->rowBins + readout->row0; - psFoo *foo = psAlloc (sizeof(psFoo)); - psMemSetDeallocator(foo, (psFreeFunc) psFooFree); - - return (foo); -} -# endif - -// pmFPA -static void pmFPAFree (pmFPA *fpa) { - - if (fpa == NULL) return; - - psFree (fpa->toSky); - psFree (fpa->toTPA); - psFree (fpa->fromTPA); - psFree (fpa->chips); - - return; + return (cellpix); } -pmFPA *pmFPAAlloc (void) { +psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout) { - pmFPA *fpa = psAlloc (sizeof(pmFPA)); - psMemSetDeallocator(fpa, (psFreeFunc) pmFPAFree); + if (readpix == NULL) { + readpix = psPlaneAlloc (); + } - fpa->toSky = NULL; - fpa->toTPA = NULL; - fpa->fromTPA = NULL; - fpa->chips = NULL; + readpix->x = (cellpix->x - readout->col0) / readout->colBins; + readpix->y = (cellpix->y - readout->row0) / readout->rowBins; - return (fpa); + return (readpix); } -// pmChip -static void pmChipFree (pmChip *chip) { +psPlane* psCoordChipToCell_EAM(psPlane* cellCoord, + const psPlane* chipCoord, + const pmCell* cell) +{ + PS_ASSERT_PTR_NON_NULL(chipCoord, NULL); + PS_ASSERT_PTR_NON_NULL(cell, NULL); + PS_ASSERT_PTR_NON_NULL(cell->parent, NULL); - if (chip == NULL) return; - - psFree (chip->toFPA); - psFree (chip->fromFPA); - psFree (chip->cells); - - return; + // XXX EAM : why was this being done? + // pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent); + PS_ASSERT_PTR_NON_NULL(cell->toChip, NULL); + psPlaneTransform *tmpChipToCell = p_psPlaneTransformLinearInvert(cell->toChip); + PS_ASSERT_PTR_NON_NULL(tmpChipToCell, NULL); + cellCoord = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord); + psFree(tmpChipToCell); + return(cellCoord); } -pmChip *pmChipAlloc (void) { +bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky) { - pmChip *chip = psAlloc (sizeof(pmChip)); - psMemSetDeallocator(chip, (psFreeFunc) pmChipFree); + bool status; - chip->toFPA = NULL; - chip->fromFPA = NULL; - chip->cells = NULL; + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + for (int j = 0; j < chip->cells->n; j++) { + pmCell *cell = chip->cells->data[j]; + for (int k = 0; k < cell->readouts->n; k++) { + pmReadout *readout = cell->readouts->data[k]; + psArray *stars = psMetadataLookupPtr (&status, readout->analysis, starlist); + if (toSky) { + psastroProjectRawstars (stars, readout); + } else { + psastroProjectRefstars (stars, readout); + } + } + } + } + return true; +} + +bool psastroProjectRawstars (psArray *stars, pmReadout *readout) { - return (chip); + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + pmFPA *fpa = chip->parent; + + for (int i = 0; i < stars->n; i++) { + pmAstromObj *star = stars->data[i]; + psCoordReadoutToCell (&star->cell, &star->pix, readout); + psCoordCellToChip (&star->chip, &star->cell, cell); + psCoordChipToFPA (&star->FP, &star->chip, chip); + psCoordFPAToTP (&star->TP, &star->FP, 0.0, 0.0, fpa); + psCoordTPToSky (&star->sky, &star->TP, fpa->projection); + } + return true; +} + +bool psastroProjectRefstars (psArray *stars, pmReadout *readout) { + + pmCell *cell = readout->parent; + pmChip *chip = cell->parent; + pmFPA *fpa = chip->parent; + + for (int i = 0; i < stars->n; i++) { + pmAstromObj *star = stars->data[i]; + psCoordSkyToTP (&star->TP, &star->sky, fpa->projection); + psCoordTPToFPA (&star->FP, &star->TP, 0.0, 0.0, fpa); + psCoordFPAToChip (&star->chip, &star->FP, chip); + psCoordChipToCell_EAM (&star->cell, &star->chip, cell); + psCoordCellToReadout (&star->pix, &star->cell, readout); + } + return true; +} + +pmFPA *pmFPACopyAstrom (pmFPA *inFPA) { + + pmFPA *fpa = pmFPAAlloc (NULL, NULL); + + // copy FPA astrometry data + fpa->toSky = psProjectionCopy (inFPA->toSky); + fpa->toTPA = psPlaneDistortCopy (inFPA->toTPA); + fpa->fromTPA = psPlaneDistortCopy (inFPA->fromTPA); + + psArrayRealloc (fpa->chips, inFPA->chips->n); + for (int i = 0; i < inFPA->chips->n; i++) { + pmChip *inChip = inFPA->chips->data[i]; + pmChip *chip = pmChipAlloc (fpa); + + // copy Chip astrometry data + chip->toFPA = psPlaneTransformCopy (inChip->toFPA); + chip->fromFPA = psPlaneTransformCopy (inChip->fromFPA); + + psArrayRealloc (chip->cells, inChip->cells->n); + for (int j = 0; j < inChip->cells->n; j++) { + pmCell *inCell = inChip->cells->data[j]; + pmCell *cell = pmCellAlloc (chip); + + // cell.header is a view on inCell.header + cell->header = psMemCopy (inCell->header); + + // copy Cell astrometry data + cell->toChip = psPlaneTransformCopy (inCell->toChip); + + psArrayRealloc (cell->readouts, inCell->readouts->n); + for (int k = 0; k < inCell->readouts->n; k++) { + pmReadout *inReadout = inCell->readouts->data[k]; + pmReadout *readout = pmReadoutAlloc (cell); + + // copy Readout data + *readout = *inReadout; + + cell->readouts->data[k] = readout; + } + chip->cells->data[j] = cell; + } + fpa->chips->data[i] = chip; + } + return (fpa); } -// pmCell -static void pmCellFree (pmCell *cell) { - - if (cell == NULL) return; - - psFree (cell->toChip); - psFree (cell->readouts); - - return; -} - -pmCell *pmCellAlloc (void) { - - pmCell *cell = psAlloc (sizeof(pmCell)); - psMemSetDeallocator(cell, (psFreeFunc) pmCellFree); - - cell->toChip = NULL; - cell->readouts = NULL; - - return (cell); -} - -// pmReadout -static void pmReadoutFree (pmReadout *readout) { - - if (readout == NULL) return; - - psFree (readout->stars); - - return; -} - -pmReadout *pmReadoutAlloc (void) { - - pmReadout *readout = psAlloc (sizeof(pmReadout)); - psMemSetDeallocator(readout, (psFreeFunc) pmReadoutFree); - - readout->col0 = 0; - readout->row0 = 0; - readout->colBins = 1; - readout->rowBins = 1; - readout->stars = NULL; - - return (readout); -} - Index: /trunk/psastro/src/psastro-mktest.c =================================================================== --- /trunk/psastro/src/psastro-mktest.c (revision 5564) +++ /trunk/psastro/src/psastro-mktest.c (revision 5565) @@ -1,7 +1,3 @@ # include "psastro.h" - -bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); -bool testWriteRef (char *filename, psArray *sources); -bool testWriteRaw (char *filename, psArray *sources); int main (int argc, char **argv) { @@ -173,92 +169,2 @@ exit (0); } - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteCMP (psMetadata *header, char *filename, psArray *sources) { - - int i; - psMetadataItem *mdi; - - psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "foo", sources->n); - - // create file, write-out header - unlink (filename); - psFits *fits = psFitsAlloc (filename); - - // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 - mdi = psMetadataLookup (header, "NAXIS"); - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; - - psFitsWriteHeader (header, fits); - psFree (fits); - - // re-open, add data to end of file - FILE *f = fopen (filename, "a+"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - fseek (f, 0, SEEK_END); - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f\n", - star->pix.x, star->pix.y, star->Mag, (int)(star->dMag*1000), 1, 2.5, 11.0, 12.0, 2.0, 0.5, 15.0); - } - fclose (f); - return true; -} - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteRef (char *filename, psArray *sources) { - - int i; - - // re-open, add data to end of file - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%10.7f %10.7f %6.3f %6.3f\n", star->sky.r, star->sky.d+10*RAD_DEG/3600, 10.0, 0.05); - } - fclose (f); - return true; -} - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteRaw (char *filename, psArray *sources) { - - int i; - - // re-open, add data to end of file - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", - star->pix.x, star->pix.y, - star->cell.x, star->cell.y, - star->chip.x, star->chip.y, - star->FP.x, star->FP.y, - star->TP.x, star->TP.y, - star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, - star->Mag, star->dMag); - } - fclose (f); - return true; -} Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 5564) +++ /trunk/psastro/src/psastro.h (revision 5565) @@ -11,51 +11,4 @@ # define toSky projection -# if (0) -typedef struct -{ - psProjection *toSky; ///< Projection from tangent plane to sky - psPlaneDistort *toTPA; ///< Transformation from focal plane to tangent plane - psPlaneDistort *fromTPA; ///< Transformation from focal plane to tangent plane - psArray *chips; ///< The chips -} -pmFPA; - -typedef struct -{ - // Astrometric transformations - psPlaneTransform *toFPA; ///< Transformation from chip to FPA coordinates - psPlaneTransform *fromFPA; ///< Transformation from FPA to chip coordinates - psArray *cells; ///< The cells - pmFPA *fpa; -} -pmChip; - -typedef struct -{ - // Astrometric transformations - psPlaneTransform *toChip; ///< Transformations from cell to chip coordinates - psArray *readouts; ///< The readouts (referred to by number) - psMetadata *header; ///< header always corresponds to a cell - pmChip *chip; -} -pmCell; - -typedef struct -{ - // Position on the cell - int col0; ///< Offset from the left of chip. - int row0; ///< Offset from the bottom of chip. - int colBins; ///< Amount of binning in x-dimension - int rowBins; ///< Amount of binning in y-dimension - psArray *stars; ///< sources detected / measured on readout - pmCell *cell; -} -pmReadout; -# endif - -pmAstromStats pmAstromGridAngle (psArray *st1, psArray *st2, psMetadata *config); -psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat); -pmAstromStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config); -int pmAstromObjSortByMag (const void **a, const void **b); bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) ; pmFPA *pmFPACopyAstrom (pmFPA *inFPA); @@ -70,22 +23,22 @@ psMetadata *testArguments (int *argc, char **argv); -# if (0) -pmFPA *pmFPAAlloc (void); -pmChip *pmChipAlloc (void); -pmCell *pmCellAlloc (void); -pmReadout *pmReadoutAlloc (void); -# endif - -# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) - -psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input); -psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input); -psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input); -psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input); -psProjection *psProjectionCopy (psProjection *input); psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout); psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout); -double psPlaneTransformGetRotation (psPlaneTransform *map); +psPlane* psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell); -psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort); -psPlaneTransform *psPlaneTransformLinearInvert(psPlaneTransform *transform); +bool testWriteCMP (psMetadata *header, char *filename, psArray *sources); +bool testWriteRef (char *filename, psArray *sources); +bool testWriteRaw (char *filename, psArray *sources); + +bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky); +bool psastroProjectRawstars (psArray *stars, pmReadout *readout); +bool psastroProjectRefstars (psArray *stars, pmReadout *readout); + +psPlane* psCoordChipToCell_EAM( + psPlane* out, ///< a plane struct to recycle. If NULL, a new struct is created + const psPlane* in, ///< the Chip coordinate + const pmCell* cell ///< the cell of interest +); + +psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform); +pmFPA *pmFPACopyAstrom (pmFPA *inFPA); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 5564) +++ /trunk/psastro/src/psastroArguments.c (revision 5565) @@ -46,5 +46,5 @@ static void usage_test (void) { - fprintf (stderr, "USAGE: psastro-mktest (config) (cmpfile) (reflist) (rawfile) (rawfile)\n"); + fprintf (stderr, "USAGE: psastro-mktest (config) (cmpfile) (catfile) (rawfile) (reffile)\n"); exit (2); } Index: /trunk/psastro/src/psastroIO.c =================================================================== --- /trunk/psastro/src/psastroIO.c (revision 5564) +++ /trunk/psastro/src/psastroIO.c (revision 5565) @@ -124,2 +124,107 @@ } +// elixir-style pseudo FITS table (header + ascii list) +bool testWriteCMP (psMetadata *header, char *filename, psArray *sources) { + + int i; + psMetadataItem *mdi; + + psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "foo", sources->n); + + // create file, write-out header + unlink (filename); + psFits *fits = psFitsAlloc (filename); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + psFitsWriteHeader (header, fits); + psFree (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f\n", + star->pix.x, star->pix.y, star->Mag, (int)(star->dMag*1000), 1, 2.5, 11.0, 12.0, 2.0, 0.5, 15.0); + } + fclose (f); + return true; +} + +// write out raw objects +bool testWriteRaw (char *filename, psArray *sources) { + + int i; + + // re-open, add data to end of file + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", + star->pix.x, star->pix.y, + star->cell.x, star->cell.y, + star->chip.x, star->chip.y, + star->FP.x, star->FP.y, + star->TP.x, star->TP.y, + star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, + star->Mag, star->dMag); + } + fclose (f); + return true; +} + +// elixir-style pseudo FITS table (header + ascii list) +bool testWriteRef (char *filename, psArray *sources) { + + int i; + + // re-open, add data to end of file + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + + pmAstromObj *star = sources->data[i]; + + fprintf (f, "%10.7f %10.7f %6.3f %6.3f\n", star->sky.r, star->sky.d+10*RAD_DEG/3600, 10.0, 0.05); + } + fclose (f); + return true; +} + +void psastroDumpStars (psArray *sources) { + + for (int i = 0; i < sources->n; i++) { + pmAstromObj *star = sources->data[i]; + + fprintf (stderr, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", + star->pix.x, star->pix.y, + star->cell.x, star->cell.y, + star->chip.x, star->chip.y, + star->FP.x, star->FP.y, + star->TP.x, star->TP.y, + star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, + star->Mag, star->dMag); + } +} Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 5564) +++ /trunk/psastro/src/psastroUtils.c (revision 5565) @@ -1,34 +1,4 @@ # include "psastro.h" # define RENORM 0 - -bool testWriteRaw (char *filename, psArray *sources); - -void psastroDumpStars (psArray *sources) { - - for (int i = 0; i < sources->n; i++) { - pmAstromObj *star = sources->data[i]; - - fprintf (stderr, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", - star->pix.x, star->pix.y, - star->cell.x, star->cell.y, - star->chip.x, star->chip.y, - star->FP.x, star->FP.y, - star->TP.x, star->TP.y, - star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, - star->Mag, star->dMag); - } -} - -// sort by mag (descending) -int pmAstromObjSortByMag (const void **a, const void **b) -{ - pmAstromObj *A = *(pmAstromObj **)a; - pmAstromObj *B = *(pmAstromObj **)b; - - psF32 diff = A->Mag - B->Mag; - if (diff > FLT_EPSILON) return (-1); - if (diff < FLT_EPSILON) return (+1); - return (0); -} bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) { @@ -62,50 +32,4 @@ } -pmFPA *pmFPACopyAstrom (pmFPA *inFPA) { - - pmFPA *fpa = pmFPAAlloc (NULL, NULL); - - // copy FPA astrometry data - fpa->toSky = psProjectionCopy (inFPA->toSky); - fpa->toTPA = psPlaneDistortCopy (inFPA->toTPA); - fpa->fromTPA = psPlaneDistortCopy (inFPA->fromTPA); - - psArrayRealloc (fpa->chips, inFPA->chips->n); - for (int i = 0; i < inFPA->chips->n; i++) { - pmChip *inChip = inFPA->chips->data[i]; - pmChip *chip = pmChipAlloc (fpa); - - // copy Chip astrometry data - chip->toFPA = psPlaneTransformCopy (inChip->toFPA); - chip->fromFPA = psPlaneTransformCopy (inChip->fromFPA); - - psArrayRealloc (chip->cells, inChip->cells->n); - for (int j = 0; j < inChip->cells->n; j++) { - pmCell *inCell = inChip->cells->data[j]; - pmCell *cell = pmCellAlloc (chip); - - // cell.header is a view on inCell.header - cell->header = psMemCopy (inCell->header); - - // copy Cell astrometry data - cell->toChip = psPlaneTransformCopy (inCell->toChip); - - psArrayRealloc (cell->readouts, inCell->readouts->n); - for (int k = 0; k < inCell->readouts->n; k++) { - pmReadout *inReadout = inCell->readouts->data[k]; - pmReadout *readout = pmReadoutAlloc (cell); - - // copy Readout data - *readout = *inReadout; - - cell->readouts->data[k] = readout; - } - chip->cells->data[j] = cell; - } - fpa->chips->data[i] = chip; - } - return (fpa); -} - // measure per-chip astrometry terms bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) { @@ -153,6 +77,6 @@ psastroProjectRefstars (refstars, readout); - testWriteRaw ("ref.inp", refstars); - testWriteRaw ("raw.inp", rawstars); + // testWriteRaw ("ref.inp", refstars); + // testWriteRaw ("raw.inp", rawstars); // fprintf (stderr, "rawstars:\n"); @@ -172,6 +96,6 @@ psastroProjectRefstars (refstars, readout); - testWriteRaw ("ref.dat", refstars); - testWriteRaw ("raw.dat", rawstars); + // testWriteRaw ("ref.dat", refstars); + // testWriteRaw ("raw.dat", rawstars); // use small radius to match stars @@ -284,63 +208,2 @@ return true; } - -// returns the rotation term, forcing positive parity -double psPlaneTransformGetRotation (psPlaneTransform *map) { - - if (map->x->nX < 1) return 0; - if (map->x->nY < 1) return 0; - - if (map->y->nX < 1) return 0; - if (map->y->nY < 1) return 0; - - double pc1_1 = map->x->coeff[1][0]; - double pc1_2 = map->x->coeff[0][1]; - double pc2_1 = map->y->coeff[1][0]; - double pc2_2 = map->y->coeff[0][1]; - - double px = SIGN (pc1_1); - double py = SIGN (pc2_2); - - // both x and y terms imply an angle. take the average - double t1 = -atan2 (px*pc1_2, px*pc1_1); - double t2 = +atan2 (py*pc2_1, py*pc2_2); - - // careful near -pi,+pi boundary... - if (t1 - t2 > M_PI/2) t2 += 2*M_PI; - if (t2 - t1 > M_PI/2) t1 += 2*M_PI; - - double theta = 0.5*(t1 + t2); - while (theta < M_PI) theta += 2*M_PI; - while (theta > M_PI) theta -= 2*M_PI; - - return (theta); -} - -// elixir-style pseudo FITS table (header + ascii list) -bool testWriteRaw (char *filename, psArray *sources) { - - int i; - - // re-open, add data to end of file - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - - pmAstromObj *star = sources->data[i]; - - fprintf (f, "%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %10.6f %10.6f %8.2f %8.2f\n", - star->pix.x, star->pix.y, - star->cell.x, star->cell.y, - star->chip.x, star->chip.y, - star->FP.x, star->FP.y, - star->TP.x, star->TP.y, - star->sky.r*DEG_RAD, star->sky.d*DEG_RAD, - star->Mag, star->dMag); - } - fclose (f); - return true; -}