Index: /trunk/psastro/src/psastro.c =================================================================== --- /trunk/psastro/src/psastro.c (revision 10854) +++ /trunk/psastro/src/psastro.c (revision 10855) @@ -53,11 +53,14 @@ } - if (!psastroChipAstrom (config, refs)) { - psErrorStackPrint(stderr, "failed to perform single chip astrometry\n"); - exit (1); + bool chipastro = psMetadataLookupBool (NULL, config->arguments, "PSASTRO.CHIP.MODE"); + if (chipastro) { + if (!psastroChipAstrom (config, refs)) { + psErrorStackPrint(stderr, "failed to perform single chip astrometry\n"); + exit (1); + } } - bool mosastro = psMetadataLookupBool (NULL, config->arguments, "MOSAIC.MODE"); - if (!mosastro) { + bool mosastro = psMetadataLookupBool (NULL, config->arguments, "PSASTRO.MOSAIC.MODE"); + if (mosastro) { if (!psastroMosaicAstrom (config, refs)) { psErrorStackPrint(stderr, "failed to perform mosaic camera astrometry\n"); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 10854) +++ /trunk/psastro/src/psastroArguments.c (revision 10855) @@ -33,13 +33,26 @@ if ((N = psArgumentGet (argc, argv, "-chip"))) { psArgumentRemove (N, &argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", psStringCopy(argv[N])); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_META_REPLACE, "", psStringCopy(argv[N])); psArgumentRemove (N, &argc, argv); } - // mosastro mode also specifies output header file + // apply mosastro mode? if ((N = psArgumentGet (argc, argv, "-mosastro"))) { psArgumentRemove (N, &argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MOSAIC.MODE", PS_DATA_BOOL, "", psStringCopy(argv[N])); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.MOSAIC.MODE", PS_META_REPLACE, "", false); + } + if ((N = psArgumentGet (argc, argv, "+mosastro"))) { psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.MOSAIC.MODE", PS_META_REPLACE, "", true); + } + + // apply chipastro mode? + if ((N = psArgumentGet (argc, argv, "-chipastro"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.CHIP.MODE", PS_META_REPLACE, "", false); + } + if ((N = psArgumentGet (argc, argv, "+chipastro"))) { + psArgumentRemove (N, &argc, argv); + psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.CHIP.MODE", PS_META_REPLACE, "", true); } Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 10854) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 10855) @@ -11,17 +11,15 @@ bool newFPA = true; bool status = false; - double RAmin = NAN; - double RAmax = NAN; + double RAmin = FLT_MAX; + double RAmax = FLT_MIN; + double DECmin = FLT_MAX; + double DECmax = FLT_MIN; + double RAminSky = NAN; double RAmaxSky = NAN; - double DECmin = NAN; - double DECmax = NAN; pmChip *chip = NULL; pmCell *cell = NULL; pmReadout *readout = NULL; - - DECmin = DECmax = -90; - RAmin = RAmax = RAminSky = RAmaxSky = 0; // select the current recipe @@ -63,6 +61,4 @@ RAminSky = fpa->toSky->R - M_PI; RAmaxSky = fpa->toSky->R + M_PI; - RAmin = RAmax = fpa->toSky->R; - DECmin = DECmax = fpa->toSky->D; } @@ -79,4 +75,5 @@ if (rawstars == NULL) { continue; } + FILE *f = fopen ("rawstars.dat", "w"); for (int i = 0; i < rawstars->n; i++) { pmAstromObj *raw = rawstars->data[i]; @@ -85,4 +82,10 @@ psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); psDeproject (raw->sky, raw->TP, fpa->toSky); + + fprintf (f, "%d %f %f %f %f %f %f %f %f\n", i, + raw->sky->r, raw->sky->d, + raw->TP->x, raw->TP->y, + raw->FP->x, raw->FP->y, + raw->chip->x, raw->chip->y); if (i < 10) { @@ -122,4 +125,5 @@ DECmax = PS_MAX (raw->sky->d, DECmax); } + fclose (f); } } Index: /trunk/psastro/src/psastroChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroChipAstrom.c (revision 10854) +++ /trunk/psastro/src/psastroChipAstrom.c (revision 10855) @@ -49,10 +49,19 @@ if (refstars == NULL) { continue; } + // the absolute minimum number of stars is 4 (for order = 1) + if ((rawstars->n < 4) || (refstars->n < 4)) { + psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", + rawstars->n, refstars->n); + continue; + } + // save WCS and analysis metadata in update header psMetadata *updates = psMetadataAlloc(); + // XXX update the header with info to reflect the failure if (!psastroOneChip (fpa, chip, refstars, rawstars, recipe, updates)) { - psError (PSASTRO_ERR_DATA, false, "failed to find a solution\n"); - return false; + psLogMsg ("psastro", 3, "failed to find a solution\n"); + psFree (updates); + continue; } pmAstromWriteWCS (updates, fpa, chip, NONLIN_TOL); Index: /trunk/psastro/src/psastroChooseRefstars.c =================================================================== --- /trunk/psastro/src/psastroChooseRefstars.c (revision 10854) +++ /trunk/psastro/src/psastroChooseRefstars.c (revision 10855) @@ -56,4 +56,6 @@ psArray *refstars = psArrayAllocEmpty (100); + FILE *f = fopen ("refstars.dat", "w"); + // select the reference objects within range of this readout // project the reference objects to this chip @@ -64,4 +66,10 @@ psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP); psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + + fprintf (f, "%d %f %f %f %f %f %f %f %f\n", i, + ref->sky->r, ref->sky->d, + ref->TP->x, ref->TP->y, + ref->FP->x, ref->FP->y, + ref->chip->x, ref->chip->y); // limit the X,Y range of the refs to the selected chip @@ -75,4 +83,5 @@ psFree (ref); } + fclose (f); psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n); Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 10854) +++ /trunk/psastro/src/psastroConvert.c (revision 10855) @@ -39,4 +39,6 @@ psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars); + psLogMsg ("psastro", 4, "loaded %ld sources\n", rawstars->n); + psFree (rawstars); Index: /trunk/psastro/src/psastroLuminosityFunction.c =================================================================== --- /trunk/psastro/src/psastroLuminosityFunction.c (revision 10854) +++ /trunk/psastro/src/psastroLuminosityFunction.c (revision 10855) @@ -109,5 +109,5 @@ } - psLogMsg ("psastro", 4, "logN %f + %f mag\n", line->coeff[0], line->coeff[1]); + psLogMsg ("psastro", 4, "logN = %f + %f mag\n", line->coeff[0], line->coeff[1]); psLogMsg ("psastro", 4, "logN vs mag residuals: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); Index: /trunk/psastro/src/psastroOneChip.c =================================================================== --- /trunk/psastro/src/psastroOneChip.c (revision 10854) +++ /trunk/psastro/src/psastroOneChip.c (revision 10855) @@ -1,3 +1,9 @@ # include "psastro.h" + +# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\ + VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \ + if (!status) { \ + psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \ + return false; } bool psastroOneChip (pmFPA *fpa, pmChip *chip, psArray *refstars, psArray *rawstars, psMetadata *recipe, psMetadata *updates) { @@ -5,8 +11,23 @@ bool status; + // supplied radius is in pixels + REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32, "Failed to lookup matching radius"); + + // correct radius to FP units (physical pixel scale in microns per pixel) + REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); + RADIUS *= pixelScale; + + // select the desired chip order + REQUIRED_RECIPE_VALUE (int order, "PSASTRO.CHIP.ORDER", S32, "failed to find single-chip fit order\n"); + + // allowed limits for valid solutions + REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32, "failed to find single-chip max allowed error\n"); + + REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n"); + // find initial offset / rotation pmAstromStats *gridStats = pmAstromGridMatch (rawstars, refstars, recipe); if (gridStats == NULL) { - psError(PSASTRO_ERR_DATA, false, "failed to find a grid match solution\n"); + psLogMsg ("psastro", 3, "failed to find a grid match solution\n"); return false; } @@ -16,5 +37,5 @@ pmAstromStats *stats = pmAstromGridTweak (rawstars, refstars, recipe, gridStats); if (stats == NULL) { - psError(PSASTRO_ERR_DATA, false, "failed to measure tweaked grid solution\n"); + psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n"); return false; } @@ -25,32 +46,25 @@ psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); - // supplied radius is in pixels - double RADIUS = psMetadataLookupF32 (&status, recipe, "PSASTRO.MATCH.RADIUS"); - if (!status) { - psError(PS_ERR_IO, false, "Failed to lookup matching radius"); + // use small radius to match stars + psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS); + if (match == NULL) { + psLogMsg ("psastro", 3, "failed to find radius-matched sources\n"); + return false; + } + + // create the output fit model, modify the order to correspond to the actual number of + // matched stars: + if ((match->n < 11) && (order >= 3)) order = 2; + if ((match->n < 7) && (order >= 2)) order = 1; + if ((match->n < 4) && (order >= 1)) order = 0; + if (order < 1) { + psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); + psFree (match); + psFree (stats); + psFree (gridStats); return false; } - // correct to FP units (physical pixel scale in microns per pixel) - double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); - if (!status) { - psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); - return false; - } - RADIUS *= pixelScale; - - // use small radius to match stars - psArray *match = pmAstromRadiusMatchFP (rawstars, refstars, RADIUS); - if (stats == NULL) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to find radius-matched sources\n"); - return false; - } - - // create the output fit model - int order = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.ORDER"); - if (!status) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n"); - return false; - } + // set masks appropriate to the Elixir DVO astrometry format psFree (chip->toFPA); chip->toFPA = psPlaneTransformAlloc (order, order); @@ -72,14 +86,18 @@ pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); if (!results) { - psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); + psLogMsg ("psastro", 3, "failed to perform the matched fit\n"); return false; } - // allowed limits for valid solutions - float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.ERROR"); - int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MIN.NSTAR"); + // toTPA converts from FPA units (microns) to TPA units (linear degrees) + float plateScale1 = hypot (fpa->toTPA->x->coeff[1][0], fpa->toTPA->x->coeff[0][1]); + float plateScale2 = hypot (fpa->toTPA->y->coeff[1][0], fpa->toTPA->y->coeff[0][1]); + float plateScale = 0.5*(plateScale1 + plateScale2)*3600.0; - // astError is the average 1D scatter in pixels ('results' are in FPA units = microns) - float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; + // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns) + float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; + + // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns) + float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; int astNstar = results->yStats->clippedNvalues; @@ -98,13 +116,14 @@ // DVO expects NASTRO = 0 if we fail to find a solution + psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR", PS_META_REPLACE, "astrometry error (pixels)", pixError); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "astrometry error (arcsec)", astError); if (validSolution) { - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", astNstar); - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar)); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar)); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", astNstar); } else { - psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "", 0); - psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0); + psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0); } - psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR", PS_META_REPLACE, "", astError); - psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars + psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars // write results Index: /trunk/psastro/src/psastroParseCamera.c =================================================================== --- /trunk/psastro/src/psastroParseCamera.c (revision 10854) +++ /trunk/psastro/src/psastroParseCamera.c (revision 10855) @@ -22,6 +22,8 @@ pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUTPUT"); - // Chip selection: turn on only the chips specified + // Chip selection: turn on only the chips specified (option is not required) char *chipLine = psMetadataLookupStr(NULL, config->arguments, "CHIP_SELECTIONS"); + if (chipLine == NULL) psErrorClear(); + psArray *chips = psStringSplitArray (chipLine, ",", false); if (chips->n > 0) { Index: /trunk/psastro/src/psastroRefstarSubset.c =================================================================== --- /trunk/psastro/src/psastroRefstarSubset.c (revision 10854) +++ /trunk/psastro/src/psastroRefstarSubset.c (revision 10855) @@ -61,2 +61,10 @@ return true; } + +/* this test is a bit sensitive to the total number of refstars or rawstars available + watch out if: + - the fitted slopes are extremely different + - the average number of stars per bin is ~1 + + skip the cut in these cases? +*/ Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 10854) +++ /trunk/psastro/src/psastroUtils.c (revision 10855) @@ -37,4 +37,5 @@ for (int i = 0; i < fpa->chips->n; i++) { pmChip *chip = fpa->chips->data[i]; + if (!chip->process || !chip->file_exists) { continue; } pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]); @@ -50,4 +51,5 @@ for (int i = 0; i < fpa->chips->n; i++) { pmChip *chip = fpa->chips->data[i]; + if (!chip->process || !chip->file_exists) { continue; } psPlaneTransform *toFPA = chip->toFPA;