Index: /trunk/psastro/Makefile =================================================================== --- /trunk/psastro/Makefile (revision 7331) +++ /trunk/psastro/Makefile (revision 7332) @@ -39,4 +39,6 @@ $(SRC)/psastroMosaicSetAstrom.$(ARCH).o \ $(SRC)/psastroMosaicSetMatch.$(ARCH).o \ +$(SRC)/psastroMosaicHeaders.$(ARCH).o \ +$(SRC)/psastroMosaicRescaleChips.$(ARCH).o \ $(SRC)/psastroWCS.$(ARCH).o Index: /trunk/psastro/src/psastro.c =================================================================== --- /trunk/psastro/src/psastro.c (revision 7331) +++ /trunk/psastro/src/psastro.c (revision 7332) @@ -26,23 +26,15 @@ psArray *refs = psastroLoadReferences (config); - psastroChipAstrom (config, refs); - psastroMosaicAstrom (config, refs); + psastroMosaicGetRefstars (config, refs); - # if 0 - // perform the desired astrometry analysis - switch (mode) { - case stack: - psastroStackAstrom (config, refs); - break - case chip: + char *mosastro = psMetadataLookupStr (NULL, config->arguments, "MOSASTRO"); + if (mosastro == NULL) { psastroChipAstrom (config, refs); - break - case mosaic: + } else { psastroMosaicAstrom (config, refs); - break; - default: - break; } - # endif + + // XXX how do we specify stack astrometry? + // psastroStackAstrom (config, refs); // write out the results @@ -60,6 +52,6 @@ pmConceptsDone (); pmConfigDone (); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); - // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); + // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro"); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro"); exit (0); Index: /trunk/psastro/src/psastro.h =================================================================== --- /trunk/psastro/src/psastro.h (revision 7331) +++ /trunk/psastro/src/psastro.h (revision 7332) @@ -23,5 +23,6 @@ bool psastroAstromGuess (pmConfig *config); -bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale); +// bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale); +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic); bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale); psPlaneDistort *psPlaneDistortIdentity (); @@ -41,6 +42,11 @@ psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicAstrom (pmConfig *config, psArray *refs); -bool psastroMosaicGetRefstars (pmFPA *fpa, psArray *refs, psMetadata *recipe); +bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs); bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe); bool psastroMosaicSetAstrom (pmFPA *fpa); +bool psastroMosaicHeaders (pmConfig *config); +bool psastroMosaicRescaleChips (pmFPA *fpa); + +bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale); +psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale); Index: /trunk/psastro/src/psastroArguments.c =================================================================== --- /trunk/psastro/src/psastroArguments.c (revision 7331) +++ /trunk/psastro/src/psastroArguments.c (revision 7332) @@ -45,4 +45,11 @@ } + // mosastro mode also specifies output header file + if ((N = psArgumentGet (*argc, argv, "-mosastro"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MOSASTRO", PS_DATA_STRING, "", psStringCopy(argv[N])); + psArgumentRemove (N, argc, argv); + } + status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list"); if (!status) { usage ();} Index: /trunk/psastro/src/psastroAstromGuess.c =================================================================== --- /trunk/psastro/src/psastroAstromGuess.c (revision 7331) +++ /trunk/psastro/src/psastroAstromGuess.c (revision 7332) @@ -3,8 +3,11 @@ // XXX this function assumes the initial astrometry arrives in the form of // XXX WCS keywords in the headers corresponding to the chips. +// XXX this function is both converting the header WCS astrometry terms to the fpa terms +// and applying the astrometry to the detected objects. bool psastroAstromGuess (pmConfig *config) { bool newFPA = true; bool status = false; + bool isMosaic = false; double RAmin, RAmax, RAminSky, RAmaxSky; double DECmin, DECmax; @@ -13,4 +16,7 @@ pmCell *cell = NULL; pmReadout *readout = NULL; + + // is this a mosaic astrometry analysis? + psMetadataLookupStr (&isMosaic, config->arguments, "MOSASTRO"); // select the current recipe @@ -30,4 +36,5 @@ double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE"); if (!status) plateScale = 1.0; + plateScale = 1.0; pmFPAview *view = pmFPAviewAlloc (0); @@ -41,5 +48,5 @@ pmHDU *hdu = pmFPAviewThisHDU (view, fpa); - pmAstromReadWCS (&fpa->projection, &fpa->toTangentPlane, &chip->toFPA, hdu->header, plateScale); + pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic); if (newFPA) { newFPA = false; @@ -49,10 +56,4 @@ DECmin = DECmax = fpa->projection->D; } - - // build projection inversions - // XXX region should be set from image data - psRegion region = {0, 0, 2000, 4500}; - fpa->fromTangentPlane = psPlaneDistortIdentity (); - chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20); // apply the new WCS guess data to all of the data in the readouts @@ -66,23 +67,45 @@ if (! readout->data_exists) { continue; } - psArray *objects = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); - if (objects == NULL) { continue; } + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) { continue; } - for (int i = 0; i < objects->n; i++) { - pmAstromObj *object = objects->data[i]; + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; - psPlaneTransformApply (object->FP, chip->toFPA, object->chip); - psPlaneDistortApply (object->TP, fpa->toTangentPlane, object->FP, 0.0, 0.0); - p_psDeproject (object->sky, object->TP, fpa->projection); + psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); + psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); + p_psDeproject (raw->sky, raw->TP, fpa->projection); + + if (i < 0) { + fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", + raw->chip->x, raw->chip->y, + raw->FP->x, raw->FP->y, + raw->TP->x, raw->TP->y, + raw->sky->r, raw->sky->d); + + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + psPlane *ch = psPlaneAlloc(); + + p_psProject (tp, raw->sky, fpa->projection); + psPlaneDistortApply (fp, fpa->fromTangentPlane, tp, 0.0, 0.0); + psPlaneTransformApply (ch, chip->fromFPA, fp); + + fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n", + ch->x, ch->y, + fp->x, fp->y, + tp->x, tp->y, + raw->sky->r, raw->sky->d); + } // rationalize ra to sky range centered on boresite - while (object->sky->r < RAminSky) object->sky->r += M_PI; - while (object->sky->r > RAmaxSky) object->sky->r -= M_PI; + while (raw->sky->r < RAminSky) raw->sky->r += M_PI; + while (raw->sky->r > RAmaxSky) raw->sky->r -= M_PI; - RAmin = PS_MIN (object->sky->r, RAmin); - RAmax = PS_MAX (object->sky->r, RAmax); + RAmin = PS_MIN (raw->sky->r, RAmin); + RAmax = PS_MAX (raw->sky->r, RAmax); - DECmin = PS_MIN (object->sky->d, DECmin); - DECmax = PS_MAX (object->sky->d, DECmax); + DECmin = PS_MIN (raw->sky->d, DECmin); + DECmax = PS_MAX (raw->sky->d, DECmax); } } Index: /trunk/psastro/src/psastroChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroChipAstrom.c (revision 7331) +++ /trunk/psastro/src/psastroChipAstrom.c (revision 7332) @@ -41,30 +41,12 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); if (rawstars == NULL) { continue; } - // the refstars is a subset within range of this chip - psArray *refstars = psArrayAlloc (100); + // select the raw objects for this readout + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } - // select the reference objects within range of this readout - // project the reference objects to this chip - for (int i = 0; i < refs->n; i++) { - pmAstromObj *ref = refs->data[i]; - - p_psProject (ref->TP, ref->sky, fpa->projection); - psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0); - psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); - - // XXX what is the X,Y range of the selected chip? - // XXX for the moment, force these to be good for Megacam - if (ref->chip->x < -500) continue; - if (ref->chip->x > 2500) continue; - if (ref->chip->y < -500) continue; - if (ref->chip->y > 5000) continue; - - psArrayAdd (refstars, 100, ref); - } psastroOneChip (fpa, chip, refstars, rawstars, recipe); - psFree (refstars); // read WCS data from the corresponding header Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 7331) +++ /trunk/psastro/src/psastroConvert.c (revision 7332) @@ -9,8 +9,8 @@ return false; - psArray *objects = pmSourceToAstromObj (sources); + psArray *rawstars = pmSourceToAstromObj (sources); - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.OBJECTS", PS_DATA_ARRAY, "astrometry objects", objects); - psFree (objects); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars); + psFree (rawstars); return true; Index: /trunk/psastro/src/psastroLoadReferences.c =================================================================== --- /trunk/psastro/src/psastroLoadReferences.c (revision 7331) +++ /trunk/psastro/src/psastroLoadReferences.c (revision 7332) @@ -62,7 +62,7 @@ // convert the Average table to the pmAstromObj entries - psArray *objects = psArrayAlloc (table->n); + psArray *refs = psArrayAlloc (table->n); for (int i = 0; i < table->n; i++) { - pmAstromObj *object = pmAstromObjAlloc (); + pmAstromObj *ref = pmAstromObjAlloc (); psMetadata *row = table->data[i]; @@ -70,15 +70,15 @@ // DVO tables are stored in degrees # if ELIXIR_MODE - object->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); - object->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); - object->Mag = 0.001*psMetadataLookupS32 (&status, row, "MAG"); // XXX ELIXIR uses millimags, PANSTARRS uses mags + ref->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); + ref->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); + ref->Mag = 0.001*psMetadataLookupS32 (&status, row, "MAG"); // XXX ELIXIR uses millimags, PANSTARRS uses mags # else - object->sky->r = RAD_DEG*psMetadataLookupF64 (&status, row, "RA"); - object->sky->d = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC"); - object->Mag = psMetadataLookupF32 (&status, row, "MAG"); + ref->sky->r = RAD_DEG*psMetadataLookupF64 (&status, row, "RA"); + ref->sky->d = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC"); + ref->Mag = psMetadataLookupF32 (&status, row, "MAG"); # endif - psArrayAdd (objects, 100, object); - psFree (object); + psArrayAdd (refs, 100, ref); + psFree (ref); } psFree (header); @@ -86,6 +86,6 @@ psTrace (__func__, 3, "loaded %d reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n", - objects->n, RAmin, DECmin, RAmax, DECmax); - return objects; + refs->n, RAmin, DECmin, RAmax, DECmax); + return refs; } Index: /trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicAstrom.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicAstrom.c (revision 7332) @@ -20,11 +20,26 @@ } - double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE"); - if (!status) plateScale = 1.0; + int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.DISTORTION.ORDER"); + if (!status) order = 3.0; pmFPA *fpa = input->fpa; - psastroMosaicGetRefstars (fpa, refs, recipe); + // XXX test point + if (0) { + pmChip *chip = fpa->chips->data[10]; + pmCell *cell = chip->cells->data[0]; + pmReadout *readout = cell->readouts->data[0]; + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + pmAstromObj *ref = refstars->data[0]; + fprintf (stderr, "test 1: %f %f %f %f %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d); + } + psastroMosaicSetMatch (fpa, recipe); + + // XXX no input distortion model yet; existing fpa distortion is identity: replace it + psFree (fpa->toTangentPlane); + psFree (fpa->fromTangentPlane); + fpa->toTangentPlane = psPlaneDistortIdentity (order); + fpa->fromTangentPlane = psPlaneDistortIdentity (order); grads = psastroMosaicGetGrads (fpa, recipe); @@ -32,4 +47,5 @@ // fit the measured gradients with the telescope distortion model (3rd order polynomial) pmAstromFitDistortion (fpa, grads, recipe); + psastroMosaicRescaleChips (fpa); // apply the new distortion terms up and down @@ -50,5 +66,5 @@ // re-fit the chips with higher-order fits psastroMosaicChipAstrom (fpa, recipe); - // psastroMosaicHeaders (fpa); + psastroMosaicHeaders (config); return true; Index: /trunk/psastro/src/psastroMosaicChipAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicChipAstrom.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicChipAstrom.c (revision 7332) @@ -23,5 +23,5 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); if (rawstars == NULL) { continue; } Index: /trunk/psastro/src/psastroMosaicGetGrads.c =================================================================== --- /trunk/psastro/src/psastroMosaicGetGrads.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicGetGrads.c (revision 7332) @@ -25,5 +25,5 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); if (rawstars == NULL) { continue; } Index: /trunk/psastro/src/psastroMosaicGetRefstars.c =================================================================== --- /trunk/psastro/src/psastroMosaicGetRefstars.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicGetRefstars.c (revision 7332) @@ -1,10 +1,30 @@ # include "psastro.h" -bool psastroMosaicGetRefstars (pmFPA *fpa, psArray *refs, psMetadata *recipe) { +bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs) { + bool status; pmChip *chip = NULL; pmCell *cell = NULL; pmReadout *readout = NULL; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n"); + exit(EXIT_FAILURE); + } + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!input) { + psErrorStackPrint(stderr, "Can't find input data!\n"); + exit(EXIT_FAILURE); + } + + float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA"); + if (!status) fieldExtra = 0.0; + pmFPAview *view = pmFPAviewAlloc (0); + pmFPA *fpa = input->fpa; // this loop selects the matched stars for all chips @@ -22,15 +42,24 @@ if (! readout->data_exists) { continue; } - // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); - if (rawstars == NULL) { continue; } + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1"); + int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2"); + + float minX = -fieldExtra*Nx; + float maxX = (1+fieldExtra)*Nx; + float minY = -fieldExtra*Ny; + float maxY = (1+fieldExtra)*Ny; + // the refstars is a subset within range of this chip psArray *refstars = psArrayAlloc (100); + + int Nst = 0; // select the reference objects within range of this readout // project the reference objects to this chip for (int i = 0; i < refs->n; i++) { - pmAstromObj *ref = refs->data[i]; + pmAstromObj *ref = pmAstromObjCopy(refs->data[i]); p_psProject (ref->TP, ref->sky, fpa->projection); @@ -38,17 +67,55 @@ psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); - // XXX what is the X,Y range of the selected chip? - // XXX for the moment, force these to be good for Megacam - if (ref->chip->x < -10) continue; - if (ref->chip->x > 2060) continue; - if (ref->chip->y < -10) continue; - if (ref->chip->y > 4610) continue; + // limit the X,Y range of the refs to the selected chip + if (ref->chip->x < minX) continue; + if (ref->chip->x > maxX) continue; + if (ref->chip->y < minY) continue; + if (ref->chip->y > maxY) continue; + if (Nst < 0) { + fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n", + ref->chip->x, ref->chip->y, + ref->FP->x, ref->FP->y, + ref->TP->x, ref->TP->y, + ref->sky->r, ref->sky->d); + + psSphere *sk = psSphereAlloc(); + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + + psPlaneTransformApply (fp, chip->toFPA, ref->chip); + psPlaneDistortApply (tp, fpa->toTangentPlane, fp, 0.0, 0.0); + p_psDeproject (sk, tp, fpa->projection); + + fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", + ref->chip->x, ref->chip->y, + fp->x, fp->y, + tp->x, tp->y, + sk->r, sk->d); + } + Nst ++; + psArrayAdd (refstars, 100, ref); } + psTrace (__func__, 4, "Added %d refstars\n", refstars->n); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); + + // XXX test point + if (0) { + pmChip *chip1 = fpa->chips->data[10]; + pmCell *cell1 = chip1->cells->data[0]; + pmReadout *readout1 = cell1->readouts->data[0]; + psArray *refstars1 = psMetadataLookupPtr (NULL, readout1->analysis, "PSASTRO.REFSTARS"); + pmAstromObj *ref1 = refstars1->data[0]; + fprintf (stderr, "test 3: %f %f %f %f %f %f\n", ref1->chip->x, ref1->chip->y, ref1->FP->x, ref1->FP->y, ref1->sky->r, ref1->sky->d); + } + + psFree (refstars); } } } + + psFree (view); return true; } Index: /trunk/psastro/src/psastroMosaicSetAstrom.c =================================================================== --- /trunk/psastro/src/psastroMosaicSetAstrom.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicSetAstrom.c (revision 7332) @@ -23,5 +23,5 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); if (rawstars == NULL) { continue; } Index: /trunk/psastro/src/psastroMosaicSetMatch.c =================================================================== --- /trunk/psastro/src/psastroMosaicSetMatch.c (revision 7331) +++ /trunk/psastro/src/psastroMosaicSetMatch.c (revision 7332) @@ -7,4 +7,10 @@ pmReadout *readout = NULL; pmFPAview *view = pmFPAviewAlloc (0); + + FILE *f; + char name[128]; + + FILE *g1 = fopen ("raw.ps.dat", "w"); + FILE *g2 = fopen ("ref.ps.dat", "w"); // this loop selects the matched stars for all chips @@ -23,5 +29,5 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); if (rawstars == NULL) { continue; } @@ -29,14 +35,52 @@ psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); if (refstars == NULL) { continue; } + psTrace (__func__, 4, "Trying %d refstars\n", refstars->n); // use small radius to match stars (assume starting astrometry is good) // XXX should this take a (double radius)? - psArray *matches = pmAstromRadiusMatch (rawstars, refstars, recipe); + psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, recipe); + psTrace (__func__, 4, "Matched %d refstars\n", matches->n); + + sprintf (name, "raw.%02d.dat", view->chip); + f = fopen (name, "w"); + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + fprintf (f, "%f %f %f %f %f %f\n", raw->chip->x, raw->chip->y, raw->FP->x, raw->FP->y, raw->sky->r, raw->sky->d); + } + fclose (f); + + sprintf (name, "ref.%02d.dat", view->chip); + f = fopen (name, "w"); + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + fprintf (f, "%f %f %f %f %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d); + } + fclose (f); + + for (int i = 0; i < matches->n; i++) { + pmAstromMatch *match = matches->data[i]; + + pmAstromObj *raw = rawstars->data[match->raw]; + fprintf (g1, "%d %f %f %f %f %f %f %f %f\n", i, + DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d, + raw->TP->x, raw->TP->y, + raw->FP->x, raw->FP->y, + raw->chip->x, raw->chip->y); + + pmAstromObj *ref = refstars->data[match->ref]; + fprintf (g2, "%d %f %f %f %f %f %f %f %f\n", i, + DEG_RAD*ref->sky->r, DEG_RAD*ref->sky->d, + ref->TP->x, ref->TP->y, + ref->FP->x, ref->FP->y, + ref->chip->x, ref->chip->y); + } // XXX drop the old one - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY, "astrometry matches", matches); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches); } } } + fclose (g1); + fclose (g2); return true; } Index: /trunk/psastro/src/psastroWCS.c =================================================================== --- /trunk/psastro/src/psastroWCS.c (revision 7331) +++ /trunk/psastro/src/psastroWCS.c (revision 7332) @@ -3,9 +3,6 @@ // interpret header WCS (only handles traditional WCS for the moment) // plateScale is nominal physical scale on tangent plane (microns / arcsecond) -bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) { - - psPlaneTransform *toFPA; - psPlaneDistort *toTPA; - psProjection *toSky; +bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) { + psProjectionType type; bool status; @@ -13,5 +10,5 @@ float pc1_1, pc1_2, pc2_1, pc2_2; - // *** interpret header data, convert to crval(i), etc + // interpret header data, convert to crval(i), etc char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2"); if (!status) { @@ -88,51 +85,109 @@ got_matrix: - // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the - // toFPA must be modified to match the crval(i), scale(i) and crpix(i) - // scale = scale(i)/scale(0) (i == chip #) - // project crval1(0),crval2(0 using projection - toFPA = psPlaneTransformAlloc (1, 1); + /***** + + For mosaic astrometry, we need to have a starting set of projection terms in which the + chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or + microns). This set of projections, coupled with an identity toTPA (ie, no distortion) will + result in substantial errors between the observed and predicted star positions on the focal + plane: this is the measurement of the optical distortion in the camera. At the same time, + we need to carry around the transformations which allow us to make an accurate calculation + of the position of the stars based on the input (per-chip) astrometry. These + transformations will allow us to match the raw and ref stars robustly. To convert the + per-chip astrometry (which may have been calculated with a different plate scale for each + chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust + the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for + each chip (eg, cdelt1). Thus, we need to carry around both the + + *****/ + + { + double rX = 1.0; + double rY = 1.0; + + // XXX free an existing toFPA? + psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1); - double cdelt = hypot (cdelt1, cdelt2) / plateScale; // degrees / micron (eg, in fact, whatever unit user chooses for focal plane) - cdelt1 /= cdelt; - cdelt2 /= cdelt; - - toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2); - toFPA->x->coeff[1][0] = +(pc1_1*cdelt1); - toFPA->x->coeff[0][1] = +(pc1_2*cdelt2); - toFPA->x->mask[1][1] = 1; - - toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2); - toFPA->y->coeff[1][0] = +(pc2_1*cdelt1); - toFPA->y->coeff[0][1] = +(pc2_2*cdelt2); - toFPA->y->mask[1][1] = 1; - *toFPAout = toFPA; - - if (*toSkyOut != NULL) { - if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined"); - - psSphere sky; - psPlane tpa; - - sky.r = crval1*RAD_DEG; - sky.d = crval2*RAD_DEG; - p_psProject (&tpa, &sky, *toSkyOut); + // basic transformation from chip to FPA + toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2); + toFPA->x->coeff[1][0] = pc1_1; + toFPA->x->coeff[0][1] = pc1_2; + toFPA->x->mask[1][1] = 1; + + toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2); + toFPA->y->coeff[1][0] = pc2_1; + toFPA->y->coeff[0][1] = pc2_2; + toFPA->y->mask[1][1] = 1; + + // projection from TPA to SKY + psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type); + + if (fpa->toSky == NULL) { + // XXX for now, use the identity for TPA <--> FPA + fpa->toTangentPlane = psPlaneDistortIdentity (1); + fpa->fromTangentPlane = psPlaneDistortIdentity (1); + fpa->toSky = toSky; + } else { + if (fpa->toTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined"); + if (fpa->fromTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined"); + + // adjust for common toSky for mosaic: + // ignore the TPA since toTPA is identity + // find the FPA coordinate of 0,0 for this chip. + psPlane *fp = psPlaneAlloc(); + psPlane *chip = psPlaneAlloc(); + psSphere *sky = psSphereAlloc(); + chip->x = chip->y = 0; - // XXX for the moment, assume toTPA is the identity transformation - toFPA->x->coeff[0][0] = tpa.x; - toFPA->y->coeff[0][0] = tpa.y; - } else { - // XXX for now, use the identity for TPA <--> FPA - toTPA = psPlaneDistortIdentity (); - - // center of projection is (0,0) coordinate of TPA - toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type); - - *toTPAout = toTPA; - *toSkyOut = toSky; + psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate + p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate + p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection + + toFPA->x->coeff[0][0] = fp->x; + toFPA->y->coeff[0][0] = fp->y; + + rX = toSky->Xs/fpa->toSky->Xs; + rY = toSky->Ys/fpa->toSky->Ys; + + // correct to common plate scale + toFPA->x->coeff[1][0] *= rX; + toFPA->x->coeff[0][1] *= rX; + toFPA->y->coeff[1][0] *= rY; + toFPA->y->coeff[0][1] *= rY; + + psFree (fp); + psFree (sky); + psFree (chip); + psFree (toSky); + } + + chip->toFPA = toFPA; + chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA); + + // remove the correction to the common plate scale + if (isMosaic) { + chip->toFPA->x->coeff[1][0] /= rX; + chip->toFPA->x->coeff[0][1] /= rX; + chip->toFPA->y->coeff[1][0] /= rY; + chip->toFPA->y->coeff[0][1] /= rY; + } + + fprintf (stderr, "toFPA: %f %f (%f,%f),(%f,%f)\n", + chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0], + chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1], + chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]); + + fprintf (stderr, "frFPA: %f %f (%f,%f),(%f,%f)\n", + chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0], + chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1], + chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]); + + fprintf (stderr, "field: %f,%f %f,%f\n", + DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D, + 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys); + } return true; } - // convert toFPA / toSky components to traditional WCS @@ -169,11 +224,11 @@ psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD); - psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD); - - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]); - psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale); // alternative representations use @@ -184,18 +239,133 @@ } -psPlaneDistort *psPlaneDistortIdentity () { +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) { + + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP"); + + // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toFPA->x->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toFPA->y->coeff[0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", 1.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", 1.0); + + if (toFPA->x->nX != toFPA->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)"); + if (toFPA->x->nX != toFPA->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)"); + if (toFPA->x->nX != toFPA->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)"); + + switch (toFPA->x->nX) { + case 3: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]); + + case 2: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]); + + case 1: + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale); + break; + + case 0: + psAbort ("psastro", "invalid tangent plane order"); + } + return true; +} + +// convert toFPA / toSky components to traditional WCS +// plateScale is nominal physical scale on tangent plane (microns / arcsecond) +psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) { + + psMetadata *header = psMetadataAlloc (); + + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS"); + psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS"); + + // XXX need to handle the plateScale?? + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", PS_META_REPLACE, "", toSky->R*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2", PS_META_REPLACE, "", toSky->D*DEG_RAD); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", 0.0); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale); + psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale); + + if (toTP->x->nX != toTP->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)"); + if (toTP->x->nX != toTP->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)"); + if (toTP->x->nX != toTP->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)"); + + switch (toTP->x->nX) { + case 3: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]); + + case 2: + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]); + + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]); + + case 1: + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]); + psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]); + break; + + case 0: + psAbort ("psastro", "invalid tangent plane order"); + } + return header; +} + +// XXX for the moment, we mimic the limited Elixir mosastro orders +psPlaneDistort *psPlaneDistortIdentity (int order) { psPlaneDistort *distort; - distort = psPlaneDistortAlloc (1, 1, 0, 0); - distort->x->coeff[0][0][0][0] = 0; + if (order < 1) psAbort ("psastro", "invalid order"); + if (order > 3) psAbort ("psastro", "invalid order"); + + // all coeffs and masks initially set to 0 + distort = psPlaneDistortAlloc (order, order, 0, 0); + + for (int i = 0; i <= order; i++) { + for (int j = 0; j <= order; j++) { + if (i + j > order) { + distort->x->mask [1][1][0][0] = 1; + distort->y->mask [1][1][0][0] = 1; + } + } + } distort->x->coeff[1][0][0][0] = 1; - distort->x->coeff[0][1][0][0] = 0; - distort->x->mask [1][1][0][0] = 1; - - distort->y->coeff[0][0][0][0] = 0; - distort->y->coeff[1][0][0][0] = 0; distort->y->coeff[0][1][0][0] = 1; - distort->y->mask [1][1][0][0] = 1; return distort;