Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 12730) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 12731) @@ -95,5 +95,7 @@ int STAR_BAD; int MEAS_BAD; -int STAR_TOOFEW; +int POS_TOOFEW; +int PM_TOOFEW; +double PM_DT_MIN; int IMAGE_TOOFEW; double IMAGE_GOOD_FRACTION; @@ -272,6 +274,6 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); int ParFactor (double *pR, double *pD, double R, double D, time_t T); -int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts); -int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts); +int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts); +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts); Image *getMosaicForImage (int N); Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 12731) @@ -25,5 +25,7 @@ GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ); - GetConfig (config, "STAR_TOOFEW", "%d", 0, &STAR_TOOFEW); + GetConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN); + GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); + GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); GetConfig (config, "IMAGE_TOOFEW", "%d", 0, &IMAGE_TOOFEW); GetConfig (config, "IMAGE_GOOD_FRACTION", "%lf", 0, &IMAGE_GOOD_FRACTION); @@ -31,7 +33,7 @@ GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); GetConfig (config, "CATDIR", "%s", 0, CATDIR); - ScanConfig (config, "CATMODE", "%s", 0, CATMODE); - ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); - ScanConfig (config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); + ScanConfig(config, "CATMODE", "%s", 0, CATMODE); + ScanConfig(config, "CATFORMAT", "%s", 0, CATFORMAT); + ScanConfig(config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); sprintf (ImageCat, "%s/Images.dat", CATDIR); @@ -44,10 +46,9 @@ } - GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); - - GetConfig (config, "RELPHOT_GRID_X", "%d", 0, &RELPHOT_GRID_X); - GetConfig (config, "RELPHOT_GRID_Y", "%d", 0, &RELPHOT_GRID_Y); - GetConfig (config, "RELPHOT_GRID_BINNING", "%d", 0, &RELPHOT_GRID_BINNING); - GetConfig (config, "CAMERA_CONFIG", "%s", 0, CameraConfig); + // GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + // GetConfig (config, "RELPHOT_GRID_X", "%d", 0, &RELPHOT_GRID_X); + // GetConfig (config, "RELPHOT_GRID_Y", "%d", 0, &RELPHOT_GRID_Y); + // GetConfig (config, "RELPHOT_GRID_BINNING", "%d", 0, &RELPHOT_GRID_BINNING); + // GetConfig (config, "CAMERA_CONFIG", "%s", 0, CameraConfig); if (*CATMODE == 0) strcpy (CATMODE, "RAW"); Index: /trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPM.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/FitPM.c (revision 12731) @@ -2,5 +2,5 @@ /* do we want an init function which does the alloc and a clear function to free? */ -int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts) { +int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) { int i; @@ -8,4 +8,5 @@ double **A, **B; double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + double chisq, Xf, Yf; /* do I need to do this as 2 2x2 matrix equations? */ @@ -69,4 +70,16 @@ array_free (B, 4); + // add up the chi square for the fit + chisq = 0.0; + for (i = 0; i < Npts; i++) { + Xf = fit[0].Ro + fit[0].uR*T[i]; + Yf = fit[0].Do + fit[0].uD*T[i]; + chisq += SQ(X[i] - Xf) / SQ(dX[i]); + chisq += SQ(Y[i] - Yf) / SQ(dY[i]); + } + fit[0].Nfit = Npts; + + // the reduced chisq is divided by (Ndof = 2*Npts - 4) + fit[0].chisq = chisq / (2.0*Npts - 4.0); return (TRUE); } Index: /trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 12731) @@ -2,5 +2,5 @@ /* do we want an init function which does the alloc and a clear function to free? */ -int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts) { +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) { int i; Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 12731) @@ -8,5 +8,5 @@ static double *pX; static double *pY; -static time_t *T; +static double *T; static double *dT; @@ -24,5 +24,5 @@ ALLOCATE (R, double, MAX (1, Nmax)); ALLOCATE (D, double, MAX (1, Nmax)); - ALLOCATE (T, time_t, MAX (1, Nmax)); + ALLOCATE (T, double, MAX (1, Nmax)); ALLOCATE (X, double, MAX (1, Nmax)); ALLOCATE (Y, double, MAX (1, Nmax)); @@ -44,7 +44,11 @@ Coords coords; PMFit fit; + time_t To; + int mode, Nave, Npm, Npar, Nskip; + double Tmin, Tmax; initObjectData (catalog, Ncatalog); + /* project coordinates to a plane centered on the object with units of arcsec */ coords.crval1 = 0; coords.crval2 = 0; @@ -57,12 +61,27 @@ strcpy (coords.ctype, "RA---SIN"); + // use J2000 as a reference time + To = date_to_sec ("2000/01/01"); + Nave = Npar = Npm = 0; + for (i = 0; i < Ncatalog; i++) { + + if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); + + Nskip = 0; for (j = 0; j < catalog[i].Naverage; j++) { /* calculate the average value of R,D for a single star */ - if (catalog[i].average[j].code & STAR_BAD) continue; + if (catalog[i].average[j].code & STAR_BAD) { + Nskip ++; + continue; + } N = 0; m = catalog[i].average[j].offset; + + Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25); + mode = FIT_MODE; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { if (catalog[i].measure[m].flags & MEAS_BAD) continue; @@ -70,5 +89,8 @@ R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); - T[N] = catalog[i].measure[m].t; + T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years + + Tmin = MIN(Tmin, T[N]); + Tmax = MAX(Tmax, T[N]); /* the astrometric errors are not being carried yet (but should be!) */ @@ -81,6 +103,9 @@ } + if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE; + if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE; + // XXX This criterion needs to be better considered: adjust to match Ndof - if (N < STAR_TOOFEW) { /* too few measurements */ + if (N < POS_TOOFEW) { /* too few measurements */ catalog[i].average[j].code |= ID_STAR_FEW; continue; @@ -95,9 +120,12 @@ /* project all of the R,D coordinates to a plane centered on this coordinate */ for (k = 0; k < N; k++) { - RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + dX[k] = dR[k]; + dY[k] = dD[k]; + // fprintf (stderr, "%d %f %f %f %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]); } /* fit the model components as needed */ - switch (FIT_MODE) { + switch (mode) { case FIT_AVERAGE: liststats (R, dR, N, &statsR); @@ -105,15 +133,22 @@ fit.Ro = statsR.mean; - fit.dRo = statsR.sigma; + fit.dRo = 3600.0*statsR.sigma; fit.Do = statsD.mean; - fit.dDo = statsD.sigma; + fit.dDo = 3600.0*statsD.sigma; fit.chisq = 0.5*(statsR.chisq + statsD.chisq); fit.Nfit = N; + Nave ++; break; case FIT_PM_ONLY: FitPM (&fit, X, dX, Y, dY, T, N); + // fprintf (stderr, "fitted: %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); + // project Ro, Do back to RA,DEC + XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords); + // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); + // continue; + Npm ++; break; @@ -123,4 +158,5 @@ } FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N); + Npar ++; break; default: @@ -130,10 +166,10 @@ if (0 && (j < 100)) { - fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", - catalog[i].average[j].R, - catalog[i].average[j].D, - fit.Ro, fit.Do, - 3600*(catalog[i].average[j].R - fit.Ro), - 3600*(catalog[i].average[j].D - fit.Do)); + fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", + catalog[i].average[j].R, + catalog[i].average[j].D, + fit.Ro, fit.Do, + 3600*(catalog[i].average[j].R - fit.Ro), + 3600*(catalog[i].average[j].D - fit.Do)); } @@ -146,21 +182,24 @@ } - catalog[i].average[j].R = fit.Ro; - catalog[i].average[j].D = fit.Do; - catalog[i].average[j].dR = fit.dRo; - catalog[i].average[j].dD = fit.dDo; - - catalog[i].average[j].uR = fit.uR; - catalog[i].average[j].uD = fit.uD; - catalog[i].average[j].duR = fit.duR; - catalog[i].average[j].duD = fit.duD; - - catalog[i].average[j].P = fit.p; - catalog[i].average[j].dP = fit.dp; + catalog[i].average[j].R = fit.Ro; // RA in degrees + catalog[i].average[j].D = fit.Do; // DEC in degrees + catalog[i].average[j].dR = fit.dRo; // RA scatter in arcsec + catalog[i].average[j].dD = fit.dDo; // DEC scatter in arcsec + + catalog[i].average[j].uR = fit.uR; // RA proper motion in arcsec/year + catalog[i].average[j].uD = fit.uD; // DEC proper motion in arcsec/year + catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year + catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year + + catalog[i].average[j].P = fit.p; // parallax in arcsec + catalog[i].average[j].dP = fit.dp; // parallax error in arcsec catalog[i].average[j].Xp = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG; - } + + if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip); } + + if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar); return (TRUE); } @@ -168,25 +207,25 @@ /* fitting proper-motion and parallax: - given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) - - let x,y be the coordinate in the linearized frame with y parallel to DEC lines - - L,B are the ecliptic longitude and latitude of the object, - dL and dB are the offsets in the L and B directions - - dL = sin(t - topp) - dB = cos(t - topp)*sin(B) - - these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion - for the source in the x,y frame is: - - x = Ro + uR * (t - to) + p * pR - y = Do + uD * (t - to) + p * pD - - the unknowns in these equations are Ro, uR, Do, uD, and p - - XXX think through the concepts for the pole a bit better. all objects near the pole - move the same way with the same phase. choose a projection center and define dL,dB relative - to that center point coordinate system? +given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) + +let x,y be the coordinate in the linearized frame with y parallel to DEC lines + +L,B are the ecliptic longitude and latitude of the object, +dL and dB are the offsets in the L and B directions + +dL = sin(t - topp) +dB = cos(t - topp)*sin(B) + +these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion +for the source in the x,y frame is: + +x = Ro + uR * (t - to) + p * pR +y = Do + uD * (t - to) + p * pD + +the unknowns in these equations are Ro, uR, Do, uD, and p + +XXX think through the concepts for the pole a bit better. all objects near the pole +move the same way with the same phase. choose a projection center and define dL,dB relative +to that center point coordinate system? */ Index: /trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- /trunk/Ohana/src/relastro/src/initialize.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/initialize.c (revision 12731) @@ -48,9 +48,4 @@ } fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF); - fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", - RELPHOT_GRID_X, - RELPHOT_GRID_Y, - RELPHOT_GRID_BINNING, - (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING)); fprintf (stderr, "MAG_LIM %lf\n", MAG_LIM); Index: /trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 12731) @@ -3,5 +3,5 @@ Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) { - int i, Nstar; + int i, j, k, m, Nstar; Catalog *catalog, *pcatalog, tcatalog; @@ -20,4 +20,5 @@ pcatalog[0].catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data pcatalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + pcatalog[0].Nsecfilt = GetPhotcodeNsecfilt (); if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) { @@ -32,4 +33,14 @@ dvo_catalog_unlock (&tcatalog); dvo_catalog_free (&tcatalog); + } else { + if (RESET) { + for (j = 0; j < catalog[i].Naverage; j++) { + catalog[i].average[j].code = 0; + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++) { + catalog[i].measure[m+k].flags = 0; + } + } + } } } Index: /trunk/Ohana/src/relastro/src/save_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/save_catalogs.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/save_catalogs.c (revision 12731) @@ -8,4 +8,5 @@ for (i = 0; i < Ncatalog; i++) { + if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename); dvo_catalog_save (&catalog[i], VERBOSE); dvo_catalog_unlock (&catalog[i]); Index: /trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- /trunk/Ohana/src/relastro/src/select_images.c (revision 12730) +++ /trunk/Ohana/src/relastro/src/select_images.c (revision 12731) @@ -66,4 +66,5 @@ if (VERBOSE) fprintf (stderr, "finding images\n"); + BuildChipMatch (timage, Ntimage); nimage = 0; @@ -90,4 +91,6 @@ } + if (!FindMosaicForImage (timage, Ntimage, i)) continue; + /* define image corners */ Xi[0] = 0; Yi[0] = 0;