Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/parallel.txt =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/parallel.txt (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/parallel.txt (revision 33422) @@ -15,6 +15,77 @@ * relastro_objects, high_speed_catalogs, and standard relastro all call UpdateObjects, one way or another. - The big question is if we can use + The big question is if / how we can use MeasureTiny (or if we need a separate astrometry version) - UpdateObjects: - * if (catalog[i].average[j].Nmeasure == 0) { + * average.d fields we are NOT using in UpdateObjects (28 / 120 bytes) + FIELD Nmissing, NMISSING, unsigned short, number of missings + FIELD Nextend, NEXTEND, unsigned short, number of extended measurements + FIELD missingOffset, OFF_MISSING, int, offset to first missing obs + FIELD extendOffset, OFF_EXTEND, int, offset to first extended measurement + FIELD photFlagsUpper, PHOTFLAGS_U, uint32_t, upper bit of 2 bit summary of per-measure photflags + FIELD photFlagsLower, PHOTFLAGS_L, uint32_t, lower bit of 2 bit summary of per-measure photflags + FIELD extID, EXT_ID, uint64_t, external ID for object (eg PSPS objID) + + * measure.d fields in use in UpdateObjects +* measure[m].t +* measure[m].dbFlags +* measure[m].dt +* measure[0].dR +* measure[0].dD +* measure[0].M +* measure[0].dM +* measure[0].photcode +* measure[0].photFlags +* measure[m].Xccd; +* measure[m].Yccd; +* measure[m].averef; +* measure[0].imageID; +* measure[0].airmass +* measure[0].Mcal (only used in PhotRelTiny) ++ measure[0].dXccd ++ measure[0].dYccd ++ measure[m].dRsys + + float dR; + float dD; + float M; + float Mcal; + float dM; + float airmass; + float Xccd; + float Yccd; + float dt; + int t; + unsigned int averef; + unsigned int imageID; + unsigned int dbFlags; + unsigned int photFlags; +- int catID; // unsigned int? + unsigned short photcode; + +maybe +FIELD Map, M_APER, float, aperture mag, mag +FIELD psfQual, PSF_QF, float, psf coverage/quality factor + +not used +FIELD dMcal, MAG_CAL_ERR, float, systematic calibration error, mag +FIELD az, AZ, float, telescope azimuth +FIELD Sky, SKY_FLUX, float, local estimate of sky flux, counts/sec +FIELD dSky, SKY_FLUX_ERR, float, local estimate of sky flux, counts/sec +FIELD t_msec, TIME_MSEC, unsigned short, time fraction of second, milliseconds +FIELD detID, DET_ID, unsigned int, detection ID +FIELD objID, OBJ_ID, unsigned int, unique ID for object in table +FIELD catID, CAT_ID, unsigned int, unique ID for table in which object was first realized +FIELD psfChisq, PSF_CHISQ, float, psf fit chisq +FIELD psfNdof, PSF_NDOF, int, psf degrees of freedom +FIELD psfNpix, PSF_NPIX, int, psf number of pixels +FIELD crNsigma, CR_NSIGMA, float, Nsigma deviation towards CR +FIELD extNsigma, EXT_NSIGMA, float, Nsigma deviation towards EXT +FIELD FWx, FWHM_MAJOR, short, object fwhm major axis, 1/100 of pixels +FIELD FWy, FWHM_MINOR, short, object fwhm minor axis, 1/100 of pixels +FIELD theta, PSF_THETA, short, angle wrt ccd X dir, (0xffff/360) deg +FIELD Mxx, MXX, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD Mxy, MXY, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD Myy, MYY, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD posangle, POSANGLE, short, position angle sky to chip, (0xffff/360) deg +FIELD pltscale, PLTSCALE, float, plate scale, arcsec/pixel + Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c (revision 33422) @@ -5,4 +5,7 @@ int FixProblemImages (SkyList *skylist) { + + fprintf (stderr, "this function is currently not used : check on change from measure -> measureT in load_catalogs\n"); + abort(); int Nbad; Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c (revision 33422) @@ -181,8 +181,8 @@ off_t idx, ID; - Measure *measure; + MeasureTiny *measure; int i, found; - measure = &catalog[cat].measure[meas]; + measure = &catalog[cat].measureT[meas]; ID = measure[0].imageID; @@ -241,7 +241,7 @@ off_t i; - Measure *measure; - - measure = &catalog[cat].measure[meas]; + MeasureTiny *measure; + + measure = &catalog[cat].measureT[meas]; /* find the image that supplied this measurement */ @@ -333,5 +333,5 @@ } -void dump_measures(Average *average, Measure *measure) { +void dump_measures(Average *average, MeasureTiny *measure) { off_t j, off; @@ -406,7 +406,9 @@ c = clist[im][i]; - X = catalog[c].measure[m].Xccd; - Y = catalog[c].measure[m].Yccd; - n = catalog[c].measure[m].averef; + MeasureTiny *measure = &catalog[c].measureT[m]; + + X = measure[0].Xccd; + Y = measure[0].Yccd; + n = measure[0].averef; if (moscoords == NULL) { @@ -429,20 +431,20 @@ if (fabs(dR) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measure); + dump_measures (&catalog[c].average[n], catalog[c].measureT); // abort (); } if (fabs(dD) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measure); + dump_measures (&catalog[c].average[n], catalog[c].measureT); // abort (); } // complain if the new location is far from the old location - if (fabs(catalog[c].measure[m].dR - dR) > DPOS_MAX_ASEC) { + if (fabs(catalog[c].measureT[m].dR - dR) > DPOS_MAX_ASEC) { fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measure); + dump_measures (&catalog[c].average[n], catalog[c].measureT); // abort(); } - if (fabs(catalog[c].measure[m].dD - dD) > DPOS_MAX_ASEC) { + if (fabs(measure[0].dD - dD) > DPOS_MAX_ASEC) { fprintf (stderr, "measurement is far from original location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); dump_measures (&catalog[c].average[n], catalog[c].measure); @@ -450,23 +452,23 @@ } - dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD); + dPos += SQ(measure[0].dR - dR) + SQ(measure[0].dD - dD); nPos ++; - catalog[c].measure[m].dR = dR; - catalog[c].measure[m].dD = dD; - - if (catalog[c].measure[m].dR > +180.0*3600.0) { + measure[0].dR = dR; + measure[0].dD = dD; + + if (measure[0].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) { + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (measure[0].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); + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); } // set the systematic error for this image: - catalog[c].measure[m].dRsys = ToShortPixels(dPosSys); + measure[0].dRsys = ToShortPixels(dPosSys); } @@ -487,4 +489,7 @@ Mosaic *mosaic; Coords *moscoords, *imcoords, *oldcoords; + + fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n"); + abort(); // check if this image is bad and should be skipped @@ -511,7 +516,10 @@ c = clist[im][i]; - X = catalog[c].measure[m].Xccd; - Y = catalog[c].measure[m].Yccd; - n = catalog[c].measure[m].averef; + // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated + Measure *measure = &catalog[c].measure[m]; + + X = measure[0].Xccd; + Y = measure[0].Yccd; + n = measure[0].averef; dR = dD = 0.0; @@ -527,16 +535,16 @@ } - catalog[c].measure[m].dR = dR; - catalog[c].measure[m].dD = dD; - - if (catalog[c].measure[m].dR > +180.0*3600.0) { + measure[0].dR = dR; + measure[0].dD = dD; + + if (measure[0].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) { + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (measure[0].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); + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); } } @@ -573,13 +581,15 @@ c = clist[im][i]; + MeasureTiny *measure = &catalog[c].measure[m]; + /* apply the current image transformation or use the current value of R+dR, D+dD? */ - raw[i].X = catalog[c].measure[m].Xccd; - raw[i].Y = catalog[c].measure[m].Yccd; - - 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].X = measure[0].Xccd; + raw[i].Y = measure[0].Yccd; + + raw[i].Mag = measure[0].M; + raw[i].dMag = measure[0].dM; + raw[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS); + + n = measure[0].averef; // an object with only one detection provides no information about the image calibration @@ -589,5 +599,5 @@ raw[i].mask |= 0x0001; } - if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) { + if (!finite(measure[0].dR) || !finite(measure[0].dD)) { raw[i].mask |= 0x0002; } @@ -646,11 +656,13 @@ n = catalog[c].measure[m].averef; + MeasureTiny *measure = &catalog[c].measure[m]; + /* apply the current image transformation or use the current value of R+dR, D+dD? */ ref[i].R = catalog[c].average[n].R; ref[i].D = catalog[c].average[n].D; - 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].Mag = measure[0].M; + ref[i].dMag = measure[0].dM; + ref[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS); ref[i].mask = FALSE; @@ -683,5 +695,6 @@ void FlagOutliers (Catalog *catalog) { - // XXX FlagOutliers is just using FlagOutliers2D + + // XXX FlagOutliers is now just using FlagOutliers2D FlagOutliers2D(catalog); return; @@ -716,15 +729,17 @@ for (j = 0; j < catalog[0].Naverage; j++) { + // pointer to this set of measurements + m = catalog[0].average[j].measureOffset; + MeasureTiny *measure = &catalog[0].measureT[m]; + /* accumulate list of valid measurements */ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - R[N] = catalog[0].measure[m].dR; - D[N] = catalog[0].measure[m].dD; - dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); - dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (!MeasFilterTest(&measure[k], FALSE)) continue; + R[N] = measure[k].dR; + D[N] = measure[k].dD; + dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA); + dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC); if (isnan(R[N]) || isnan(D[N])) continue; N++; @@ -741,19 +756,18 @@ /* compare per-object distance to this standard deviation, and flag outliers*/ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { //reset flag on each invocation - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; + if (!MeasFilterTest(&measure[k], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta)))) { - catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; Ndel++; } @@ -809,19 +823,21 @@ for (j = 0; j < catalog[0].Naverage; j++) { + // pointer to this set of measurements + m = catalog[0].average[j].measureOffset; + MeasureTiny *measure = &catalog[0].measureT[m]; + /* accumulate list of valid measurements */ - m = catalog[0].average[j].measureOffset; N = 0; - - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - - //reset flag on each invocation - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { + + // reset flag on each invocation + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - R[N] = catalog[0].measure[m].dR; - D[N] = catalog[0].measure[m].dD; - dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); - dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (!MeasFilterTest(&measure[k], FALSE)) continue; + R[N] = measure[k].dR; + D[N] = measure[k].dD; + dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA); + dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC); if (isnan(R[N]) || isnan(D[N])) continue; N++; @@ -837,15 +853,14 @@ /* calculate deviations of all points*/ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { + // skip bad measurements + if (!MeasFilterTest(&measure[k], FALSE)) continue; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta))); - index[N] = m; + index[N] = k; N++; } @@ -856,10 +871,10 @@ // recalculate image center, sigma based on closest 50% of points - for(k = 0; k < N; k++) { + for (k = 0; k < N; k++) { off_t ind = (off_t) index[k]; - R[k] = catalog[0].measure[ind].dR; - D[k] = catalog[0].measure[ind].dD; - dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA); - dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC); + R[k] = measure[ind].dR; + D[k] = measure[ind].dD; + dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA); + dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC); } liststats (R, dR, N, &statsR); @@ -869,26 +884,25 @@ // use these new statistics to flag outliers - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + if (!MeasFilterTest(&measure[k], FALSE)) continue; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta))); if ((d2[N]) > 1) { - catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; Ndel ++; } N++; Nave++; - } //done rejecting outliers + } // done rejecting outliers // examine results // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); - } //done looping over objects + } // done looping over objects if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); @@ -905,5 +919,5 @@ /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ // we only optionally apply the sigma limit: for object averages, this should not be used (should it?) -int MeasFilterTest(Measure *measure, int applySigmaLim) { +int MeasFilterTest(MeasureTiny *measure, int applySigmaLim) { int found, k; long mask; @@ -958,5 +972,5 @@ /* select measurements by mag limit */ if (ImagSelect) { - mag = PhotInst (measure); + mag = PhotInstTiny (measure); if (mag < ImagMin || mag > ImagMax) return FALSE; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c (revision 33422) @@ -68,9 +68,11 @@ for (i = 0; i < catalog[0].Nmeasure; i++) { - N = getImageByID(catalog[0].measure[i].imageID); + MeasureTiny *measure = &catalog[0].measureT[i] + + N = getImageByID(measure[0].imageID); if (N < 0) continue; - xbin = catalog[0].measure[i].Xccd / starmap[N].Nx; - ybin = catalog[0].measure[i].Yccd / starmap[N].Ny; + xbin = measure[0].Xccd / starmap[N].Nx; + ybin = measure[0].Yccd / starmap[N].Ny; xbin = MAX(0, MIN(NX_MAP-1, xbin)); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33422) @@ -40,4 +40,5 @@ UpdateMeasures (&catalog, 1); + XXX: note that this version of UpdateObjects is / should be expecting catalog.measure, not catalog.measureT UpdateObjects (&catalog, 1); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 33422) @@ -86,13 +86,4 @@ XVERB = FALSE; - // skip objects which are known to be problematic - // XXX include this code or not? -# if (0) - if (catalog[i].average[j].code & STAR_BAD) { - Nskip ++; - continue; - } -# endif - if (catalog[i].average[j].Nmeasure == 0) { continue; @@ -101,40 +92,33 @@ N = 0; m = catalog[i].average[j].measureOffset; - Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25); + MeasureTiny *measure = &catalog[i].measureT[m]; + + Tmin = Tmax = (measure[0].t - T2000) / (86400*365.25); mode = FIT_MODE; // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) Tmean = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - - //does the measurement pass the supplied filtering constraints? - if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) { - catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { + + // does the measurement pass the supplied filtering constraints? + if (!MeasFilterTest(&measure[k], FALSE)) { + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; continue; } //outlier rejection - if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) { - catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; continue; } - // exclude measurements by previous outlier detection - // XXX include this code or not? -# if (0) - if (catalog[i].measure[m].dbFlags & MEAS_BAD) { - catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; - continue; - } -# endif - - catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ; - - R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + measure[k].dbFlags |= ID_MEAS_USED_OBJ; + + R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects // to get Time in years relative to UNIX Tzero (1970/01/01) - T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years + T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years Tmin = MIN(Tmin, T[N]); @@ -143,11 +127,11 @@ // dX, dY : error in arcsec -- - dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); - dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); + dX[N] = GetAstromError (&measure[k], ERROR_MODE_RA); + dY[N] = GetAstromError (&measure[k], ERROR_MODE_DEC); // add systematic error in quadrature, if desired // only do this after the fit has converged (or you will never improve the poor images) // if (INCLUDE_SYS_ERR) { - // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys); + // float dRsys = FromShortPixels(measure[k].dRsys); // dX[N] = hypot(dX[N], dRsys); // dY[N] = hypot(dY[N], dRsys); @@ -157,5 +141,5 @@ // dY[N] = 0.1; - dT[N] = catalog[i].measure[m].dt; + dT[N] = measure[k].dt; // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in @@ -204,5 +188,4 @@ XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC); XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST); - // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY); // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq @@ -311,7 +294,8 @@ // the measure fields must be updated before the average fields m = catalog[i].average[j].measureOffset; + MeasureTiny *measure = &catalog[i].measureT[m]; for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c (revision 33421) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c (revision 33422) @@ -12,9 +12,10 @@ /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ + memset(subcatalog, 0, sizeof(Catalog)); NAVERAGE = 50; NMEASURE = 1000; - ALLOCATE (subcatalog[0].average, Average, NAVERAGE); - ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); - ALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + ALLOCATE (subcatalog[0].average, Average, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); Nmeasure = Naverage = 0; @@ -70,13 +71,12 @@ // allowed. - // CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); - - subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; - subcatalog[0].measure[Nmeasure].averef = Naverage; + CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } @@ -85,5 +85,5 @@ if (Nmeasure == NMEASURE) { NMEASURE += 1000; - REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); } } @@ -96,7 +96,7 @@ } } - REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); - REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1)); - REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt); + REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt); subcatalog[0].Naverage = Naverage; subcatalog[0].Nmeasure = Nmeasure; @@ -168,7 +168,7 @@ // allocate the output data - ALLOCATE (tmpcatalog.average, Average, Nmax); - ALLOCATE (tmpcatalog.measure, Measure, NMEASURE); - ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); off_t Nmeasure = 0; @@ -181,7 +181,10 @@ for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { off_t offset = subcatalog[0].average[ave].measureOffset + j; - tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset]; - tmpcatalog.measure[Nmeasure].averef = i; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; } } @@ -193,9 +196,9 @@ free (subcatalog[0].average); - free (subcatalog[0].measure); + free (subcatalog[0].measureT); free (subcatalog[0].secfilt); subcatalog[0].average = tmpcatalog.average; - subcatalog[0].measure = tmpcatalog.measure; + subcatalog[0].measureT = tmpcatalog.measureT; subcatalog[0].secfilt = tmpcatalog.secfilt; subcatalog[0].Naverage = Nmax;