Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 38040) +++ trunk/psastro/src/psastro.h (revision 39926) @@ -70,5 +70,5 @@ bool psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars); -psArray *pmSourceToAstromObj (psArray *sources); +psArray *pmSourceToAstromObj (psArray *sources, float MagOffset); bool psastroAstromGuess (int *nStars, pmConfig *config); bool psastroAstromGuessCheck (pmConfig *config); Index: trunk/psastro/src/psastroAstromGuess.c =================================================================== --- trunk/psastro/src/psastroAstromGuess.c (revision 38040) +++ trunk/psastro/src/psastroAstromGuess.c (revision 39926) @@ -75,5 +75,6 @@ pmFPA *fpa = input->fpa; - if (DEBUG) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa); + // this call only works if we have loaded a model : seg fault if not + if (DEBUG && useModel) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa); // load mosaic-level astrometry? @@ -91,4 +92,10 @@ if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue; } + + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } if (newFPA) { @@ -312,4 +319,10 @@ while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } + + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } // XXX we are currently inconsistent with marking the good vs the bad data Index: trunk/psastro/src/psastroChipAstrom.c =================================================================== --- trunk/psastro/src/psastroChipAstrom.c (revision 38040) +++ trunk/psastro/src/psastroChipAstrom.c (revision 39926) @@ -85,9 +85,9 @@ for (int nn = 0; nn < refstars->n; nn++) { pmAstromObj *ref = refstars->data[nn]; - fprintf (outfile, "%lf %lf %lf %lf %lf %lf %lf %lf\n", + fprintf (outfile, "%lf %lf %lf %lf %lf %lf %lf %lf : %f %f\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y, - ref->chip->x, ref->chip->y); + ref->chip->x, ref->chip->y, ref->Mag, ref->magCal); } fclose (outfile); @@ -101,9 +101,9 @@ for (int nn = 0; nn < gridrawstars->n; nn++) { pmAstromObj *ref = gridrawstars->data[nn]; - fprintf (outfile, "%lf %lf %lf %lf %lf %lf %lf %lf\n", + fprintf (outfile, "%lf %lf %lf %lf %lf %lf %lf %lf : %f %f\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y, - ref->chip->x, ref->chip->y); + ref->chip->x, ref->chip->y, ref->Mag, ref->magCal); } fclose (outfile); Index: trunk/psastro/src/psastroConvert.c =================================================================== --- trunk/psastro/src/psastroConvert.c (revision 38040) +++ trunk/psastro/src/psastroConvert.c (revision 39926) @@ -51,4 +51,29 @@ bool status; + // XXX I want to make this optional + float MagOffset = 0.0; + if (1) { + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "failed to find PSASTRO.INPUT\n"); + return false; + } + pmFPA *fpa = input->fpa; + + float zeropt, exptime; + + // really error-out here? or just skip? + if (!psastroZeroPointFromRecipe (&zeropt, &exptime, NULL, fpa, recipe)) { + psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + zeropt = 0.0; + exptime = 1.0; + } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + MagOffset = zeropt + 2.5*log10(exptime); + } + // PSPHOT.SOURCES carries the pmSource objects (from psphot analysis or loaded externally) pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -59,5 +84,5 @@ // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?) - psArray *inStars = pmSourceToAstromObj (sources); + psArray *inStars = pmSourceToAstromObj (sources, MagOffset); // apply Koppenhoefer correction if needed @@ -175,5 +200,5 @@ // select a magnitude range? -psArray *pmSourceToAstromObj (psArray *sources) { +psArray *pmSourceToAstromObj (psArray *sources, float MagOffset) { psArray *objects = psArrayAllocEmpty (sources->n); @@ -211,4 +236,5 @@ obj->dMag = source->psfMagErr; obj->SBinst = source->psfMag + 5.0*log10(axes.major); + obj->magCal = obj->Mag + MagOffset; // XXX do we have the information giving the readout and cell offset? Index: trunk/psastro/src/psastroDemoDump.c =================================================================== --- trunk/psastro/src/psastroDemoDump.c (revision 38040) +++ trunk/psastro/src/psastroDemoDump.c (revision 39926) @@ -234,4 +234,10 @@ if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; } + if (!chip->toFPA || !chip->fromFPA) { + char *name = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + fprintf (stderr, "no astrom model for %s, skipping\n", name); + continue; + } + // XXX write out the four corners for a test psRegion *region = pmChipPixels (chip); Index: trunk/psastro/src/psastroFixChips.c =================================================================== --- trunk/psastro/src/psastroFixChips.c (revision 38040) +++ trunk/psastro/src/psastroFixChips.c (revision 39926) @@ -306,4 +306,5 @@ psRegion *region = pmChipPixels (obsChip); obsChip->toFPA = toFPA; + // NOTE: when we call psPlaneTransformInvert here, we do not increase the order as in other places obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50); psFree (region); Index: trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- trunk/psastro/src/psastroLoadRefstars.c (revision 38040) +++ trunk/psastro/src/psastroLoadRefstars.c (revision 39926) @@ -265,4 +265,5 @@ ref->Color = 0.0; } + ref->magCal = ref->Mag; // XXX VERY temporary hack to avoid M31 bulge @@ -306,4 +307,5 @@ ref->Color = 0.0; } + ref->magCal = ref->Mag; // XXX VERY temporary hack to avoid M31 bulge Index: trunk/psastro/src/psastroModelAdjust.c =================================================================== --- trunk/psastro/src/psastroModelAdjust.c (revision 38040) +++ trunk/psastro/src/psastroModelAdjust.c (revision 39926) @@ -28,4 +28,6 @@ } + + // if we have not measured the boresite position, no adjustment is needed bool fitBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.FIT.BORESITE"); @@ -52,5 +54,14 @@ return false; } - + float refChipAngleNominal = PS_RAD_DEG*psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.REF.CHIP.ANGLE"); + if (!refChipName) { + psError(PS_ERR_IO, true, "reference chip is missing from recipe"); + return false; + } + + int rotatorParity = psMetadataLookupS32(&status, recipe, "PSASTRO.MODEL.ROT.PARITY"); + if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.ROT.PARITY"); + psMetadataAddS32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.ROT_PARITY", PS_META_REPLACE, "rotator parity parameter", rotatorParity); + // get reference chip from name pmChip *refChip = pmConceptsChipFromName (output->fpa, refChipName); @@ -81,11 +92,18 @@ psFree (PT); } - + // rotate the chip-to-FPA transforms to have 0.0 posangle for refChip; // compensate by rotating fpa to TPA transform // get the current posangle of the ref chip - float chipAngle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]); - fprintf (stderr, "chipAngle: %f\n", chipAngle*PS_DEG_RAD); + // this should have a negative sign : chipAngle = -atan2 (dM/dy, dM/dx) + float chipAngle = -atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]); + + // chipAngle should be refChipAngleNominal @ POSANGLE = 0.0 + float posAngleOffset = rotatorParity * (chipAngle - refChipAngleNominal); + + fprintf (stderr, "chipAngle is: %f, at PA = 0.0, it should be %f, rotating model by %f (parity %d)\n", + chipAngle*PS_DEG_RAD, refChipAngleNominal*PS_DEG_RAD, posAngleOffset*PS_DEG_RAD, rotatorParity); + // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle); @@ -99,9 +117,11 @@ psRegion *region = pmChipPixels (chip); - psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, chipAngle); + // this should ALSO have a negative sign: this rotates by +posAngleOffset + psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posAngleOffset); psFree (chip->toFPA); chip->toFPA = toFPA; // invert the new fromFPA transform to get the new toFPA transform + // NOTE: when we call psPlaneTransformInvert here, we do not increase the order as in other places psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); psFree (chip->fromFPA); @@ -125,4 +145,5 @@ psFree (output->fpa->fromTPA); + // NOTE: when we call psPlaneTransformInvert here, we do not increase the order as in other places output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50); @@ -131,4 +152,5 @@ psMetadata *header = output->fpa->hdu->header; + pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL); @@ -172,4 +194,5 @@ // invert the new fromFPA transform to get the new toFPA transform // the region used here is the region covered by the chip in the FPA + // NOTE: when we call psPlaneTransformInvert here, we do not increase the order as in other places psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); psFree (chip->fromFPA); Index: trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- trunk/psastro/src/psastroMosaicAstrom.c (revision 38040) +++ trunk/psastro/src/psastroMosaicAstrom.c (revision 39926) @@ -15,4 +15,5 @@ bool psastroMosaicFit (pmFPA *fpa, psMetadata *recipe, const char *rootname, int pass); +bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe); // XXX require this fpa to have multiple chip extensions and a PHU? @@ -43,27 +44,36 @@ if (!status) psAbort ("missing config value"); + // if projection is not TAN, fit the measured projection here and combine toTPA/fromTPA functions + bool fitMosaicDistortion = true; + if ((fpa->toSky->type != PS_PROJ_TAN) && (fpa->toSky->type != PS_PROJ_DIS)) { + if (!psastroProjectionRefit (fpa, recipe)) psAbort ("failed to refit distortion"); + fitMosaicDistortion = false; + } + // this should be in a loop with nIter = - for (int iter = 0; iter < nIter; iter++) { + for (int iter = 0; fitMosaicDistortion && (iter < nIter); iter++) { if (!psastroMosaicFit (fpa, recipe, outroot, iter)) return false; } - // now fit the chips under the common distortion with higher-order terms - // first, re-perform the match with a slightly tighter circle - if (!psastroMosaicSetMatch (fpa, recipe, nIter)) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass)"); - return false; - } - if (!psastroMosaicChipAstrom (fpa, recipe, nIter)) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (4th pass)"); - return false; - } - - if (psTraceGetLevel("psastro.dump") > 0) { + for (int iter = 0; (iter < nIter); iter++) { + // now fit the chips under the common distortion with higher-order terms + // first, re-perform the match with a slightly tighter circle + if (!psastroMosaicSetMatch (fpa, recipe, iter)) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass %d)", iter); + return false; + } + if (!psastroMosaicChipAstrom (fpa, recipe, iter)) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (pass %d)", iter); + return false; + } + + if (psTraceGetLevel("psastro.dump") > 0) { // the last filename (see filenames in psastroMosaicFit) char filename[256]; snprintf (filename, 256, "%s.%d.dat", outroot, 2*nIter + 2); psastroDumpMatches (fpa, filename); - } - + } + } + // save WCS and analysis metadata in update header. // (pull or create local view to entry on readout->analysis) @@ -154,2 +164,143 @@ return true; } + +// we have a fpa->toSky projection which is NOT of type TAN along with a toTPA which is +// the Identity transform. we want to convert this to TAN projection plus a non-identity +// transformation to take the non-TAN components. + +// we are going to generate a set of FP->(x,y) using the current projection + fromTPA +// transformations along with a set of TP->(x,y) values using the desired TAN projection, +// then we will fit TP vs FP + +bool psastroProjectionRefit (pmFPA *fpa, psMetadata *recipe) { + + bool status; + + // allocate mosaic-level polynomial transformation and set masks needed by DVO + int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER"); + if (!status) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic distortion fit order\n"); + return false; + } + + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + pmFPAview *view = pmFPAviewAlloc (0); + + psVector *L = psVectorAllocEmpty (1000, PS_TYPE_F32); + psVector *M = psVectorAllocEmpty (1000, PS_TYPE_F32); + + psVector *P = psVectorAllocEmpty (1000, PS_TYPE_F32); + psVector *Q = psVectorAllocEmpty (1000, PS_TYPE_F32); + + psProjection *toSkyTan = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, fpa->toSky->Xs, fpa->toSky->Ys, PS_PROJ_TAN); + + float xMin = NAN; + float xMax = NAN; + float yMin = NAN; + float yMax = NAN; + + bool firstObject = true; + + // this loop selects the matched stars for all chips + 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) continue; + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) continue; + + // process each of the readouts + // XXX there can only be one readout per chip, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) continue; + + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); + if (refstars == NULL) continue; + psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); + + psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); + if (matches == NULL) continue; + + // we are looking over the matched refstars only to be sure we are + // covering the valid space of the projection + transformation, and not + // beyond + + for (int i = 0; i < matches->n; i++) { + pmAstromMatch *match = matches->data[i]; + pmAstromObj *ref = refstars->data[match->ref]; + + psPlane fpOld, tpNew, tpOld; + + psProject (&tpOld, ref->sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + psPlaneTransformApply (&fpOld, fpa->fromTPA, &tpOld); + + psProject (&tpNew, ref->sky, toSkyTan); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + + psVectorAppend (L, fpOld.x); + psVectorAppend (M, fpOld.y); + + psVectorAppend (P, tpNew.x); + psVectorAppend (Q, tpNew.y); + + if (firstObject) { + xMin = xMax = fpOld.x; + yMin = yMax = fpOld.y; + firstObject = false; + } + + xMin = PS_MIN (xMin, fpOld.x); + xMax = PS_MAX (xMax, fpOld.x); + yMin = PS_MIN (yMin, fpOld.y); + yMax = PS_MAX (yMax, fpOld.y); + } + } + } + } + + // the original transforms are (1, 1), but we will need higher order to fit the distortions + psFree (fpa->toTPA); + + // the original transforms are (1, 1), but we will need higher order to fit the distortions + fpa->toTPA = psPlaneTransformAlloc (order, order); + for (int i = 0; i <= fpa->toTPA->x->nX; i++) { + for (int j = 0; j <= fpa->toTPA->x->nY; j++) { + if (i + j > order) { + fpa->toTPA->x->coeffMask[i][j] = PS_POLY_MASK_SET; + fpa->toTPA->y->coeffMask[i][j] = PS_POLY_MASK_SET; + } + } + } + + psVectorFitPolynomial2D (fpa->toTPA->x, NULL, 0, P, NULL, L, M); + psVectorFitPolynomial2D (fpa->toTPA->y, NULL, 0, Q, NULL, L, M); + + psRegion fitRegion = psRegionSet (xMin, xMax, yMin, yMax); + + // psPlaneTransformInvert will generate a new fromTPA with order to match toTPA + psFree (fpa->fromTPA); + psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4); + fpa->fromTPA = psPlaneTransformInvert (myPT, fpa->toTPA, fitRegion, 100); + psFree (myPT); + + psFree (fpa->toSky); + fpa->toSky = toSkyTan; + + psFree (L); + psFree (M); + psFree (P); + psFree (Q); + + psFree (view); + + if (!psastroMosaicSetAstrom (fpa)) { + psError(PSASTRO_ERR_UNKNOWN, false, "failed to apply mosaic distortion terms\n"); + return false; + } + + return true; +} + Index: trunk/psastro/src/psastroMosaicCorrectDistortion.c =================================================================== --- trunk/psastro/src/psastroMosaicCorrectDistortion.c (revision 38040) +++ trunk/psastro/src/psastroMosaicCorrectDistortion.c (revision 39926) @@ -62,5 +62,7 @@ psFree (fpa->fromTPA); - fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50); + psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4); + fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, *region, 50); + psFree (myPT); psFree (region); Index: trunk/psastro/src/psastroMosaicGradients.c =================================================================== --- trunk/psastro/src/psastroMosaicGradients.c (revision 38040) +++ trunk/psastro/src/psastroMosaicGradients.c (revision 39926) @@ -79,5 +79,5 @@ int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER"); if (!status) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n"); + psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic distortion fit order\n"); psFree (gradients); psFree (view); Index: trunk/psastro/src/psastroMosaicOneChip.c =================================================================== --- trunk/psastro/src/psastroMosaicOneChip.c (revision 38040) +++ trunk/psastro/src/psastroMosaicOneChip.c (revision 39926) @@ -100,19 +100,19 @@ } - // XXX allow statitic to be set by the user + // XXX allow statistic to be set by the user // only clip if we are fitting the chip parameters. psStats *fitStats = NULL; -// if (order == 0) { -// fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); -// fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); -// fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); -// } else { - fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); - fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); -// } + if (FALSE && (order == 0)) { + fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); + } else { + fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NSIGMA"); + fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); + } // need to pass in an update header, sent in from above - pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + pmAstromFitResults *results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe); if (!results) { psError(PSASTRO_ERR_DATA, false, "failed to perform the matched fit\n"); Index: trunk/psastro/src/psastroOneChipFit.c =================================================================== --- trunk/psastro/src/psastroOneChipFit.c (revision 38040) +++ trunk/psastro/src/psastroOneChipFit.c (revision 39926) @@ -124,5 +124,5 @@ // improved fit for astrometric terms - results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats); + results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats, recipe); if (!results) { psLogMsg ("psastro", 3, "failed to perform the matched fit\n");