Index: trunk/psastro/src/Makefile.am =================================================================== --- trunk/psastro/src/Makefile.am (revision 15538) +++ trunk/psastro/src/Makefile.am (revision 15562) @@ -33,4 +33,5 @@ psastroLuminosityFunction.c \ psastroRefstarSubset.c \ + psastroEnforceChips.c \ psastroMosaicAstrom.c \ psastroMosaicGradients.c \ Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 15538) +++ trunk/psastro/src/psastro.h (revision 15562) @@ -27,4 +27,5 @@ bool psastroDataSave (pmConfig *config); bool psastroDefineFiles (pmConfig *config, pmFPAfile *input); +bool psastroDefineFile (pmConfig *config, pmFPA *input, char *filerule, char *argname, pmFPAfileType fileType, pmDetrendType detrendType); bool psastroAnalysis (pmConfig *config); @@ -74,7 +75,13 @@ bool psastroDumpStars (psArray *stars, char *filename); bool psastroDumpMatchedStars (char *filename, psArray *rawstars, psArray *refstars, psArray *match); +bool psastroDumpGradients (psArray *gradients, char *filename); bool psastroMosaicSetAstrom_tmp (pmFPA *fpa); int psastroSortByMag (const void *a, const void *b); +bool psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe); + +psArray *psastroReadGetstarCatalog (psFits *fits); +psArray *psastroReadGetstar_PS1_DEV_0 (psFits *fits); + # endif /* PSASTRO_H */ Index: trunk/psastro/src/psastroAnalysis.c =================================================================== --- trunk/psastro/src/psastroAnalysis.c (revision 15538) +++ trunk/psastro/src/psastroAnalysis.c (revision 15562) @@ -16,7 +16,4 @@ return true; } - - // XXX test retrun - // return true; // load the reference stars overlapping the data stars Index: trunk/psastro/src/psastroArguments.c =================================================================== --- trunk/psastro/src/psastroArguments.c (revision 15538) +++ trunk/psastro/src/psastroArguments.c (revision 15562) @@ -36,4 +36,16 @@ psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, &argc, argv); + } + + // apply the chip correction based on the reference astrometry? + if ((N = psArgumentGet (argc, argv, "-fixchips"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); + } + // define the reference astrometry file + status = pmConfigFileSetsMD (config->arguments, &argc, argv, "REF.ASTROM", "-refastrom", "-refastromlist"); + if (status) { + // if supplied, assume -fixchips is desired + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); } Index: trunk/psastro/src/psastroDefineFiles.c =================================================================== --- trunk/psastro/src/psastroDefineFiles.c (revision 15538) +++ trunk/psastro/src/psastroDefineFiles.c (revision 15562) @@ -4,7 +4,14 @@ // these calls bind the I/O handle to the specified fpa - pmFPAfile *output = NULL; + bool status; - output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT"); + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, false, "Can't find PSASTRO recipe!\n"); + return false; + } + + pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT"); if (!output) { psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA from PSASTRO.INPUT"); @@ -12,4 +19,13 @@ } 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"); + return NULL; + } + } # if (0) @@ -25,2 +41,52 @@ return true; } + +bool psastroDefineFile (pmConfig *config, pmFPA *input, char *filerule, char *argname, pmFPAfileType fileType, pmDetrendType detrendType) { + + bool status; + pmFPAfile *file; + + // look for the file on the argument list + file = pmFPAfileDefineFromArgs (&status, config, filerule, argname); + if (!status) { + psError (PS_ERR_UNKNOWN, false, "failed to load find definition"); + return false; + } + if (file) { + if (file->type != fileType) { + psError(PS_ERR_IO, true, "%s is not of type %s", filerule, pmFPAfileStringFromType (fileType)); + return false; + } + return true; + } + + // look for the file in the camera config table + file = pmFPAfileDefineFromConf (&status, config, filerule); + if (!status) { + psError (PS_ERR_UNKNOWN, false, "failed to load find definition"); + return false; + } + if (file) { + if (file->type != fileType) { + psError(PS_ERR_IO, true, "%s is not of type %s", filerule, pmFPAfileStringFromType (fileType)); + return false; + } + return true; + } + + // look for the file to be loaded from the detrend database + file = pmFPAfileDefineFromDetDB (&status, config, filerule, input, detrendType); + if (!status) { + psError (PS_ERR_UNKNOWN, false, "failed to load file definition"); + return false; + } + if (file) { + if (file->type != fileType) { + psError(PS_ERR_IO, true, "%s is not of type %s", filerule, pmFPAfileStringFromType (fileType)); + return false; + } + return true; + } + return false; +} + Index: trunk/psastro/src/psastroDemoDump.c =================================================================== --- trunk/psastro/src/psastroDemoDump.c (revision 15538) +++ trunk/psastro/src/psastroDemoDump.c (revision 15562) @@ -173,2 +173,22 @@ return true; } + +// this function is used for test purposes (-trace psastro.dump 1) +bool psastroDumpGradients (psArray *gradients, char *filename) { + + FILE *f = fopen (filename, "w"); + + for (int i = 0; i < gradients->n; i++) { + pmAstromGradient *gradient = gradients->data[i]; + + // write out the refstar data + fprintf (f, "%d %f %f %f %f %f %f\n", i, + gradient->FP.x, gradient->FP.y, + gradient->dTPdL.x, gradient->dTPdL.y, + gradient->dTPdM.x, gradient->dTPdM.y); + } + + fclose (f); + return true; +} + Index: trunk/psastro/src/psastroEnforceChip.c =================================================================== --- trunk/psastro/src/psastroEnforceChip.c (revision 15538) +++ (revision ) @@ -1,75 +1,0 @@ -# include "psastroInternal.h" - -bool psastroEnforceChips (pmFPA *fpa, psMetadata *recipe) { - - pmChip *chip = NULL; - // identify reference astrometry table, stored as an fpa? better format? - - pmFPAfile *refAstrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.REF.ASTROM"); - - pmFPAview *view = pmFPAviewAlloc (0); - - // physical pixel scale in microns per pixel - double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE"); - if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); - return false; - } - double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE"); - if (!status) { - psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); - return false; - } - - // loop over all chips - while ((obsChip = pmFPAviewNextChip (view, 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 refAstrom file - pmChip *refChip = pmFPAviewThisChip (view, refAstrom->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); - - 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); - - // 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); - } - return true; -} - -// 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 -// 3) allow for this operation to be optional Index: trunk/psastro/src/psastroEnforceChips.c =================================================================== --- trunk/psastro/src/psastroEnforceChips.c (revision 15562) +++ trunk/psastro/src/psastroEnforceChips.c (revision 15562) @@ -0,0 +1,80 @@ +# include "psastroInternal.h" + +bool psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe) { + + bool status; + + // 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; + + pmFPAview *view = pmFPAviewAlloc (0); + + // physical pixel scale in microns per pixel + double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE"); + if (!status) { + psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); + return false; + } + double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE"); + if (!status) { + psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); + return false; + } + + pmChip *obsChip = NULL; + + // loop over all chips + while ((obsChip = pmFPAviewNextChip (view, 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 refAstrom file + pmChip *refChip = pmFPAviewThisChip (view, refAstrom->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); + + // 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); + } + return true; +} + +// 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 +// 3) allow for this operation to be optional Index: trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- trunk/psastro/src/psastroLoadRefstars.c (revision 15538) +++ trunk/psastro/src/psastroLoadRefstars.c (revision 15562) @@ -15,5 +15,4 @@ float DECmin = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MIN"); float DECmax = DEG_RAD*psMetadataLookupF32(NULL, recipe, "DEC_MAX"); - // float MAGmax = psMetadataLookupF32(NULL, recipe, "MAG_MAX"); // XXX CATDIR needs to look up abstracted name from psastro.config @@ -25,4 +24,12 @@ char *getstarCommand = psStringCopy(psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR")); PS_ASSERT (getstarCommand, NULL); + + char *outformat = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.OUTFORMAT"); + PS_ASSERT (outformat, NULL); + + char *photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE"); + PS_ASSERT (photcode, NULL); + + float MAGmax = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MAG.MAX"); // issue the following command: @@ -39,9 +46,9 @@ psTimerStart ("psastro"); - // use fork to add timeout capability + // supply a known output format (for CATALOG output) so the code below knows what to read if (ELIXIR_MODE) { - psStringAppend (&getstarCommand, " -D CATFORMAT elixir"); + psStringAppend (&getstarCommand, " -D CATMODE mef -D CATFORMAT elixir"); } else { - psStringAppend (&getstarCommand, " -D CATFORMAT panstarrs"); + psStringAppend (&getstarCommand, " -D CATMODE mef -D CATFORMAT panstarrs"); } @@ -52,11 +59,13 @@ psFree(CATDIR); - // psStringAppend (&getstarCommand, " -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); - // XXX TEST : no magnitude limit - // XXX need to specify the desired photometry system - psStringAppend (&getstarCommand, " -D CATMODE mef -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile); + // supply the max magnitude, the output format, and the photcode + psStringAppend (&getstarCommand, " -maglim %f -format %s -photcode %s", MAGmax, outformat, photcode); + + // add region and output filename + psStringAppend (&getstarCommand, " -region %f %f %f %f -o %s", RAmin, DECmin, RAmax, DECmax, tempFile); psTrace ("psastro", 3, "%s\n", getstarCommand); // XXX use psPipe: catch stderr, stdout, allow for Nsec timeout... + // use fork to add timeout capability status = system (getstarCommand); if (status) { @@ -71,4 +80,41 @@ psFits *fits = psFitsOpen (tempFile, "r"); + psTimerStart ("psastro"); + + psArray *refstars = NULL; + if (!strcmp (outformat, "CATALOG")) { + refstars = psastroReadGetstarCatalog (fits); + } + if (!strcmp (outformat, "PS1_DEV_0")) { + refstars = psastroReadGetstar_PS1_DEV_0 (fits); + } + if (refstars == NULL) { + psError(PSASTRO_ERR_REFSTARS, true, "error reading reference data\n"); + return NULL; + } + psLogMsg ("psastro", 3, "loaded %ld reference stars : %f sec\n", refstars->n, psTimerMark ("psastro")); + + psTrace ("psastro", 3, "loaded %ld reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n", + refstars->n, RAmin, DECmin, RAmax, DECmax); + + psFitsClose (fits); + unlink (tempFile); + + // dump or plot the available refstars + if (psTraceGetLevel("psastro.dump") > 0) { + psastroDumpRefstars (refstars, "refstars.dat"); + } + + if (psTraceGetLevel("psastro.plot") > 0) { + psastroPlotRefstars (refstars, recipe); + } + + return refstars; +} + +psArray *psastroReadGetstarCatalog (psFits *fits) { + + bool status; + if (ELIXIR_MODE) { psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR"); @@ -77,5 +123,4 @@ } - psTimerStart ("psastro"); long numSources = psFitsTableSize(fits); // Number of sources in table @@ -102,21 +147,30 @@ psFree (row); } - psLogMsg ("psastro", 3, "loaded %ld reference stars : %f sec\n", refstars->n, psTimerMark ("psastro")); - - psTrace ("psastro", 3, "loaded %ld reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n", - refstars->n, RAmin, DECmin, RAmax, DECmax); - - psFitsClose (fits); - unlink (tempFile); - - // dump or plot the available refstars - if (psTraceGetLevel("psastro.dump") > 0) { - psastroDumpRefstars (refstars, "refstars.dat"); - } - - if (psTraceGetLevel("psastro.plot") > 0) { - psastroPlotRefstars (refstars, recipe); - } - return refstars; } + +psArray *psastroReadGetstar_PS1_DEV_0 (psFits *fits) { + + bool status; + + psFitsMoveExtName (fits, "GETSTAR_PS1_DEV_0"); + + long numSources = psFitsTableSize(fits); // Number of sources in table + + // convert the Average table to the pmAstromObj entries + psArray *refstars = psArrayAllocEmpty (numSources); + for (int i = 0; i < numSources; i++) { + pmAstromObj *ref = pmAstromObjAlloc (); + + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + + ref->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); + ref->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); + ref->Mag = psMetadataLookupS32 (&status, row, "MAG"); + + psArrayAdd (refstars, 100, ref); + psFree (ref); + psFree (row); + } + return refstars; +} Index: trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- trunk/psastro/src/psastroMosaicAstrom.c (revision 15538) +++ trunk/psastro/src/psastroMosaicAstrom.c (revision 15562) @@ -21,7 +21,10 @@ pmFPA *fpa = input->fpa; - // XXX before we do object matches, we need to fix failed chips - // compare chips with supplied mosaic model - // adjust significant outliers to match model + // 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)) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips"); + return false; + } // given the existing per-chip astrometry, determine matches between raw and ref stars Index: trunk/psastro/src/psastroMosaicGradients.c =================================================================== --- trunk/psastro/src/psastroMosaicGradients.c (revision 15538) +++ trunk/psastro/src/psastroMosaicGradients.c (revision 15562) @@ -1,3 +1,4 @@ # include "psastroInternal.h" +static int nPass = 0; bool psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) { @@ -50,4 +51,12 @@ } + // if desired, dump the gradients to a file + if (psTraceGetLevel("psastro.dump") > 0) { + char name[80]; + sprintf (name, "gradients.%d.dat", nPass); + psastroDumpGradients (gradients, name); + nPass ++; + } + // allocate mosaic-level polynomial transformation and set masks needed by DVO int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER"); Index: trunk/psastro/src/psastroMosaicOneChip.c =================================================================== --- trunk/psastro/src/psastroMosaicOneChip.c (revision 15538) +++ trunk/psastro/src/psastroMosaicOneChip.c (revision 15562) @@ -55,4 +55,6 @@ return false; } + + psLogMsg ("psastro", PS_LOG_DETAIL, "mosaic fit chip order %d", order); // create output toFPA; set masks appropriate to the Elixir DVO astrometry format if we are Index: trunk/psastro/src/psastroMosaicSetMatch.c =================================================================== --- trunk/psastro/src/psastroMosaicSetMatch.c (revision 15538) +++ trunk/psastro/src/psastroMosaicSetMatch.c (revision 15562) @@ -15,6 +15,15 @@ if (!status) { psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord); - return NULL; + return false; } + + if (RADIUS <= 0.0) { + if (iteration == 0) { + psError(PS_ERR_IO, false, "Invalid match radius for first iteration: %s", radiusWord); + return false; + } + psWarning ("skipping match for iteration %d\n", iteration); + return true; + } // this loop selects the matched stars for all chips