Index: trunk/psastro/src/Makefile.am =================================================================== --- trunk/psastro/src/Makefile.am (revision 15890) +++ trunk/psastro/src/Makefile.am (revision 15891) @@ -3,5 +3,5 @@ libpsastro_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) -bin_PROGRAMS = psastro psastroModel +bin_PROGRAMS = psastro psastroModel gpcModel psastro_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) @@ -12,4 +12,8 @@ psastroModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS) psastroModel_LDADD = libpsastro.la + +gpcModel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) +gpcModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS) +gpcModel_LDADD = libpsastro.la psastro_SOURCES = \ @@ -29,5 +33,9 @@ psastroModelBoresite.c \ psastroModelFitBoresite.c \ + psastroModelAdjust.c \ + psastroModelDataSave.c \ psastroCleanup.c + +gpcModel_SOURCES = gpcModel.c ## move DataSave to psastro_SOURCES? @@ -48,5 +56,6 @@ psastroLuminosityFunction.c \ psastroRefstarSubset.c \ - psastroEnforceChips.c \ + psastroFixChips.c \ + psastroUseModel.c \ psastroMosaicAstrom.c \ psastroMosaicGradients.c \ Index: trunk/psastro/src/gpcModel.c =================================================================== --- trunk/psastro/src/gpcModel.c (revision 15891) +++ trunk/psastro/src/gpcModel.c (revision 15891) @@ -0,0 +1,165 @@ +# include "psastroStandAlone.h" + +int main (int argc, char **argv) { + + // USAGE: gpcModel (input) (output) + // generate model for gpc based on input model for all but chips + + if (argc != 3) { + fprintf (stderr, "USAGE: gpcModel (input) (output)\n"); + exit (2); + } + + psFits *input = psFitsOpen (argv[1], "r"); + psFits *output = psFitsOpen (argv[2], "w"); + + { // PHU + psMetadata *header = psFitsReadHeader (NULL, input); + psFitsWriteBlank (output, header, ""); + } + + { // CHIPS + + // read initial chips + if (!psFitsMoveExtName (input, "CHIPS")) { + psError(PS_ERR_IO, false, "missing CHIPS extension in astrometry table\n"); + return false; + } + psArray *chips = psFitsReadTable (input); + if (!chips) psAbort("cannot read chips"); + fprintf (stderr, "read %ld rows from CHIPS\n", chips->n); + + psMetadata *header = psMetadataAlloc(); + psMetadataAddStr(header, PS_LIST_TAIL, "COORD", PS_META_REPLACE, "name of this layer", "CHIPS"); + psMetadataAddStr(header, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "FOCAL_PLANE"); + psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region", "RECTANGLE"); + psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent", "POLYNOMIAL"); + + float coeff[2][2]; + char coeffMask[2][2]; + coeff[0][0] = 0.0; + coeff[1][0] = 1.0; + coeff[0][1] = 0.0; + coeff[1][1] = 0.0; + + coeffMask[0][0] = 0; + coeffMask[1][0] = 0; + coeffMask[0][1] = 0; + coeffMask[1][1] = 1; + + psArray *table = psArrayAllocEmpty (1); + for (int i = 0; i < 8; i++) { + for (int j = 0; j < 8; j++) { + if ((i == 0) && (j == 0)) continue; + if ((i == 0) && (j == 7)) continue; + if ((i == 7) && (j == 0)) continue; + if ((i == 7) && (j == 7)) continue; + + float Xo; + float Yo; + if (i < 4) { + Xo = (3 - i)*4970.0 + 60.0; + Yo = (3 - j)*5133.0 + 125.0; + coeff[1][0] = +1.0; + } else { + Xo = (4 - i)*4970.0 - 60.0; + Yo = (4 - j)*5133.0 - 150.0; + coeff[1][0] = -1.0; + } + + for (int ix = 0; ix < 2; ix++) { + for (int iy = 0; iy < 2; iy++) { + + psMetadata *row = psMetadataAlloc (); + + char chipname[80]; + sprintf (chipname, "XY%d%d", i, j); + psMetadataAddStr(row, PS_LIST_TAIL, "SEGMENT", PS_META_REPLACE, "name of this segment", chipname); + psMetadataAddStr(row, PS_LIST_TAIL, "PARENT", PS_META_REPLACE, "next layer up", "FOCAL_PLANE"); + psMetadataAddF32(row, PS_LIST_TAIL, "MINX", PS_META_REPLACE, "range", 0.0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXX", PS_META_REPLACE, "range", 4846.0); + psMetadataAddF32(row, PS_LIST_TAIL, "MINY", PS_META_REPLACE, "range", 0.0); + psMetadataAddF32(row, PS_LIST_TAIL, "MAXY", PS_META_REPLACE, "range", 4868.0); + + psMetadataAddS32(row, PS_LIST_TAIL, "XORDER", PS_META_REPLACE, "", ix); + psMetadataAddS32(row, PS_LIST_TAIL, "YORDER", PS_META_REPLACE, "", iy); + psMetadataAddS32(row, PS_LIST_TAIL, "NXORDER", PS_META_REPLACE, "", 1); + psMetadataAddS32(row, PS_LIST_TAIL, "NYORDER", PS_META_REPLACE, "", 1); + if ((ix == 0) && (iy == 0)) { + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", Xo); + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", Yo); + } else { + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_X", PS_META_REPLACE, "", coeff[ix][iy]); + psMetadataAddF32(row, PS_LIST_TAIL, "POLY_Y", PS_META_REPLACE, "", coeff[iy][ix]); + } + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_X", PS_META_REPLACE, "", 0.0); + psMetadataAddF32(row, PS_LIST_TAIL, "ERROR_Y", PS_META_REPLACE, "", 0.0); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_X", PS_META_REPLACE, "", coeffMask[ix][iy]); + psMetadataAddU8 (row, PS_LIST_TAIL, "MASK_Y", PS_META_REPLACE, "", coeffMask[ix][iy]); + psArrayAdd (table, 100, row); + psFree (row); + } + } + } + } + + if (!psFitsWriteTable (output, header, table, "CHIPS")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + } + + { // FP + if (!psFitsMoveExtName (input, "FP")) { + psError(PS_ERR_IO, false, "missing FP extension in astrometry table\n"); + return false; + } + psMetadata *header = psFitsReadHeader (NULL, input); + psArray *table = psFitsReadTable (input); + if (!table) psAbort("cannot read fp"); + fprintf (stderr, "read %ld rows from FP\n", table->n); + + if (!psFitsWriteTable (output, header, table, "FP")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + } + + { // TP + if (!psFitsMoveExtName (input, "TP")) { + psError(PS_ERR_IO, false, "missing TP extension in astrometry table\n"); + return false; + } + psMetadata *header = psFitsReadHeader (NULL, input); + psArray *table = psFitsReadTable (input); + if (!table) psAbort("cannot read tp"); + fprintf (stderr, "read %ld rows from TP\n", table->n); + if (!psFitsWriteTable (output, header, table, "TP")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(table); + return false; + } + } + + { // SKY + if (!psFitsMoveExtName (input, "SKY")) { + psError(PS_ERR_IO, false, "missing SKY extension in astrometry table\n"); + return false; + } + psMetadata *header = psFitsReadHeader (NULL, input); + psArray *sky = psFitsReadTable (input); + if (!sky) psAbort("cannot read sky"); + fprintf (stderr, "read %ld rows from SKY\n", sky->n); + if (!psFitsWriteTable (output, header, sky, "SKY")) { + psError(PS_ERR_IO, false, "writing sky data\n"); + psFree(sky); + return false; + } + } + + psFitsClose (input); + psFitsClose (output); + exit (0); +} Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 15890) +++ trunk/psastro/src/psastro.h (revision 15891) @@ -80,8 +80,12 @@ int psastroSortByMag (const void *a, const void *b); -bool psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe); +bool psastroFixChips (pmConfig *config, psMetadata *recipe); +bool psastroUseModel (pmConfig *config, psMetadata *recipe); -psArray *psastroReadGetstarCatalog (psFits *fits); -psArray *psastroReadGetstar_PS1_DEV_0 (psFits *fits); +psArray *psastroReadGetstarCatalog (psFits *fits); +psArray *psastroReadGetstar_PS1_DEV_0 (psFits *fits); + +bool psastroAstromGuessSetChip (pmFPA *fpa, pmChip *chip, pmFPAview *view, double pixelScale, bool bilevelAstrometry); +bool psastroAstromGuessSetFPA (pmFPA *fpa, bool *bilevelAstrometry); # endif /* PSASTRO_H */ Index: trunk/psastro/src/psastroAnalysis.c =================================================================== --- trunk/psastro/src/psastroAnalysis.c (revision 15890) +++ trunk/psastro/src/psastroAnalysis.c (revision 15891) @@ -5,4 +5,16 @@ bool status; int nStars; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + if (!psastroUseModel (config, recipe)) { + psError (PSASTRO_ERR_UNKNOWN, false, "failed to set model astrometry\n"); + return false; + } // interpret the available initial astrometric information @@ -26,11 +38,4 @@ if (!psastroChooseRefstars (config, refs)) { psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n"); - return false; - } - - // select the current recipe - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); - if (!recipe) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); return false; } Index: trunk/psastro/src/psastroArguments.c =================================================================== --- trunk/psastro/src/psastroArguments.c (revision 15890) +++ trunk/psastro/src/psastroArguments.c (revision 15891) @@ -43,10 +43,16 @@ psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); } + // no valid header WCS: supply from astrom model + if ((N = psArgumentGet (argc, argv, "-use-astrommodel"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.USE.MODEL", PS_META_REPLACE, "", true); + } // define the reference astrometry file - status = pmConfigFileSetsMD (config->arguments, &argc, argv, "REF.ASTROM", "-refastrom", "-refastromlist"); + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist"); if (status) { // if supplied, assume -fixchips is desired psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); } + // apply mosastro mode? Index: trunk/psastro/src/psastroAstromGuess.c =================================================================== --- trunk/psastro/src/psastroAstromGuess.c (revision 15890) +++ trunk/psastro/src/psastroAstromGuess.c (revision 15891) @@ -32,4 +32,10 @@ } + // have we already supplied the astrometry from the model? + bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); + if (!status) { + useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + } + // select the input data sources pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); @@ -46,36 +52,19 @@ } - pmFPAview *view = pmFPAviewAlloc (0); pmFPA *fpa = input->fpa; // load mosaic-level astrometry? bool bilevelAstrometry = false; - pmHDU *phu = pmFPAviewThisPHU (view, fpa); - if (phu) { - char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); - if (ctype) { - bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } + if (!useModel) { + psastroAstromGuessSetFPA (fpa, &bilevelAstrometry); } - if (bilevelAstrometry) { - pmAstromReadBilevelMosaic (fpa, phu->header); - } + pmFPAview *view = pmFPAviewAlloc (0); while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } - // read WCS data from the corresponding header - pmHDU *hdu = pmFPAviewThisHDU (view, fpa); - if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (chip, hdu->header)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - continue; - } - } else { - if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { - psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); - continue; - } + if (!useModel) { + if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; } @@ -160,2 +149,41 @@ sky (ra, dec) */ + +bool psastroAstromGuessSetChip (pmFPA *fpa, pmChip *chip, pmFPAview *view, double pixelScale, bool bilevelAstrometry) { + + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } + } else { + if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) { + psWarning("Could not get WCS information from header for chip %d, skipping", view->chip); + return false; + } + } + return true; +} + +bool psastroAstromGuessSetFPA (pmFPA *fpa, bool *bilevelAstrometry) { + + pmFPAview *view = pmFPAviewAlloc (0); + pmHDU *phu = pmFPAviewThisPHU (view, fpa); + + *bilevelAstrometry = false; + + // load mosaic-level astrometry? + if (phu) { + char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1"); + if (ctype) { + *bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (*bilevelAstrometry) { + pmAstromReadBilevelMosaic (fpa, phu->header); + } + psFree (view); + return true; +} Index: trunk/psastro/src/psastroDefineFiles.c =================================================================== --- trunk/psastro/src/psastroDefineFiles.c (revision 15890) +++ trunk/psastro/src/psastroDefineFiles.c (revision 15891) @@ -20,20 +20,17 @@ output->save = true; - bool fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); - if (fixChips) { - if (!psastroDefineFile (config, input->fpa, "PSASTRO.REF.ASTROM", "REF.ASTROM", PM_FPA_FILE_ASTROM, PM_DETREND_TYPE_ASTROM)) { - psError (PS_ERR_IO, false, "Can't find a reference astrometry file"); + bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); + if (!status) { + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); + } + bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); + if (!status) { + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + } + if (fixChips || useModel) { + if (!psastroDefineFile (config, input->fpa, "PSASTRO.MODEL", "ASTROM.MODEL", PM_FPA_FILE_ASTROM, PM_DETREND_TYPE_ASTROM)) { + psError (PS_ERR_IO, false, "Can't find an astrometry model file"); return NULL; } - } - - bool saveAstrom = psMetadataLookupBool (&status, recipe, "SAVE.REF.ASTROM"); - if (saveAstrom) { - pmFPAfile *outAstrom = pmFPAfileDefineOutputFromFile (config, input, "PSASTRO.OUT.ASTROM"); - if (!outAstrom) { - psError (PS_ERR_IO, false, "Can't find a reference astrometry file"); - return NULL; - } - outAstrom->save = true; } Index: trunk/psastro/src/psastroEnforceChips.c =================================================================== --- trunk/psastro/src/psastroEnforceChips.c (revision 15890) +++ trunk/psastro/src/psastroEnforceChips.c (revision 15891) @@ -1,14 +1,82 @@ # include "psastroInternal.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ -bool psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe) { +XXX add a function to supply for all entries, as opposed to fixing the poor astrometry ones. +bool psastroEnforceChips (pmConfig *config, psMetadata *recipe) { bool status; + bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); + if (!status) { + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); + } + if (!fixChips) return true; + // identify reference astrometry table // if not defined, correction was not requested; skip step - pmFPAfile *refAstrom = psMetadataLookupPtr (NULL, config->files, "REF.ASTROM"); - if (!refAstrom) return true; + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL"); + if (!astrom) return true; + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + // make sure the astrometry model is loaded + // de-activate all files except PSASTRO.MODEL + // XXX do I need to supply the input->fpa->concepts to the astrom->fpa->concepts? + psFree (astrom->fpa->concepts); + astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts); + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.MODEL"); pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError (PS_ERR_IO, false, "Can't load the astrometry model file"); + return false; + } + + // loop over all chips, replace input astrometry elements with those from astrom + pmChip *obsChip = NULL; + while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // XXX can we do a straight copy? I think so, but we may need to apply a scaling factor based + // on the toSky or toTPA values + obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); + obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); + + // XXX test to write out adjusted header + pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL); + } + + psFree (input->fpa->toSky); + psFree (input->fpa->toTPA); + psFree (input->fpa->fromTPA); + input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + + psMetadata *updates = psMetadataAlloc(); + pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL); + psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); + psFree (updates); + + psFree (view); + return true; +} + +# if (0) // physical pixel scale in microns per pixel @@ -33,6 +101,6 @@ if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } - // set the chip astrometry using the refAstrom file - pmChip *refChip = pmFPAviewThisChip (view, refAstrom->fpa); + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); // bad Astrometry test: ref pixel or angle outside nominal @@ -80,4 +148,5 @@ // XXX for this function to work, I need to: // 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?) -// 2) load the refAstrom fpa in the pmFPAfileIO stage +// 2) load the astrom fpa in the pmFPAfileIO stage // 3) allow for this operation to be optional +# endif Index: trunk/psastro/src/psastroFixChips.c =================================================================== --- trunk/psastro/src/psastroFixChips.c (revision 15891) +++ trunk/psastro/src/psastroFixChips.c (revision 15891) @@ -0,0 +1,112 @@ +# include "psastroInternal.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ + +// XXX I think the 'badAstrom' tests may need to be adjusted: see eg the nominal rotation of +// the chips +bool psastroFixChips (pmConfig *config, psMetadata *recipe) { + + bool status; + + bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS"); + if (!status) { + fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS"); + } + if (!fixChips) return true; + + // identify reference astrometry table + // if not defined, correction was not requested; skip step + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL"); + if (!astrom) return true; + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + // acceptable limits + double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE"); + if (!status) psAbort ("missing recipe value"); + + double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE"); + if (!status) psAbort ("missing recipe value"); + + // make sure the astrometry model is loaded. de-activate all files except PSASTRO.MODEL. + // supply the input concepts so the model is defined using the RA, DEC, POSANGLE of the input + // image. + psFree (astrom->fpa->concepts); + astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts); + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.MODEL"); + + pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError (PS_ERR_IO, false, "Can't load the astrometry model file"); + return false; + } + + // loop over all chips, replace input astrometry elements with those from astrom + pmChip *obsChip = NULL; + while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + // bad Astrometry test: ref pixel or angle outside nominal + + psPlane refPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsCoord, refCoord; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel); + + psPlane offPixel = {0.0, 0.0, 0.0, 0.0}; + psPlane obsOffPt, refOffPt; + + // find location of 0,0 pixel in focal plane coords for this chip + psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel); + + // find location of 0,0 pixel in focal plane coords for ref chip + psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel); + + double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x); + double refAngle = atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x); + + bool badAstrom = false; + badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol; + badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol; + badAstrom |= fabs(obsAngle - refAngle) > angleTol; + + if (!badAstrom) continue; + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // supply astrometry from model + obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); + obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); + } + + psFree (input->fpa->toSky); + psFree (input->fpa->toTPA); + psFree (input->fpa->fromTPA); + input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + + psMetadata *updates = psMetadataAlloc(); + pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL); + psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); + psFree (updates); + + psFree (view); + return true; +} Index: trunk/psastro/src/psastroModel.c =================================================================== --- trunk/psastro/src/psastroModel.c (revision 15890) +++ trunk/psastro/src/psastroModel.c (revision 15891) @@ -31,8 +31,20 @@ // run the full astrometry analysis (chip and/or mosaic) if (!psastroModelAnalysis (config)) { - psErrorStackPrint(stderr, "failure in psastro analysis\n"); + psErrorStackPrint(stderr, "failure in psastro model analysis\n"); exit (1); } + // run the full astrometry analysis (chip and/or mosaic) + if (!psastroModelAdjust (config)) { + psErrorStackPrint(stderr, "failure in psastro model adjust\n"); + exit (1); + } + + // save the model + if (!psastroModelDataSave (config)) { + psErrorStackPrint(stderr, "error saving output data\n"); + exit (1); + } + psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete")); Index: trunk/psastro/src/psastroModelAdjust.c =================================================================== --- trunk/psastro/src/psastroModelAdjust.c (revision 15891) +++ trunk/psastro/src/psastroModelAdjust.c (revision 15891) @@ -0,0 +1,136 @@ +# include "psastroStandAlone.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ + +bool psastroModelAdjust (pmConfig *config) { + + bool status; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL"); + if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL"); + if (!output->fpa) psAbort ("no existing input fpa contains the reference chip"); + + // physical pixel scale in microns per pixel + char *refChipName = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP"); + if (!refChipName) { + psError(PS_ERR_IO, true, "reference chip is missing from recipe"); + return false; + } + + // get reference chip from name + pmChip *refChip = pmConceptsChipFromName (output->fpa, refChipName); + if (!refChip) psAbort ("invalid chip name for reference"); + if (!refChip->toFPA) psAbort ("invalid astrometry for reference chip"); + + psPlane *boreCH = psPlaneAlloc(); + psPlane *boreFP = psPlaneAlloc(); + psPlane *boreTP = psPlaneAlloc(); + psSphere *boreSky = psSphereAlloc(); + + // correct Xo,Yo to Lo,Mo using the ref chip toFPA + // ref chip position of the true boresite center + boreCH->x = psMetadataLookupF32 (&status, output->fpa->concepts, "FPA.BORE.X0"); + boreCH->y = psMetadataLookupF32 (&status, output->fpa->concepts, "FPA.BORE.Y0"); + psPlaneTransformApply (boreFP, refChip->toFPA, boreCH); + + // adjust the reference pixel for all chips + for (int i = 0; i < output->fpa->chips->n; i++) { + pmChip *chip = output->fpa->chips->data[i]; + if (!chip->fromFPA) continue; + // skip the chips without astrometry + + psPlaneTransform *fromFPA = psPlaneTransformSetCenter (NULL, chip->fromFPA, boreFP->x, boreFP->y); + psFree (chip->fromFPA); + chip->fromFPA = fromFPA; + + // invert the new fromFPA transform to get the new toFPA transform + psRegion *region = pmChipPixels (chip); + psFree (chip->toFPA); + chip->toFPA = psPlaneTransformInvert(NULL, chip->fromFPA, *region, 50); + psFree (region); + } + + // XXX we have now shifted the (0,0) pixel of the focal plane to the true boresite from the + // reported boresite. At this point, we need to shift the tangent plane to use the new + // center as well. if the toTPA transform were not linear, we would need to modify fromFPA + // to yield 0,0 at the new boresite location (ie, find Po,Qo = toTPA(Lo,Mo), probably could modify the + // toTPA/fromTPA transforms to use the new ref pixel, but this would only give us a tranf + + // save the old (L,M) = (0,0) TP coordinate + float Po = output->fpa->toTPA->x->coeff[0][0]; + float Qo = output->fpa->toTPA->y->coeff[0][0]; + + // the new toTPA yields the same TP coordinates for FP coordinates offset by Lo,Mo + psPlaneTransform *toTPA = psPlaneTransformSetCenter (NULL, output->fpa->toTPA, boreFP->x, boreFP->y); + psFree (output->fpa->toTPA); + output->fpa->toTPA = toTPA; + + // we are going to shift P,Q to have toTPA(0,0) = Po, Qo. + // find the sky coordinates of the 0,0 pixel for the new frame + boreTP->x = output->fpa->toTPA->x->coeff[0][0] - Po; + boreTP->y = output->fpa->toTPA->y->coeff[0][0] - Qo; + psDeproject (boreSky, boreTP, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate + + // adjust the new TP frame to yield the same old (L,M) = 0,0 TP coordinate: + output->fpa->toTPA->x->coeff[0][0] = Po; + output->fpa->toTPA->y->coeff[0][0] = Qo; + + // set the projection to use the new (P,Q) = (0,0) position + output->fpa->toSky->R = boreSky->r; + output->fpa->toSky->D = boreSky->d; + + psRegion *fpaRegion = pmAstromFPAExtent (output->fpa); + + psFree (output->fpa->fromTPA); + output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50); + + // rotate the chip to FPA transforms to have 0.0 posangle; + // compensate by rotating fpa to TPA transforms + + // get the current posangle of the ref chip + float posangle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle); + + // rotate the chip transforms + for (int i = 0; i < output->fpa->chips->n; i++) { + pmChip *chip = output->fpa->chips->data[i]; + if (!chip->toFPA) continue; + // skip chips without astrometry + + psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posangle); + psFree (chip->toFPA); + chip->toFPA = toFPA; + + // invert the new fromFPA transform to get the new toFPA transform + psRegion *region = pmChipPixels (chip); + psFree (chip->fromFPA); + chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); + psFree (region); + + // XXX temporary output dump + psMetadata *header = chip->hdu->header; + // XXX to make the output single-level, this needs to be in a loop *after* the fromTPA rotation below + // pmAstromWriteWCS (header, output->fpa, chip, NONLIN_TOL); + pmAstromWriteBilevelChip (header, chip, NONLIN_TOL); + } + + psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, output->fpa->fromTPA, posangle); + psFree (output->fpa->fromTPA); + output->fpa->fromTPA = fromTPA; + + psFree (output->fpa->toTPA); + output->fpa->toTPA = psPlaneTransformInvert(NULL, output->fpa->fromTPA, *fpaRegion, 50); + + psMetadata *header = output->fpa->hdu->header; + pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL); + + psFree (fpaRegion); + + return true; +} Index: trunk/psastro/src/psastroModelAnalysis.c =================================================================== --- trunk/psastro/src/psastroModelAnalysis.c (revision 15890) +++ trunk/psastro/src/psastroModelAnalysis.c (revision 15891) @@ -12,4 +12,7 @@ } + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL"); + if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL"); + // physical pixel scale in microns per pixel char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP"); @@ -19,22 +22,18 @@ } - /* model analysis steps: - - 1) determine POS_ZERO: - * choose a reference chip - * what is posangle of chip - * posangle = FPA.POSANGLE - POS_ZERO - * POS_ZERO = FPA.POSANGLE - posangle - - ** average over input images - - 2) determine boresite model - - * find CRPIX1,2 for a series of images - * fit ellipse - - */ - - /***** find the POSANGLE offset *****/ + /* model analysis: + * + * determine POS_ZERO via comparison of measured and reported posangles + * POS_ZERO = FPA.POSANGLE - posangle + * + * determine boresite model: + * X = Xo + R_X cos(FPA.POSANGLE - T_0) cos(P_0) + R_Y sin(FPA.POSANGLE - T_0) sin(P_0) + * Y = Yo + R_Y sin(FPA.POSANGLE - T_0) cos(P_0) - R_X cos(FPA.POSANGLE - T_0) sin(P_0) + * position of reported boresite in reference chip pixels + * Xo, Yo : true coordinate of boresite (rotator center) in reference chip pixels + * R_X, R_Y : amplitude of boresite offset + * T_0 : reference angle for rotator + * P_0 : orientation of boresite ellipse + */ // select the input pmFPAfile pointers @@ -44,25 +43,39 @@ psArray *files = psListToArray (item->data.list); - // temp data vectors + // data storage vectors for measurements psVector *posZero = psVectorAlloc (files->n, PS_TYPE_F32); psVector *Po = psVectorAlloc (files->n, PS_TYPE_F32); - psVector *Lo = psVectorAlloc (files->n, PS_TYPE_F32); - psVector *Mo = psVectorAlloc (files->n, PS_TYPE_F32); + psVector *Xo = psVectorAlloc (files->n, PS_TYPE_F32); + psVector *Yo = psVectorAlloc (files->n, PS_TYPE_F32); + // counter for accepted measured values int n = 0; - // convert the headers for the input file into fpa astrometry terms + // Re-select the output fpa. The output fpa needs to have valid astrometry for the + // refChip. output->fpa is a copy of the pointer to one of the input->fpa, but the choice + // is arbitrary. select a new one that has an existing ref chip + psFree (output->fpa); + output->fpa = NULL; + + // extract the relevant measured and reported values from the reference chip for (int i = 0; i < files->n; i++) { psMetadataItem *file = files->data[i]; pmFPAfile *input = file->data.V; + // reported rotator position angle (this should perhaps be ROT, not POS) double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE"); if (!status) psAbort ("missing FPA.POSANGLE"); - // get chip from name + // get reference chip from name pmChip *chip = pmConceptsChipFromName (input->fpa, refChip); if (!chip) psAbort ("invalid chip name for reference"); + fprintf (stderr, "input %d : %zx : %zx : %zx\n", i, (size_t) input, (size_t) chip, (size_t) chip->toFPA); if (!chip->toFPA) continue; + + if (!output->fpa) { + // this one matches + output->fpa = psMemIncrRefCounter(input->fpa); + } // we have two measurements of the posangle (may be parity flipped to different quadrants) @@ -72,13 +85,13 @@ while (posZero->data.F32[n] < 0.0) posZero->data.F32[n] += 360.0; - Po->data.F32[n] = POSANGLE * PM_RAD_DEG; - Lo->data.F32[n] = chip->toFPA->x->coeff[0][0]; - Mo->data.F32[n] = chip->toFPA->y->coeff[0][0]; - fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Lo->data.F32[n], Mo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle); + Po->data.F32[n] = POSANGLE * PM_RAD_DEG; // reported position angle + Xo->data.F32[n] = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates + Yo->data.F32[n] = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates + fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle); n ++; } - Lo->n = n; - Mo->n = n; + Xo->n = n; + Yo->n = n; Po->n = n; posZero->n = n; @@ -88,6 +101,16 @@ fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POS_ZERO", PS_META_REPLACE, "offset between obs and meas posangle", stats->sampleMedian); - psastroModelFitBoresite (Lo, Mo, Po); + psVector *params = psastroModelFitBoresite (Xo, Yo, Po); + if (params->n != 6) psAbort ("error"); + + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_X0]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_Y0]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RX]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RY]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_T0]); + psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_P0]); + psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip); return true; Index: trunk/psastro/src/psastroModelBoresite.c =================================================================== --- trunk/psastro/src/psastroModelBoresite.c (revision 15890) +++ trunk/psastro/src/psastroModelBoresite.c (revision 15891) @@ -1,6 +1,8 @@ # include "psastroStandAlone.h" -// chisq = sum ((L_obs - L_fit(t))^2 + (M_obs - M_fit(t))^2) -// coord[0] = measured L or measured M +// the full chisq is built of two associated sums over coordinates: +// chisq = sum ((X_obs - X_fit(t))^2 + (Y_obs - Y_fit(t))^2) +// we use split this into a 2x long vector and use coord[1] to distinguish the X and Y terms: +// coord[0] = measured X or measured Y // coord[1] = 0 or 1 psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord) { @@ -15,33 +17,33 @@ float sntht = sin(dtheta); - // value is L + // value is X if (coord->data.F32[1] == 0) { - float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi + PAR[PAR_RM]*sntht*snphi; + float value = PAR[PAR_X0] + PAR[PAR_RX]*cstht*csphi + PAR[PAR_RY]*sntht*snphi; if (deriv) { psF32 *dPAR = deriv->data.F32; - dPAR[PAR_L0] = 1.0; - dPAR[PAR_M0] = 0.0; - dPAR[PAR_RL] = +cstht*csphi; - dPAR[PAR_RM] = +sntht*snphi; - dPAR[PAR_P0] = -PAR[PAR_RL]*cstht*snphi + PAR[PAR_RM]*sntht*csphi; - dPAR[PAR_T0] = PAR[PAR_RL]*sntht*csphi - PAR[PAR_RM]*cstht*snphi; + dPAR[PAR_X0] = 1.0; + dPAR[PAR_Y0] = 0.0; + dPAR[PAR_RX] = +cstht*csphi; + dPAR[PAR_RY] = +sntht*snphi; + dPAR[PAR_P0] = -PAR[PAR_RX]*cstht*snphi + PAR[PAR_RY]*sntht*csphi; + dPAR[PAR_T0] = PAR[PAR_RX]*sntht*csphi - PAR[PAR_RY]*cstht*snphi; } return (value); } - // value is M + // value is Y if (coord->data.F32[1] == 1) { - float value = PAR[PAR_M0] + PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi; + float value = PAR[PAR_Y0] + PAR[PAR_RY]*sntht*csphi - PAR[PAR_RX]*cstht*snphi; if (deriv) { psF32 *dPAR = deriv->data.F32; - dPAR[PAR_L0] = 0.0; - dPAR[PAR_M0] = 1.0; - dPAR[PAR_RL] = -cstht*snphi; - dPAR[PAR_RM] = +sntht*csphi; - dPAR[PAR_P0] = -PAR[PAR_RM]*sntht*snphi - PAR[PAR_RL]*cstht*csphi; - dPAR[PAR_T0] = -PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi; + dPAR[PAR_X0] = 0.0; + dPAR[PAR_Y0] = 1.0; + dPAR[PAR_RX] = -cstht*snphi; + dPAR[PAR_RY] = +sntht*csphi; + dPAR[PAR_P0] = -PAR[PAR_RY]*sntht*snphi - PAR[PAR_RX]*cstht*csphi; + dPAR[PAR_T0] = -PAR[PAR_RY]*sntht*csphi - PAR[PAR_RX]*cstht*snphi; } return (value); @@ -49,53 +51,2 @@ psAbort ("programming error: invalid coordinate"); } - -# undef PAR_L0 -# undef PAR_RL -# undef PAR_PHI -# define PAR_L0 0 -# define PAR_RL 1 -# define PAR_PHI 2 - -psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord) { - - psF32 *PAR = params->data.F32; - - float csphi = cos(PAR[PAR_PHI]); - float snphi = sin(PAR[PAR_PHI]); - - float theta = coord->data.F32[0]; - float cstht = cos(theta); - float sntht = sin(theta); - - float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi - PAR[PAR_RL]*sntht*snphi; - - if (deriv) { - psF32 *dPAR = deriv->data.F32; - dPAR[PAR_L0] = 1.0; - dPAR[PAR_RL] = cstht*csphi - sntht*snphi; - dPAR[PAR_PHI] = -PAR[PAR_RL]*cstht*snphi - PAR[PAR_RL]*sntht*csphi; - } - return (value); -} - -psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord) { - - psF32 *PAR = params->data.F32; - - float csphi = cos(PAR[PAR_PHI]); - float snphi = sin(PAR[PAR_PHI]); - - float theta = coord->data.F32[0]; - float cstht = cos(theta); - float sntht = sin(theta); - - float value = PAR[PAR_L0] + PAR[PAR_RL]*sntht*csphi + PAR[PAR_RL]*cstht*snphi; - - if (deriv) { - psF32 *dPAR = deriv->data.F32; - dPAR[PAR_L0] = 1.0; - dPAR[PAR_RL] = +sntht*csphi + cstht*snphi; - dPAR[PAR_PHI] = -PAR[PAR_RL]*sntht*snphi + PAR[PAR_RL]*cstht*csphi; - } - return (value); -} Index: trunk/psastro/src/psastroModelDataLoad.c =================================================================== --- trunk/psastro/src/psastroModelDataLoad.c (revision 15890) +++ trunk/psastro/src/psastroModelDataLoad.c (revision 15891) @@ -33,6 +33,6 @@ // select the input data sources - pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL"); - if (!output) psAbort ("PSASTRO.MODEL not listed in config->files"); + pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUT.MODEL"); + if (!output) psAbort ("PSASTRO.OUT.MODEL not listed in config->files"); pmFPAview *view = pmFPAviewAlloc (0); Index: trunk/psastro/src/psastroModelDataSave.c =================================================================== --- trunk/psastro/src/psastroModelDataSave.c (revision 15891) +++ trunk/psastro/src/psastroModelDataSave.c (revision 15891) @@ -0,0 +1,36 @@ +# include "psastroStandAlone.h" + +# define ESCAPE { \ + psError(PS_ERR_UNKNOWN, false, "Failure in psastroModelDataSave"); \ + psFree (view); \ + return false; \ +} + +bool psastroModelDataSave (pmConfig *config) { + + bool status; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL"); + if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL"); + + pmFPAview *view = pmFPAviewAlloc (0); + + pmChip *chip = NULL; + while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; + + psLogMsg ("psastro", 3, "save headers : %f sec\n", psTimerMark ("psastro")); + + return true; +} Index: trunk/psastro/src/psastroModelFitBoresite.c =================================================================== --- trunk/psastro/src/psastroModelFitBoresite.c (revision 15890) +++ trunk/psastro/src/psastroModelFitBoresite.c (revision 15891) @@ -2,33 +2,33 @@ // we now have a set of observed L,M values. fit these to the boresite model -bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) { +psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po) { - assert (Lo->n > 2); - assert (Lo->n == Mo->n); - assert (Lo->n == Po->n); + assert (Xo->n > 2); + assert (Xo->n == Yo->n); + assert (Xo->n == Po->n); // arrays to hold the data to be fitted - psArray *x = psArrayAlloc(2*Lo->n); - psVector *y = psVectorAlloc(2*Lo->n, PS_TYPE_F32); + psArray *x = psArrayAlloc(2*Xo->n); + psVector *y = psVectorAlloc(2*Xo->n, PS_TYPE_F32); int n = 0; - for (int i = 0; i < Lo->n; i++) { + for (int i = 0; i < Xo->n; i++) { psVector *coord = NULL; - // L coordinate value + // X coordinate value coord = psVectorAlloc (2, PS_TYPE_F32); coord->data.F32[1] = 0.0; coord->data.F32[0] = Po->data.F32[i]; x->data[n] = coord; - y->data.F32[n] = Lo->data.F32[i]; + y->data.F32[n] = Xo->data.F32[i]; n++; - // M coordinate value + // Y coordinate value coord = psVectorAlloc (2, PS_TYPE_F32); coord->data.F32[1] = 1.0; coord->data.F32[0] = Po->data.F32[i]; x->data[n] = coord; - y->data.F32[n] = Mo->data.F32[i]; + y->data.F32[n] = Yo->data.F32[i]; n++; } @@ -47,13 +47,13 @@ psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN); - // center (Lo) = mean(Lo), RL = range / 2 - psVectorStats (stats, Lo, NULL, NULL, 0); - params->data.F32[PAR_L0] = stats->sampleMean; - params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0; + // center (Xo) = mean(Xo), RX = range / 2 + psVectorStats (stats, Xo, NULL, NULL, 0); + params->data.F32[PAR_X0] = stats->sampleMean; + params->data.F32[PAR_RX] = (stats->max - stats->min) / 2.0; - // center (Mo) = mean(Mo), RM = range / 2 - psVectorStats (stats, Mo, NULL, NULL, 0); - params->data.F32[PAR_M0] = stats->sampleMean; - params->data.F32[PAR_RM] = (stats->max - stats->min) / 2.0; + // center (Yo) = mean(Yo), RY = range / 2 + psVectorStats (stats, Yo, NULL, NULL, 0); + params->data.F32[PAR_Y0] = stats->sampleMean; + params->data.F32[PAR_RY] = (stats->max - stats->min) / 2.0; params->data.F32[PAR_P0] = 0.0; @@ -63,11 +63,11 @@ psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - fprintf (stderr, "guess values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "Mo: %f\n", params->data.F32[PAR_M0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "RM: %f\n", params->data.F32[PAR_RM]); - fprintf (stderr, "P0: %f\n", params->data.F32[PAR_P0]); - fprintf (stderr, "T0: %f\n", params->data.F32[PAR_T0]); + // fprintf (stderr, "guess values:\n"); + // fprintf (stderr, "Xo: %f\n", params->data.F32[PAR_X0]); + // fprintf (stderr, "Yo: %f\n", params->data.F32[PAR_Y0]); + // fprintf (stderr, "RX: %f\n", params->data.F32[PAR_RX]); + // fprintf (stderr, "RY: %f\n", params->data.F32[PAR_RY]); + // fprintf (stderr, "P0: %f\n", params->data.F32[PAR_P0]); + // fprintf (stderr, "T0: %f\n", params->data.F32[PAR_T0]); // XXX skip the weights for now @@ -75,136 +75,11 @@ fprintf (stderr, "fitted values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "Mo: %f\n", params->data.F32[PAR_M0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "RM: %f\n", params->data.F32[PAR_RM]); + fprintf (stderr, "Xo: %f\n", params->data.F32[PAR_X0]); + fprintf (stderr, "Yo: %f\n", params->data.F32[PAR_Y0]); + fprintf (stderr, "RX: %f\n", params->data.F32[PAR_RX]); + fprintf (stderr, "RY: %f\n", params->data.F32[PAR_RY]); fprintf (stderr, "P0: %f\n", params->data.F32[PAR_P0]); fprintf (stderr, "T0: %f\n", params->data.F32[PAR_T0]); - return true; + return params; } - -# if (0) -# undef PAR_L0 -# undef PAR_RL -# undef PAR_PHI -# define PAR_L0 0 -# define PAR_RL 1 -# define PAR_PHI 2 - -// we now have a set of observed L,M values. fit these to the boresite model -bool psastroModelFitBoresite_L (psVector *Lo, psVector *Mo, psVector *Po) { - - assert (Lo->n > 2); - assert (Lo->n == Mo->n); - assert (Lo->n == Po->n); - - // arrays to hold the data to be fitted - psArray *x = psArrayAlloc(Lo->n); - psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32); - - for (int i = 0; i < Lo->n; i++) { - - psVector *coord = NULL; - - // L coordinate value - coord = psVectorAlloc (1, PS_TYPE_F32); - coord->data.F32[0] = Po->data.F32[i]; - x->data[i] = coord; - y->data.F32[i] = Lo->data.F32[i]; - } - - psVector *params = psVectorAlloc (3, PS_TYPE_F32); - - // create the minimization constraints - psMinConstraint *constraint = psMinConstraintAlloc(); - - // XXX for now, no parameter masks, skip checkLimits - // constraint->checkLimits = psastroModelBoresiteLimits; - - // make an initial guess: - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN); - - // center (Lo) = mean(Lo), RL = range / 2 - psVectorStats (stats, Lo, NULL, NULL, 0); - params->data.F32[PAR_L0] = stats->sampleMean; - params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0; - params->data.F32[PAR_PHI] = 0.0; - - psMinimization *myMin = psMinimizationAlloc (15, 0.001); - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - - fprintf (stderr, "guess values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]); - - // XXX skip the weights for now - psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteL); - - fprintf (stderr, "fitted values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]); - - return true; -} - -// we now have a set of observed L,M values. fit these to the boresite model -bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) { - - assert (Lo->n > 2); - assert (Lo->n == Mo->n); - assert (Lo->n == Po->n); - - // arrays to hold the data to be fitted - psArray *x = psArrayAlloc(Lo->n); - psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32); - - for (int i = 0; i < Lo->n; i++) { - - psVector *coord = NULL; - - // L coordinate value - coord = psVectorAlloc (1, PS_TYPE_F32); - coord->data.F32[0] = Po->data.F32[i]; - x->data[i] = coord; - y->data.F32[i] = Mo->data.F32[i]; - } - - psVector *params = psVectorAlloc (3, PS_TYPE_F32); - - // create the minimization constraints - psMinConstraint *constraint = psMinConstraintAlloc(); - - // XXX for now, no parameter masks, skip checkLimits - // constraint->checkLimits = psastroModelBoresiteLimits; - - // make an initial guess: - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN); - - // center (Lo) = mean(Lo), RL = range / 2 - psVectorStats (stats, Mo, NULL, NULL, 0); - params->data.F32[PAR_L0] = stats->sampleMean; - params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0; - params->data.F32[PAR_PHI] = 0.0; - - psMinimization *myMin = psMinimizationAlloc (15, 0.001); - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); - - fprintf (stderr, "guess values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]); - - // XXX skip the weights for now - psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteM); - - fprintf (stderr, "fitted values:\n"); - fprintf (stderr, "Lo: %f\n", params->data.F32[PAR_L0]); - fprintf (stderr, "RL: %f\n", params->data.F32[PAR_RL]); - fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]); - - return true; -} -# endif Index: trunk/psastro/src/psastroModelParseCamera.c =================================================================== --- trunk/psastro/src/psastroModelParseCamera.c (revision 15890) +++ trunk/psastro/src/psastroModelParseCamera.c (revision 15891) @@ -18,7 +18,7 @@ // set up an output fpa structure & file based on the last input file - pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.MODEL"); + pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUT.MODEL"); if (!output) { - psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA for PSASTRO.MODEL from input"); + psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA for PSASTRO.OUT.MODEL from input"); return false; } Index: trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- trunk/psastro/src/psastroMosaicAstrom.c (revision 15890) +++ trunk/psastro/src/psastroMosaicAstrom.c (revision 15891) @@ -23,5 +23,5 @@ // before we do object matches, we need to (optionally) fix failed chips. We compare chips with // the supplied mosaic model. Adjust significant outliers to match model. - if (!psastroEnforceChips (config, fpa, recipe)) { + if (!psastroFixChips (config, recipe)) { psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips"); return false; Index: trunk/psastro/src/psastroStandAlone.h =================================================================== --- trunk/psastro/src/psastroStandAlone.h (revision 15890) +++ trunk/psastro/src/psastroStandAlone.h (revision 15891) @@ -17,20 +17,19 @@ bool psastroDataLoad (pmConfig *config); -pmConfig *psastroModelArguments (int argc, char **argv); -bool psastroModelParseCamera (pmConfig *config); -bool psastroModelDataLoad (pmConfig *config); -bool psastroModelAnalysis (pmConfig *config); +pmConfig *psastroModelArguments (int argc, char **argv); +bool psastroModelParseCamera (pmConfig *config); +bool psastroModelDataLoad (pmConfig *config); +bool psastroModelAnalysis (pmConfig *config); +bool psastroModelAdjust (pmConfig *config); +bool psastroModelDataSave (pmConfig *config); -bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po); -psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord); - -psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord); -psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord); +psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po); +psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord); // these are used to define the boresite model parameters -# define PAR_L0 0 // Lo = params[0] -# define PAR_M0 1 // Mo = params[1] -# define PAR_RL 2 // RL = params[2] -# define PAR_RM 3 // RM = params[3] +# define PAR_X0 0 // Xo = params[0] +# define PAR_Y0 1 // Yo = params[1] +# define PAR_RX 2 // RX = params[2] +# define PAR_RY 3 // RY = params[3] # define PAR_P0 4 // P0 = params[4] # define PAR_T0 5 // phi = params[4] Index: trunk/psastro/src/psastroUseModel.c =================================================================== --- trunk/psastro/src/psastroUseModel.c (revision 15891) +++ trunk/psastro/src/psastroUseModel.c (revision 15891) @@ -0,0 +1,78 @@ +# include "psastroInternal.h" +# define NONLIN_TOL 0.001 /* tolerance in pixels */ + +// apply the generic astrometry model to this image. this assumes the WCS terms either do not +// exist or are invalid +bool psastroUseModel (pmConfig *config, psMetadata *recipe) { + + bool status; + + bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL"); + if (!status) { + useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL"); + } + if (!useModel) return true; + + // identify reference astrometry table. + // if not defined, correction was not requested; skip step + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL"); + if (!astrom) psAbort ("programming error: model was not supplied, though requested"); + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) psAbort ("programming error: no input data"); + + // make sure the astrometry model is loaded. de-activate all files except PSASTRO.MODEL. + pmFPAfileActivate (config->files, false, NULL); + pmFPAfileActivate (config->files, true, "PSASTRO.MODEL"); + + pmFPAview *view = pmFPAviewAlloc (0); + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError (PS_ERR_IO, false, "Can't load the astrometry model file"); + return false; + } + + // set the model using the RA, DEC, POSANGLE of the input image. + pmAstromSetTP (astrom, input->fpa->concepts); + + // loop over all chips, replace input astrometry elements with those from astrom + pmChip *obsChip = NULL; + while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process); + if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; } + + // set the chip astrometry using the astrom file + pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa); + + psFree (obsChip->toFPA); + psFree (obsChip->fromFPA); + + // supply astrometry from model + obsChip->toFPA = psMemIncrRefCounter (refChip->toFPA); + obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA); + + // XXX if we want to write out the result, update the header here. this needs to be + // updated with the correct HDU selection. obsChip->hdu may not exist. + // pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL); + } + + psFree (input->fpa->toSky); + psFree (input->fpa->toTPA); + psFree (input->fpa->fromTPA); + input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky); + input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA); + input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA); + + // XXX if we want to write out the result, update the header here. this needs to be + // updated with the correct HDU selection. obsChip->hdu may not exist. + // psMetadata *updates = psMetadataAlloc(); + // pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL); + // psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); + // psFree (updates); + + psFree (view); + return true; +} +