Index: /branches/eam_branch_20080223/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/Makefile (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/Makefile (revision 16798) @@ -37,4 +37,6 @@ $(SRC)/UpdateObjects.$(ARCH).o \ $(SRC)/UpdateSimple.$(ARCH).o \ +$(SRC)/UpdateMeasures.$(ARCH).o \ +$(SRC)/GetAstromError.$(ARCH).o \ $(SRC)/args.$(ARCH).o \ $(SRC)/bcatalog.$(ARCH).o \ Index: /branches/eam_branch_20080223/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/include/relastro.h (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/include/relastro.h (revision 16798) @@ -10,4 +10,6 @@ } CoordMode; +typedef enum {ERROR_MODE_RA, ERROR_MODE_DEC, ERROR_MODE_POS} ErrorMode; + typedef struct { double R, D; /* Sky Coords - degrees */ @@ -15,5 +17,7 @@ double L, M; /* Focal Plane - pixels */ double X, Y; /* Chip Coords - pixels */ - double Mag, dMag; + double Mag; + double dMag; + double dPos; int mask; } StarData; @@ -263,9 +267,10 @@ void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt); void fit_eval (CoordFit *fit); +void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1); double **poly2d_dx (double **poly, int Nx, int Ny); double **poly2d_dy (double **poly, int Nx, int Ny); double **poly2d_copy (double **poly, int Nx, int Ny); double poly2d_eval (double **poly, int Nx, int Ny, double x, double y); -void fit_apply_coords (CoordFit *fit, Coords *coords); +CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords); int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo); CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo); @@ -278,4 +283,7 @@ int UpdateChips (Catalog *catalog, int Ncatalog); int UpdateMosaic (Catalog *catalog, int Ncatalog); +int UpdateMeasures (Catalog *catalog, int Ncatalog); +void fixImageRaw (Catalog *catalog, int Ncatalog, int im); + int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); int ParFactor (double *pR, double *pD, double R, double D, time_t T); Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/FitMosaic.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/FitMosaic.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/FitMosaic.c (revision 16798) @@ -5,9 +5,23 @@ int i; CoordFit *fit; + double dP, dQ, dR; fit = fit_init (coords[0].Npolyterms); for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; + + // require radius of XXX arcsec + dP = raw[i].P - ref[i].P; + dQ = raw[i].Q - ref[i].Q; + dR = 3600.0 * hypot (dP, dQ); + + // XXX the value needs to be set in a more intelligent way + if (dR > 0.15) continue; + fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0); + } + if (fit[0].Npts == 0) { + fit_free (fit); + return; } fit_eval (fit); Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/FitSimple.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/FitSimple.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/FitSimple.c (revision 16798) @@ -5,9 +5,23 @@ int i; CoordFit *fit; + double dP, dQ, dR; fit = fit_init (coords[0].Npolyterms); for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; + + // require radius of XXX arcsec + dP = raw[i].P - ref[i].P; + dQ = raw[i].Q - ref[i].Q; + dR = 3600.0 * hypot (dP, dQ); + + // XXX the value needs to be set in a more intelligent way + if (dR > 0.15) continue; + fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0); + } + if (fit[0].Npts == 0) { + fit_free (fit); + return; } fit_eval (fit); @@ -35,2 +49,4 @@ /* XXX I'm not using the errors at all : this could at least be done with the dMag values */ + +/* XXX See notes in FitChips.c */ Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/GetAstromError.c (revision 16798) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/GetAstromError.c (revision 16798) @@ -0,0 +1,36 @@ +# include "relastro.h" + +float GetAstromError (Measure *measure, int mode) { + + PhotCode *code; + float dPobs, dPsys, dPtotal, dM, AS, MS; + + switch (mode) { + case ERROR_MODE_RA: + dPobs = measure[0].dXccd; // need to redefine this as RAerr + break; + case ERROR_MODE_DEC: + dPobs = measure[0].dYccd; // need to redefine this as RAerr + break; + case ERROR_MODE_POS: + dPobs = hypot (measure[0].dXccd, measure[0].dYccd); // need to redefine this as RAerr + break; + default: + abort(); + } + + /* the astrometric errors are not being carried yet (but should be!) */ + /* we use the photometric mag error as a weighting term */ + + code = GetPhotcodebyCode (measure[0].photcode); + AS = code[0].astromErrScale; + MS = code[0].astromErrMagScale; + dPsys = code[0].astromErrSys; + dM = measure[0].dM; + + dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM)); + dPtotal = MAX (dPtotal, MIN_ERROR); + + return (dPtotal); +} + Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/ImageOps.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/ImageOps.c (revision 16798) @@ -33,6 +33,6 @@ for (i = 0; i < Nimage; i++) { - start[i] = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1); - stop[i] = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1); + start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); + stop[i] = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); } } @@ -81,4 +81,5 @@ int i, j; + char *name; for (i = 0; i < Ncatalog; i++) { @@ -87,4 +88,10 @@ } } + + for (i = 0; VERBOSE && (i < Nimage); i++) { + name = GetPhotcodeNamebyCode (image[i].photcode); + fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], + ohana_sec_to_date(image[i].tzero), name); + } } @@ -182,7 +189,75 @@ // return StarData values for detections in the specified image, converting coordinates from the // chip positions: X,Y -> L,M -> P,Q -> R,D +void fixImageRaw (Catalog *catalog, int Ncatalog, int im) { + + int i, m, c, n; + double X, Y, L, M, P, Q, R, D, dR, dD; + + Mosaic *mosaic; + Coords *moscoords, *imcoords; + + moscoords = NULL; + mosaic = getMosaicForImage (im); + if (mosaic != NULL) { + moscoords = &mosaic[0].coords; + } + imcoords = &image[im].coords; + + for (i = 0; i < Nlist[im]; i++) { + m = mlist[im][i]; + c = clist[im][i]; + + X = catalog[c].measure[m].Xccd; + Y = catalog[c].measure[m].Yccd; + n = catalog[c].measure[m].averef; + + if (moscoords == NULL) { + // this is a Simple image (not a mosaic) + // note that for a Simple image, L,M = P,Q + XY_to_LM (&L, &M, X, Y, imcoords); + LM_to_RD (&R, &D, L, M, imcoords); + } else { + XY_to_LM (&L, &M, X, Y, imcoords); + XY_to_LM (&P, &Q, L, M, moscoords); + LM_to_RD (&R, &D, P, Q, moscoords); + } + + // new dR, dD : test + dR = 3600.0*(catalog[c].average[n].R - R); + dD = 3600.0*(catalog[c].average[n].D - D); + + if (fabs(catalog[c].measure[m].dR - dR) > 10.0) { + // XXXXX running into this still for last megacam exposure: wrong mosaic? + // ???? inconsistently hitting this???? + fprintf (stderr, "!"); + abort (); + } + if (fabs(catalog[c].measure[m].dD - dD) > 10.0) { + fprintf (stderr, "*"); + abort (); + } + + catalog[c].measure[m].dR = dR; + catalog[c].measure[m].dD = dD; + + if (catalog[c].measure[m].dR > +180.0*3600.0) { + // average on high end of boundary, move star up + R += 360.0; + catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (catalog[c].measure[m].dR < -180.0*3600.0) { + // average on low end of boundary, move star down + R -= 360.0; + catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + } + } + return; +} + +// return StarData values for detections in the specified image, converting coordinates from the +// chip positions: X,Y -> L,M -> P,Q -> R,D StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) { - int i, m, c; + int i, m, c, n; Mosaic *mosaic; @@ -213,6 +288,12 @@ raw[i].Mag = catalog[c].measure[m].M; raw[i].dMag = catalog[c].measure[m].dM; + raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS); + + n = catalog[c].measure[m].averef; raw[i].mask = FALSE; + if (catalog[c].average[n].Nmeasure < 2) { + raw[i].mask = TRUE; + } switch (mode) { @@ -273,4 +354,5 @@ ref[i].Mag = catalog[c].measure[m].M; ref[i].dMag = catalog[c].measure[m].dM; + ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS); ref[i].mask = FALSE; @@ -286,5 +368,5 @@ case MODE_MOSAIC: RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords); - LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords); + LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords); LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); break; Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateMeasures.c (revision 16798) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateMeasures.c (revision 16798) @@ -0,0 +1,21 @@ +# include "relastro.h" + +int UpdateMeasures (Catalog *catalog, int Ncatalog) { + + int i, Nimage; + Image *image; + + image = getimages (&Nimage); + + for (i = 0; i < Nimage; i++) { + + /* skip DIS images (since they have no associated detections) */ + if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + + /* convert measure coordinates to raw entries */ + fixImageRaw (catalog, Ncatalog, i); + } + + return (TRUE); +} + Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateObjects.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/UpdateObjects.c (revision 16798) @@ -37,37 +37,4 @@ ALLOCATE (pY, double, MAX (1, Nmax)); } - -enum {ERROR_MODE_RA, ERROR_MODE_DEC}; - -float GetAstromError (Measure *measure, int mode) { - - PhotCode *code; - float dPobs, dPsys, dPtotal, dM, AS, MS; - - switch (mode) { - case ERROR_MODE_RA: - dPobs = measure[0].dXccd; // need to redefine this as RAerr - break; - case ERROR_MODE_DEC: - dPobs = measure[0].dYccd; // need to redefine this as RAerr - break; - default: - abort(); - } - - /* the astrometric errors are not being carried yet (but should be!) */ - /* we use the photometric mag error as a weighting term */ - - code = GetPhotcodebyCode (measure[0].photcode); - AS = code[0].astromErrScale; - MS = code[0].astromErrMagScale; - dPsys = code[0].astromErrSys; - dM = measure[0].dM; - - dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM)); - dPtotal = MAX (dPtotal, MIN_ERROR); - - return (dPtotal); -} int UpdateObjects (Catalog *catalog, int Ncatalog) { Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c (revision 16798) @@ -128,4 +128,14 @@ matrix[i][j] = fit[0].sum[ix+jx][iy+jy]; } + + // mask the terms not represented by the Coords terms + if (ix + iy > fit[0].Norder) { + for (j = 0; j < fit[0].Nelems; j++) { + matrix[i][j] = 0.0; + } + vector[i][0] = 0.0; + vector[i][1] = 0.0; + matrix[i][i] = 1.0; + } } @@ -133,6 +143,6 @@ ix = i % fit[0].Nterms; iy = i / fit[0].Nterms; - fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", - ix, iy, vector[i][0], ix, iy, vector[i][1]); + // fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + // ix, iy, vector[i][0], ix, iy, vector[i][1]); } @@ -142,6 +152,6 @@ ix = i % fit[0].Nterms; iy = i / fit[0].Nterms; - fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", - ix, iy, vector[i][0], ix, iy, vector[i][1]); + // fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + // ix, iy, vector[i][0], ix, iy, vector[i][1]); } @@ -158,4 +168,28 @@ } +void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) { + + int ix, iy; + double xterm, yterm, term; + double x, y; + + x = 0.0; + y = 0.0; + + xterm = 1; + for (ix = 0; ix < fit[0].Nterms; ix++) { + yterm = 1; + for (iy = 0; iy < fit[0].Nterms; iy++) { + term = xterm*yterm; + x += fit[0].xfit[ix][iy]*term; + y += fit[0].yfit[ix][iy]*term; + yterm *= y1; + } + xterm *= x1; + } + *x2 = x; + *y2 = y; +} + // Nx, Ny is the number of terms in x and in y double **poly2d_dx (double **poly, int Nx, int Ny) { @@ -209,5 +243,5 @@ for (i = 0; i < Nx; i++) { for (j = 0; j < Ny; j++) { - out[i][j] = poly[i][j+1]; + out[i][j] = poly[i][j]; } } @@ -237,9 +271,9 @@ /* this should only apply to the polynomial, not the projection terms */ /* compare with psastro supporting code */ -void fit_apply_coords (CoordFit *fit, Coords *coords) { +CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) { double c11, c12; double c21, c22; - double Xo, Yo, R; + double Xo, Yo, R1, R2; CoordFit *modfit; @@ -256,48 +290,41 @@ modfit = CoordsSetCenter (fit, Xo, Yo); + /* we do not modify crval1,2: these are kept at the default values */ + + // set cdelt1, cdelt2 + coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]); + coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]); + R1 = 1 / coords[0].cdelt1; + R2 = 1 / coords[0].cdelt2; + // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0) - coords[0].pc1_1 = modfit[0].xfit[1][0]; - coords[0].pc1_2 = modfit[0].xfit[0][1]; - coords[0].pc2_1 = modfit[0].yfit[1][0]; - coords[0].pc2_2 = modfit[0].yfit[0][1]; + coords[0].pc1_1 = modfit[0].xfit[1][0] * R1; + coords[0].pc1_2 = modfit[0].xfit[0][1] * R2; + coords[0].pc2_1 = modfit[0].yfit[1][0] * R1; + coords[0].pc2_2 = modfit[0].yfit[0][1] * R2; // set the polyterm elements if (coords->Npolyterms > 1) { - coords[0].polyterms[0][0] = modfit[0].xfit[2][0]; - coords[0].polyterms[1][0] = modfit[0].xfit[1][1]; - coords[0].polyterms[2][0] = modfit[0].xfit[0][2]; - - coords[0].polyterms[0][1] = modfit[0].yfit[2][0]; - coords[0].polyterms[1][1] = modfit[0].yfit[1][1]; - coords[0].polyterms[2][1] = modfit[0].yfit[0][2]; + coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1; + coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2; + coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2; + + coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1; + coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2; + coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2; } // I need to validate Norder if (coords->Npolyterms > 2) { - coords[0].polyterms[3][0] = modfit[0].xfit[3][0]; - coords[0].polyterms[4][0] = modfit[0].xfit[2][1]; - coords[0].polyterms[5][0] = modfit[0].xfit[1][2]; - coords[0].polyterms[6][0] = modfit[0].xfit[0][3]; - - coords[0].polyterms[3][1] = modfit[0].yfit[3][0]; - coords[0].polyterms[4][1] = modfit[0].yfit[2][1]; - coords[0].polyterms[5][1] = modfit[0].yfit[1][2]; - coords[0].polyterms[6][1] = modfit[0].yfit[0][3]; - } - - /* we do not modify crval1,2: these are kept at the default values */ - - /* normalize pc11,etc */ - c11 = coords[0].pc1_1; - c21 = coords[0].pc1_1; - c12 = coords[0].pc1_1; - c22 = coords[0].pc1_1; - coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21)); - R = 1 / coords[0].cdelt1; - - coords[0].pc1_1 = c11*R; - coords[0].pc2_1 = c21*R; - coords[0].pc1_2 = c12*R; - coords[0].pc2_2 = c22*R; + coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1; + coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2; + coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2; + coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2; + + coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1; + coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2; + coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2; + coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2; + } fit_free (modfit); Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/mkpolyterm.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/mkpolyterm.c (revision 16798) @@ -20,5 +20,5 @@ int i, Nx, Ny; - double R, Xo, Yo, dPos; + double R, Xo, Yo, dPos, dPosRef; double **XdX, **XdY, **YdX, **YdY, **alpha, **beta; double **xfit, **yfit; @@ -51,4 +51,5 @@ */ dPos = tol + 1; + dPosRef = dPos; for (i = 0; (dPos > tol) && (i < 20); i++) { // NOTE: order for alpha is: [y][x] @@ -66,4 +67,10 @@ Yo -= beta[1][0]; dPos = hypot(beta[0][0], beta[1][0]); + if (i == 0) { + dPosRef = dPos; + } + } + if (dPos > dPosRef) { + fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n"); } array_free (alpha, 2); Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/relastro.c =================================================================== --- /branches/eam_branch_20080223/Ohana/src/relastro/src/relastro.c (revision 16797) +++ /branches/eam_branch_20080223/Ohana/src/relastro/src/relastro.c (revision 16798) @@ -61,7 +61,23 @@ if (!UPDATE) exit (0); - save_catalogs (catalog, Ncatalog); + if (FIT_TARGET == TARGET_OBJECTS) { + save_catalogs (catalog, Ncatalog); + } else { + // XXX for now, reload all catalogs at once; as memory use gets large, reload one-at-a-time - if (FIT_TARGET != TARGET_OBJECTS) { + // load catalog data from region files : do not subselect + catalog = load_catalogs (skylist, &Ncatalog, FALSE); + + // match measurements with images + initImageBins (catalog, Ncatalog); + findImages (catalog, Ncatalog); + + // update the detection coordinates using the new image parameters + UpdateMeasures (catalog, Ncatalog); + + // write the updated detections to disk + save_catalogs (catalog, Ncatalog); + + // save the updated image parameters dvo_image_update (&db, VERBOSE); dvo_image_unlock (&db);