Index: branches/eam_branches/ipp-20110213/Ohana/src/addstar/src/ReadStarsFITS.c =================================================================== --- branches/eam_branches/ipp-20110213/Ohana/src/addstar/src/ReadStarsFITS.c (revision 30762) +++ branches/eam_branches/ipp-20110213/Ohana/src/addstar/src/ReadStarsFITS.c (revision 30812) @@ -45,4 +45,8 @@ } if (!strcmp (type, "PS1_V2")) { + if (table.header[0].Naxis[0] == 20) { + // skip the invalid DETEFF tables which were mistakenly labeled as PS1_V2 + return (NULL); + } stars = Convert_PS1_V2 (&table, &Nstars); } @@ -51,5 +55,5 @@ } if (stars == NULL) { - fprintf (stderr, "ERROR: invalid table type %s\n", type); + fprintf (stderr, "invalid table type %s\n", type); return (NULL); } @@ -212,4 +216,85 @@ Stars *stars; CMF_PS1_V1 *ps1data; + + // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec. this was a + // mistake in two ways: a new format should have been defined (eg, CMF_PS1_V2), and the + // layout used did not have clean byte-boundaries for the corresponding structure. The + // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the + // autocode tools do not work to read in the new version, even if we recognize it. Here we + // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and + // call a special conversion function if it is found. + + if (table[0].header[0].Naxis[0] == 136) { + stars = Convert_PS1_V1_Alt (table, nstars); + return (stars); + } + + ps1data = gfits_table_get_CMF_PS1_V1 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } + ZeroPt = GetZeroPoint(); + + ALLOCATE (stars, Stars, Nstars); + for (i = 0; i < Nstars; i++) { + InitStar (&stars[i]); + stars[i].measure.Xccd = ps1data[i].X; + stars[i].measure.Yccd = ps1data[i].Y; + stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); + stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); + + stars[i].measure.posangle = ToShortDegrees(ps1data[i].posangle); + stars[i].measure.pltscale = ps1data[i].pltscale; + + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + } + stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = ps1data[i].dMcal; + stars[i].measure.Map = ps1data[i].Map + ZeroPt; + + stars[i].measure.Sky = ps1data[i].sky; + stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQual = ps1data[i].psfQual; + stars[i].measure.psfNdof = ps1data[i].psfNdof; + stars[i].measure.psfNpix = ps1data[i].psfNpix; + stars[i].measure.crNsigma = ps1data[i].crNsigma; + stars[i].measure.extNsigma = ps1data[i].extNsigma; + + stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); + stars[i].measure.FWy = ToShortPixels(ps1data[i].fy); + stars[i].measure.theta = ToShortDegrees(ps1data[i].df); + + stars[i].measure.Mxx = ToShortPixels(ps1data[i].Mxx); + stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); + stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); + + stars[i].measure.photFlags = ps1data[i].flags; + + // this is may optionally be replaced by the internal sequence (see FilterStars.c) + stars[i].measure.detID = ps1data[i].detID; + + // the Average fields and the following Measure fields are set in FilterStars after + // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID, + // averef is set in find_matches, dbFlags is zero on ingest. + + // the following fields are currently not being set anywhere: t_msec + } + *nstars = Nstars; + return (stars); +} + +Stars *Convert_PS1_V1_Alt (FTable *table, unsigned int *nstars) { + + off_t Nstars; + unsigned int i; + double ZeroPt; + Stars *stars; + CMF_PS1_V2 *ps1data; // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec. this was a @@ -221,10 +306,5 @@ // call a special conversion function if it is found. - if (table[0].header[0].Naxis[0] == 136) { - stars = Convert_PS1_V1_Alt (table, nstars); - return (stars); - } - - ps1data = gfits_table_get_CMF_PS1_V1 (table, &Nstars, NULL); + ps1data = gfits_table_get_CMF_PS1_V1_Alt (table, &Nstars, NULL); if (!ps1data) { fprintf (stderr, "skipping inconsistent entry\n"); @@ -286,5 +366,5 @@ } -Stars *Convert_PS1_V1_Alt (FTable *table, unsigned int *nstars) { +Stars *Convert_PS1_V2 (FTable *table, unsigned int *nstars) { off_t Nstars; @@ -294,13 +374,5 @@ CMF_PS1_V2 *ps1data; - // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec. this was a - // mistake in two ways: a few format should have been defined (eg, CMF_PS1_V2), and the - // layout used did not have clean byte-boundaries for the corresponding structure. The - // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the - // autocode tools do not work to read in the new version, even if we recognize it. Here we - // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and - // call a special conversion function if it is found. - - ps1data = gfits_table_get_CMF_PS1_V1_Alt (table, &Nstars, NULL); + ps1data = gfits_table_get_CMF_PS1_V2 (table, &Nstars, NULL); if (!ps1data) { fprintf (stderr, "skipping inconsistent entry\n"); @@ -362,5 +434,5 @@ } -Stars *Convert_PS1_V2 (FTable *table, unsigned int *nstars) { +Stars *Convert_PS1_V3 (FTable *table, unsigned int *nstars) { off_t Nstars; @@ -368,7 +440,7 @@ double ZeroPt; Stars *stars; - CMF_PS1_V2 *ps1data; - - ps1data = gfits_table_get_CMF_PS1_V2 (table, &Nstars, NULL); + CMF_PS1_V3 *ps1data; + + ps1data = gfits_table_get_CMF_PS1_V3 (table, &Nstars, NULL); if (!ps1data) { fprintf (stderr, "skipping inconsistent entry\n"); @@ -429,70 +501,2 @@ return (stars); } - -Stars *Convert_PS1_V3 (FTable *table, unsigned int *nstars) { - - off_t Nstars; - unsigned int i; - double ZeroPt; - Stars *stars; - CMF_PS1_V3 *ps1data; - - ps1data = gfits_table_get_CMF_PS1_V3 (table, &Nstars, NULL); - if (!ps1data) { - fprintf (stderr, "skipping inconsistent entry\n"); - return (NULL); - } - ZeroPt = GetZeroPoint(); - - ALLOCATE (stars, Stars, Nstars); - for (i = 0; i < Nstars; i++) { - InitStar (&stars[i]); - stars[i].measure.Xccd = ps1data[i].X; - stars[i].measure.Yccd = ps1data[i].Y; - stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); - stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); - - stars[i].measure.posangle = ToShortDegrees(ps1data[i].posangle); - stars[i].measure.pltscale = ps1data[i].pltscale; - - if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; - } - stars[i].measure.dM = ps1data[i].dM; - stars[i].measure.dMcal = ps1data[i].dMcal; - stars[i].measure.Map = ps1data[i].Map + ZeroPt; - - stars[i].measure.Sky = ps1data[i].sky; - stars[i].measure.dSky = ps1data[i].dSky; - - stars[i].measure.psfChisq = ps1data[i].psfChisq; - stars[i].measure.psfQual = ps1data[i].psfQual; - stars[i].measure.psfNdof = ps1data[i].psfNdof; - stars[i].measure.psfNpix = ps1data[i].psfNpix; - stars[i].measure.crNsigma = ps1data[i].crNsigma; - stars[i].measure.extNsigma = ps1data[i].extNsigma; - - stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); - stars[i].measure.FWy = ToShortPixels(ps1data[i].fy); - stars[i].measure.theta = ToShortDegrees(ps1data[i].df); - - stars[i].measure.Mxx = ToShortPixels(ps1data[i].Mxx); - stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); - stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - - stars[i].measure.photFlags = ps1data[i].flags; - - // this is may optionally be replaced by the internal sequence (see FilterStars.c) - stars[i].measure.detID = ps1data[i].detID; - - // the Average fields and the following Measure fields are set in FilterStars after - // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID, - // averef is set in find_matches, dbFlags is zero on ingest. - - // the following fields are currently not being set anywhere: t_msec - } - *nstars = Nstars; - return (stars); -} Index: branches/eam_branches/ipp-20110213/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- branches/eam_branches/ipp-20110213/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 30762) +++ branches/eam_branches/ipp-20110213/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 30812) @@ -219,8 +219,14 @@ } if (fabs(output[0].measure[Nmeas].dR) > 10*RADIUS) { - fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", - output[0].average[n].R, output[0].average[n].D, Rin, Din, - X1[i], X2[Jmin], Y1[i], Y2[Jmin]); - // XXX abort on this? -- this is a bad failure... + // ok take declination into account and check again. + double cosD = cos(RAD_DEG*Din); + if (fabs(output[0].measure[Nmeas].dR*cosD) > 10*RADIUS) { + + fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", + output[0].average[n].R, output[0].average[n].D, Rin, Din, + X1[i], X2[Jmin], Y1[i], Y2[Jmin]); + + // XXX abort on this? -- this is a bad failure... + } } input[0].found[N] = Nmeas; Index: branches/eam_branches/ipp-20110213/Ohana/src/getstar/src/MatchCoords.c =================================================================== --- branches/eam_branches/ipp-20110213/Ohana/src/getstar/src/MatchCoords.c (revision 30762) +++ branches/eam_branches/ipp-20110213/Ohana/src/getstar/src/MatchCoords.c (revision 30812) @@ -3,4 +3,6 @@ # define FLT_MAX 1e32 # endif +int ComputeLMLimits(double xmin, double xmax, double ymin, double ymax, + double *Lmin, double *Lmax, double *Mmin, double *Mmax, Coords *coords); /* given coordinate, find images in list that contain the point */ @@ -12,4 +14,5 @@ double Xi[4], Yi[4]; /* image and original corners */ double xmin, xmax, ymin, ymax; + double Lmin, Lmax, Mmin, Mmax; @@ -20,5 +23,7 @@ also define the vtable entries for the images we keep */ + int haveLimits = 0; for (i = 0; i < NdbImages; i++) { + haveLimits = 0; if (!WITH_PHU && !strcmp (&dbImages[i].coords.ctype[4], "-DIS")) continue; @@ -30,5 +35,5 @@ /* define image corners */ SetImageCorners (Xi, Yi, &dbImages[i]); - // Xi[4] = Xi[0]; Yi[4] = Yi[0]; + // Xi[3] = Xi[0]; Yi[3] = Yi[0]; ymin = xmin = +FLT_MAX; @@ -52,9 +57,38 @@ continue; } - - if (VERBOSE) fprintf(stderr, "%ld %8.1f %8.1f %s\n", i, x, y, dbImages[i].name); + if (VERBOSE) fprintf(stderr, OFF_T_FMT" %8.1f %8.1f %s\n", i, x, y, dbImages[i].name); if ((x >= xmin && x <= xmax) && (y >= ymin && y <= ymax)) { + if (VERBOSE) fprintf(stderr, "\tpotential overlap for point %d image %i\n", j, (int) i); + + // We may be using the astrometry far outside it's range of validity. + // As a sanity check make sure that the focal plane coords are in bounds. + // See ticket 1381 + if (!haveLimits) { + // We defer computing the focal plane limits until we need them + status = ComputeLMLimits(xmin, xmax, ymin, ymax, + &Lmin, &Lmax, &Mmin, &Mmax, &dbImages[i].coords); + if (!status) { + continue; + } + haveLimits = 1; + } + + double L, M; + status = RD_to_LM(&L, &M, pt->ra, pt->dec, &dbImages[i].coords); + if (!status) { + // This won't happen because RD_to_XY succeeded above + continue; + } + + if (VERBOSE) fprintf(stderr, "L: %f M: %f\n", L, M); + + if (L < Lmin || L > Lmax || M < Mmin || M > Mmax) { + if (VERBOSE) fprintf(stderr, "\trejecting match for because focal plane coordinates are out of range\n"); + continue; + } + if (VERBOSE) fprintf(stderr, "\tmatch confirmed\n"); + totalMatches++; @@ -88,2 +122,36 @@ } } + +// compute bounding box of this image in focal plane coordinates +int ComputeLMLimits(double xmin, double xmax, double ymin, double ymax, + double *Lmin, double *Lmax, double *Mmin, double *Mmax, Coords *coords) +{ + double L0, L1, M0, M1; + int status = XY_to_LM(&L0, &M0, xmin, ymin, coords); + if (!status) { + fprintf(stderr, "failed to transform image corner to LM\n"); + return 0; + } + status = XY_to_LM(&L1, &M1, xmax, ymax, coords); + if (!status) { + fprintf(stderr, "failed to transform image corner to LM\n"); + return 0; + } + + if (L1 > L0) { + *Lmin = L0; + *Lmax = L1; + } else { + *Lmin = L1; + *Lmax = L0; + } + if (M1 > M0) { + *Mmin = M0; + *Mmax = M1; + } else { + *Mmin = M1; + *Mmax = M0; + } + if (VERBOSE) fprintf(stderr, "Lmin: %f Lmax: %f Mmin: %f Mmax: %f\n", *Lmin, *Lmax, *Mmin, *Mmax); + return 1; +} Index: branches/eam_branches/ipp-20110213/Ohana/src/libohana/include/ohana_allocate.h =================================================================== --- branches/eam_branches/ipp-20110213/Ohana/src/libohana/include/ohana_allocate.h (revision 30762) +++ branches/eam_branches/ipp-20110213/Ohana/src/libohana/include/ohana_allocate.h (revision 30812) @@ -49,5 +49,5 @@ # define REALLOCATE(PTR,TYPE,SIZE) { \ - PTR = (TYPE *) realloc(PTR,(unsigned)(MAX(((SIZE)*((int)sizeof(TYPE))),1))); \ + PTR = (TYPE *) realloc(PTR,(size_t)(MAX(((SIZE)*((int)sizeof(TYPE))),1))); \ if (PTR == NULL) { \ fprintf(stderr,"failed realloc at %d in %s\n", __LINE__, __FILE__); \ @@ -57,5 +57,5 @@ if ((NCURR) >= (SIZE)) { \ SIZE += DELTA; \ - PTR = (TYPE *) realloc(PTR,(unsigned)(MAX(((SIZE)*((int)sizeof(TYPE))),1))); \ + PTR = (TYPE *) realloc(PTR,(size_t)(MAX(((SIZE)*((int)sizeof(TYPE))),1))); \ if (PTR == NULL) { \ fprintf(stderr,"failed realloc increment at %d in %s\n", __LINE__, __FILE__); \