Index: /branches/eam_branches/ipp-20140813/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140813/Ohana/src/relastro/Makefile (revision 37305) +++ /branches/eam_branches/ipp-20140813/Ohana/src/relastro/Makefile (revision 37306) @@ -59,4 +59,5 @@ $(SRC)/save_catalogs.$(ARCH).o \ $(SRC)/CoordOps.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ $(SRC)/FixProblemImages.$(ARCH).o \ $(SRC)/StarMaps.$(ARCH).o \ @@ -117,4 +118,5 @@ $(SRC)/save_catalogs.$(ARCH).o \ $(SRC)/CoordOps.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ $(SRC)/high_speed_catalogs.$(ARCH).o \ $(SRC)/high_speed_objects.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140813/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20140813/Ohana/src/relastro/include/relastro.h (revision 37305) +++ /branches/eam_branches/ipp-20140813/Ohana/src/relastro/include/relastro.h (revision 37306) @@ -133,4 +133,12 @@ int Nmeas; } StatType; + +typedef struct { + off_t Nave; + off_t Npm; + off_t Npar; + off_t Nskip; + off_t Noffset; +} FitStats; /* global variables set in parameter file */ @@ -535,2 +543,6 @@ int areImagesLoaded (); int areImagesMatched (); + +int isGPC1chip (int photcode); +int isGPC1stack (int photcode); +int isGPC1warp (int photcode); Index: /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/ParFactor.c =================================================================== --- /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/ParFactor.c (revision 37305) +++ /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/ParFactor.c (revision 37306) @@ -63,5 +63,6 @@ double jd, lambda, beta, epsilon, Radius; - /* given a time T in UNIX seconds, determine the solar longitude S */ + /* given a Time relative to Tmean, Tmean in years since J2000, determine the solar + longitude S */ // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); Index: /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/UpdateObjects.c (revision 37305) +++ /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/UpdateObjects.c (revision 37306) @@ -1,3 +1,7 @@ # include "relastro.h" +# define PAR_TOOFEW 5 + +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m); +int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); static off_t Nmax; @@ -13,4 +17,26 @@ static double *C_red; +static Coords coords; + +static time_t T2000; + +void initFitStats (FitStats *fitStats) { + fitStats->Nave = 0; + fitStats->Npm = 0; + fitStats->Npar = 0; + fitStats->Nskip = 0; + fitStats->Noffset = 0; + return; +} + +void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) { + tgtFitStats->Nave += srcFitStats->Nave ; + tgtFitStats->Npm += srcFitStats->Npm ; + tgtFitStats->Npar += srcFitStats->Npar ; + tgtFitStats->Nskip += srcFitStats->Nskip ; + tgtFitStats->Noffset += srcFitStats->Noffset ; + return; +} + void initObjectData (Catalog *catalog, int Ncatalog) { @@ -41,48 +67,4 @@ ALLOCATE (C_blue, double, MAX (1, Nmax)); ALLOCATE (C_red, double, MAX (1, Nmax)); -} - -void freeObjectData () { - - free (R); - free (D); - free (T); - free (X); - free (Y); - - free (dR); - free (dD); - free (dT); - free (dX); - free (dY); - - free (pX); - free (pY); - - free (C_blue); - free (C_red); -} - -// This function operates on both Measure and MeasureTiny. In the big stages, this should -// be called with just MeasureTiny set and Measure == NULL -int UpdateObjects (Catalog *catalog, int Ncatalog) { - - off_t j, k, m; - int i, N, NcBlue, NcRed, Nsecfilt, mode, result, status, XVERB; - StatType statsR, statsD; - Coords coords; - PMFit fitAve, fitPM, fitPAR, fit; - time_t T2000; - off_t Nave, Npm, Npar, Nskip, Noffset; - off_t NaveSum, NpmSum, NparSum, NskipSum, NoffSum; - double Tmin, Tmax, Tmean, Trange; - - memset (&fit, 0, sizeof(fit)); - memset (&fitAve, 0, sizeof(fitAve)); - memset (&fitPM, 0, sizeof(fitPM)); - memset (&fitPAR, 0, sizeof(fitPAR)); - initObjectData (catalog, Ncatalog); - - int setRefColor = areImagesMatched(); /* project coordinates to a plane centered on the object with units of arcsec */ @@ -97,407 +79,605 @@ strcpy (coords.ctype, "DEC--SIN"); - XVERB = FALSE; - // use J2000 as a reference time - T2000 = ohana_date_to_sec ("2000/01/01"); + T2000 = ohana_date_to_sec ("2000/01/01,12:00:00"); +} + +void freeObjectData () { + + free (R); + free (D); + free (T); + free (X); + free (Y); + + free (dR); + free (dD); + free (dT); + free (dX); + free (dY); + + free (pX); + free (pY); + + free (C_blue); + free (C_red); +} + +// This function operates on both Measure and MeasureTiny. In the big stages, this should +// be called with just MeasureTiny set and Measure == NULL +int UpdateObjects (Catalog *catalog, int Ncatalog) { + + initObjectData (catalog, Ncatalog); + // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? - - Nsecfilt = GetPhotcodeNsecfilt (); + int Nsecfilt = GetPhotcodeNsecfilt (); if (Ncatalog) { assert (catalog[0].Nsecfilt == Nsecfilt); } - NaveSum = NparSum = NpmSum = NoffSum = NskipSum = 0; + FitStats sumStatsChips; initFitStats (&sumStatsChips); + FitStats sumStatsStack; initFitStats (&sumStatsStack); + + int i; for (i = 0; i < Ncatalog; i++) { if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); - Nave = Npar = Npm = Nskip = Noffset = 0; + FitStats fitStatsChips; initFitStats (&fitStatsChips); + FitStats fitStatsStack; initFitStats (&fitStatsStack); + + off_t j; for (j = 0; j < catalog[i].Naverage; j++) { /* calculate the average value of R,D for a single star */ - - XVERB = FALSE; - fitAve.chisq = NAN; - fitPM.chisq = NAN; - fitPAR.chisq = NAN; - - // if we fail to fit the astrometry for some reason, we need to set/reset these - catalog[i].average[j].flags |= ID_STAR_NO_ASTROM; - catalog[i].average[j].ChiSqAve = NAN; - catalog[i].average[j].ChiSqPM = NAN; - catalog[i].average[j].ChiSqPar = NAN; - catalog[i].average[j].Npos = 0; - - if (catalog[i].average[j].Nmeasure == 0) { + off_t m = catalog[i].average[j].measureOffset; + MeasureTiny *measure = &catalog[i].measureT[m]; + Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL; + Average *average = &catalog[i].average[j]; + SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt]; + + UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack); + UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m); + } + if (VERBOSE) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsChips.Nave, fitStatsChips.Npm, fitStatsChips.Npar, fitStatsChips.Nskip, fitStatsChips.Noffset); + if (VERBOSE) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsStack.Nave, fitStatsStack.Npm, fitStatsStack.Npar, fitStatsStack.Nskip, fitStatsStack.Noffset); + sumFitStats (&fitStatsChips, &sumStatsChips); + sumFitStats (&fitStatsStack, &sumStatsStack); + } + freeObjectData (); + + if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar), sumStatsChips.Nave, sumStatsChips.Npm, sumStatsChips.Npar, sumStatsChips.Nskip, sumStatsChips.Noffset); + if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar), sumStatsStack.Nave, sumStatsStack.Npm, sumStatsStack.Npar, sumStatsStack.Nskip, sumStatsStack.Noffset); + return (TRUE); +} + +// This function operates on both Measure and MeasureTiny. In the big stages, this should +// be called with just MeasureTiny set and Measure == NULL +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m) { + + int setRefColor = areImagesMatched(); + + /* calculate the average value of R,D for a single star */ + + PMFit fit; memset (&fit, 0, sizeof(fit)); + PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN; + PMFit fitPM; memset (&fitPM, 0, sizeof(fitPM)); fitPM.chisq = NAN; + PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN; + + // if we fail to fit the astrometry for some reason, we need to set/reset these + average[0].flags |= ID_STAR_NO_ASTROM; + average[0].ChiSqAve = NAN; + average[0].ChiSqPM = NAN; + average[0].ChiSqPar = NAN; + average[0].Npos = 0; + + // an object with no measurements is externally supplied + if (average[0].Nmeasure == 0) return TRUE; + + int NcBlue = 0; + int NcRed = 0; + int N = 0; + + int mode = FIT_MODE; // start with the globally-defined fit mode + + int XVERB = FALSE; + XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC); + XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST); + + // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) + off_t k; + for (k = 0; k < average[0].Nmeasure; k++) { + + if (XVERB) { + char *date = ohana_sec_to_date (measure[k].t); + int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; + fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); + free (date); + } + + // SKIP gpc1 stack data + if (isGPC1stack(measure[k].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(measure[k].photcode)) continue; + + // reset the bit to note that a detection was used (or not) + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } + + // does the measurement pass the supplied filtering constraints? + // MeasFilterTestTiny does not test psfQF + // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst + int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); + if (!keepMeasure) { + continue; + } + + double Ri = getMeanR (&measure[k], average, secfilt); + double Di = getMeanD (&measure[k], average, secfilt); + + // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius + if (ExcludeBogus) { + coords.crval1 = average[0].R; + coords.crval2 = average[0].D; + double Xi, Yi; + RD_to_XY (&Xi, &Yi, Ri, Di, &coords); + double radius = hypot(Xi, Yi); + if (radius > ExcludeBogusRadius) { + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; } continue; } - - NcBlue = 0; - NcRed = 0; - N = 0; - m = catalog[i].average[j].measureOffset; - MeasureTiny *measure = &catalog[i].measureT[m]; - Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL; - // when we update the output measure values, we need to do it here - - mode = FIT_MODE; - - // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e); - 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); - - // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) - for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { - - if (XVERB) { - char *date = ohana_sec_to_date (measure[k].t); - int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; - fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); - free (date); - } - - // reset the bit to note that a detection was used (or not) - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } - - // does the measurement pass the supplied filtering constraints? - // MeasFilterTestTiny does not test psfQF - // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst - int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); - if (!keepMeasure) { - continue; - } - - // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius - if (ExcludeBogus) { - double Ri = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - double Di = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - coords.crval1 = catalog[i].average[j].R; - coords.crval2 = catalog[i].average[j].D; - double Xi, Yi; - RD_to_XY (&Xi, &Yi, Ri, Di, &coords); - double radius = hypot(Xi, Yi); - if (radius > ExcludeBogusRadius) { - measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; - if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; } - continue; - } - measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; } - } - - // outlier rejection - if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { - continue; - } - - 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] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years - - // dX, dY : error in arcsec -- - dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); - dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); - - // allow a given photcode or measurement to be - // ignored if the error is NAN (for photcode, set astromErrSys to NaN) - if (isnan(dX[N])) continue; - if (isnan(dY[N])) continue; - - // 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(measure[k].dRsys); - // dX[N] = hypot(dX[N], dRsys); - // dY[N] = hypot(dY[N], dRsys); - // } - - // dX[N] = 0.1; - // dY[N] = 0.1; - - dT[N] = measure[k].dt; - - // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in - // arcseconds. dR, dD are the errors in those directions in degrees. IF we have - // non-circular errors (different values for X and Y), then dR and dD will be - // incorrect: they would need to be rotated to take out the position angle - dR[N] = dX[N] / 3600.0; - dD[N] = dY[N] / 3600.0; - - if (setRefColor) { - float colorBlue = getColorBlue (m+k, i); - if (!isnan(colorBlue)) { - C_blue[NcBlue] = colorBlue; - NcBlue++; - } - float colorRed = getColorRed (m+k, i); - if (!isnan(colorRed)) { - C_red[NcRed] = colorRed; - NcRed++; - } - } - - measure[k].dbFlags |= ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } - - N++; - } // loop over measurements : catalog[i].average[j].Nmeasure - - // if we have too few good detections for the desired fit, or too limited a - // baseline, use a fit with fewer parameters. XXX if we have too few measurements - // for even the average position, consider including the lower-quality detections? - - // find Tmin & Tmax from the list of accepted measurements - Tmean = 0; - Tmin = Tmax = T[0]; - for (k = 0; k < N; k++) { - Tmin = MIN(Tmin, T[k]); - Tmax = MAX(Tmax, T[k]); - Tmean += T[k]; + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; } + } + + // outlier rejection + if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { + continue; + } + + R[N] = Ri; + D[N] = Di; + + // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000. + // T[] is time in years since J2000 (jd = 2451545) + T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years + + // dX, dY : error in arcsec -- + dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); + dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); + + // allow a given photcode or measurement to be + // ignored if the error is NAN (for photcode, set astromErrSys to NaN) + if (isnan(dX[N])) continue; + if (isnan(dY[N])) continue; + + // 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(measure[k].dRsys); + // dX[N] = hypot(dX[N], dRsys); + // dY[N] = hypot(dY[N], dRsys); + // } + + // dX[N] = 0.1; + // dY[N] = 0.1; + + dT[N] = measure[k].dt; + + // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in + // arcseconds. dR, dD are the errors in those directions in degrees. IF we have + // non-circular errors (different values for X and Y), then dR and dD will be + // incorrect: they would need to be rotated to take out the position angle + dR[N] = dX[N] / 3600.0; + dD[N] = dY[N] / 3600.0; + + if (setRefColor) { + float colorBlue = getColorBlue (m+k, i); + if (!isnan(colorBlue)) { + C_blue[NcBlue] = colorBlue; + NcBlue++; } - // XXX add the parallax factor range as a criterion as well - Trange = Tmax - Tmin; - if (Trange < PM_DT_MIN) mode = FIT_AVERAGE; - if (((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) && (N <= PM_TOOFEW)) mode = FIT_AVERAGE; - - if (RELASTRO_OP == OP_HIGH_SPEED) { - Tmean = 0.5*(Tmax - Tmin); - } else { - Tmean /= (float) N; + float colorRed = getColorRed (m+k, i); + if (!isnan(colorRed)) { + C_red[NcRed] = colorRed; + NcRed++; } - - // too few measurements for average position (require 2 values) - if (N < SRC_MEAS_TOOFEW) { - if (N < 2) continue; + } + + measure[k].dbFlags |= ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } + + N++; + } // loop over measurements : average[0].Nmeasure + + if (N < 1) { + if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) { + average[0].R = average[0].Rstk; + average[0].D = average[0].Dstk; + average[0].dR = average[0].dRstk; + average[0].dD = average[0].dDstk; + } + return FALSE; + } + + // if we have too few good detections for the desired fit, or too limited a + // baseline, use a fit with fewer parameters. XXX if we have too few measurements + // for even the average position, consider including the lower-quality detections? + + // find Tmin & Tmax from the list of accepted measurements + double Tmean = 0.0; + double Tmin = T[0]; + double Tmax = T[0]; + for (k = 0; k < N; k++) { + Tmin = MIN(Tmin, T[k]); + Tmax = MAX(Tmax, T[k]); + Tmean += T[k]; + } + double Trange = Tmax - Tmin; + + if (RELASTRO_OP == OP_HIGH_SPEED) { + Tmean = 0.5*(Tmax - Tmin); + } else { + Tmean /= (float) N; + } + + /* we need to do the fit in a locally linear space; choose a ref coordinate */ + coords.crval1 = R[0]; + coords.crval2 = D[0]; + + // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq + + // *** first fit for the proper motion (skip fit if Trange or Npts is too small) *** + if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { + if (Trange < PM_DT_MIN) { + mode = FIT_AVERAGE; + goto skipPM; + } + if (N <= PM_TOOFEW) { + mode = FIT_AVERAGE; + goto skipPM; + } + + // project all of the R,D coordinates to a plane centered on this coordinate. set + // the times to be relative to Tmean (this is required for parallax as well) + for (k = 0; k < N; k++) { + RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + T[k] -= Tmean; + if (XVERB) { + fprintf (stderr, OFF_T_FMT" %f %f %f %f %f +/- %f %f\n", k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); } - - /* we need to do the fit in a locally linear space; choose a ref coordinate */ - coords.crval1 = R[0]; - coords.crval2 = D[0]; - - // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq - - if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { - // project all of the R,D coordinates to a plane centered on this coordinate. set - // the times to be relative to Tmean (this is required for parallax as well) - for (k = 0; k < N; k++) { - RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); - T[k] -= Tmean; - if (XVERB) { - fprintf (stderr, OFF_T_FMT" %f %f %f %f %f +/- %f %f\n", k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); - } - } - - FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB); - - if (XVERB) fprintf (stderr, "fitted PM: %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq); - - // project Ro, Do back to RA,DEC - XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords); - if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p); - - fitPM.p = fitPM.dp = 0.0; - catalog[i].average[j].flags |= ID_STAR_FIT_PM; - Npm ++; - - // XXX a hard-wired hack... - if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) { - mode = FIT_AVERAGE; - catalog[i].average[j].flags |= ID_STAR_BAD_PM; - } - } - - if (mode == FIT_PM_AND_PAR) { - float pXmin = +2.0; - float pXmax = -2.0; - float pYmin = +2.0; - float pYmax = -2.0; - for (k = 0; k < N; k++) { - ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean); - pXmin = MIN (pXmin, pX[k]); - pXmax = MAX (pXmax, pX[k]); - pYmin = MIN (pYmin, pY[k]); - pYmax = MAX (pYmax, pY[k]); - } - float dXRange = pXmax - pXmin; - float dYRange = pYmax - pYmin; - float parRange = hypot (dXRange, dYRange); + } + + FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB); + + if (XVERB) fprintf (stderr, "fitted PM: %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq); + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords); + if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p); + + fitPM.p = fitPM.dp = 0.0; + average[0].flags |= ID_STAR_FIT_PM; + fitStats->Npm ++; + + // XXX a hard-wired hack... + if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) { + mode = FIT_AVERAGE; + average[0].flags |= ID_STAR_BAD_PM; + } + } + +skipPM: + // fit the parallax + proper-motion model + // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail + // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only + if (mode == FIT_PM_AND_PAR) { + if (Trange < PM_DT_MIN) { + mode = FIT_PM_ONLY; + goto skipPAR; + } + if (N <= PAR_TOOFEW) { + mode = FIT_PM_ONLY; + goto skipPAR; + } + float pXmin = +2.0; + float pXmax = -2.0; + float pYmin = +2.0; + float pYmax = -2.0; + for (k = 0; k < N; k++) { + ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean); + pXmin = MIN (pXmin, pX[k]); + pXmax = MAX (pXmax, pX[k]); + pYmin = MIN (pYmin, pY[k]); + pYmax = MAX (pYmax, pY[k]); + } + float dXRange = pXmax - pXmin; + float dYRange = pYmax - pYmin; + float parRange = hypot (dXRange, dYRange); -# define PAR_TOOFEW 5 - if ((parRange >= PAR_FACTOR_MIN) && (N > PAR_TOOFEW)) { - FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB); - if (XVERB) fprintf (stderr, "fitted PM+PAR: %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq); - - XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); - catalog[i].average[j].flags |= ID_STAR_FIT_PAR; - Npar ++; - - // XXX a hard-wired hack... - if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { - mode = FIT_AVERAGE; - catalog[i].average[j].flags |= ID_STAR_BAD_PM; - } - } else { - // need to set mode = FIT_PM_ONLY if we do not fit for parallax - mode = FIT_PM_ONLY; - } - } - - // fit the average model - if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { - liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) - liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) - - fitAve.Ro = statsR.mean; - fitAve.dRo = 3600.0*statsR.sigma; - - fitAve.Do = statsD.mean; - fitAve.dDo = 3600.0*statsD.sigma; - - fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq); - fitAve.Nfit = N; - - fitAve.uR = fitAve.duR = 0.0; - fitAve.uD = fitAve.duD = 0.0; - fitAve.p = fitAve.dp = 0.0; - catalog[i].average[j].flags |= ID_STAR_FIT_AVE; - Nave ++; - } - - if (setRefColor) { - float colorMedian; - dsort (C_blue, NcBlue); - colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN; - catalog[i].average[j].refColorBlue = colorMedian; - dsort (C_red, NcRed); - colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN; - catalog[i].average[j].refColorRed = colorMedian; - } - - /* choose the result based on the chisq values */ - // XXXX for now, just use the mode as the result: - result = mode; - - switch (result) { - case FIT_AVERAGE: - catalog[i].average[j].flags |= ID_STAR_USE_AVE; - fit = fitAve; - break; - case FIT_PM_ONLY: - catalog[i].average[j].flags |= ID_STAR_USE_PM; - fit = fitPM; - break; - case FIT_PM_AND_PAR: - catalog[i].average[j].flags |= ID_STAR_USE_PAR; - fit = fitPAR; - break; - } - if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%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), - fit.uR, fit.uD, fit.p, fit.dp); - - // make sure that the fit succeeded - status = TRUE; - status &= finite(fit.Ro); - status &= finite(fit.Do); - status &= finite(fit.dRo); - status &= finite(fit.dDo); - status &= finite(fit.uR); - status &= finite(fit.uD); - status &= finite(fit.duR); - status &= finite(fit.duD); - status &= finite(fit.p); - status &= finite(fit.dp); - if (!status) { - Nskip ++; - continue; - } - - // what is the offset relative to the mean fit position? - coords.crval1 = catalog[i].average[j].R; - coords.crval2 = catalog[i].average[j].D; - - double dXoff, dYoff; - RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords); - float dPos = hypot (dXoff, dYoff); - if (dPos > MaxMeanOffset) { - if (Noffset < 100) { - fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff); - } - Noffset ++; - continue; - } - - - // the measure fields must be updated before the average fields - for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { - 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]); - if (measureBig) { - setMeanR_Big (fit.Ro, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - setMeanD_Big (fit.Do, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - } - } - - 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].ChiSqAve = fitAve.chisq; - catalog[i].average[j].ChiSqPM = fitPM.chisq; - catalog[i].average[j].ChiSqPar = fitPAR.chisq; - catalog[i].average[j].Tmean = (Tmean * 86400 * 365.25) + T2000; - catalog[i].average[j].Trange = (Trange * 86400 * 365.25); - catalog[i].average[j].Npos = fit.Nfit; - - // XXX EAM 20140812: for a test, set average.Rstk,Dstk, etc to match R,D - catalog[i].average[j].Rstk = fit.Ro; // RA in degrees - catalog[i].average[j].Dstk = fit.Do; // DEC in degrees - catalog[i].average[j].dRstk = fit.dRo; // RA scatter in arcsec - catalog[i].average[j].dDstk = fit.dDo; // DEC scatter in arcsec - - // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) - catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM; - - if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%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), - catalog[i].average[j].uR, - catalog[i].average[j].uD, - fitAve.chisq, fitPM.chisq, fitPAR.chisq); - } - - NaveSum += Nave; - NpmSum += Npm; - NparSum += Npar; - NskipSum += Nskip; - NoffSum += Noffset; - if (VERBOSE) fprintf (stderr, "catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, Nave, Npm, Npar, Nskip, Noffset); - } - - freeObjectData (); - - if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (NaveSum + NpmSum + NparSum), NaveSum, NpmSum, NparSum, NskipSum, NoffSum); + if (parRange < PAR_FACTOR_MIN) { + mode = FIT_PM_ONLY; + goto skipPAR; + } + + FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB); + if (XVERB) fprintf (stderr, "fitted PM+PAR: %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq); + + XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); + average[0].flags |= ID_STAR_FIT_PAR; + fitStats->Npar ++; + + // XXX a hard-wired hack... + if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { + mode = FIT_PM_ONLY; + } + } + +skipPAR: + { + // ALWAYS fit the average model + StatType statsR, statsD; + liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) + liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) + + fitAve.Ro = statsR.mean; + fitAve.dRo = 3600.0*statsR.sigma; + + fitAve.Do = statsD.mean; + fitAve.dDo = 3600.0*statsD.sigma; + + fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN; + fitAve.Nfit = N; + + fitAve.uR = fitAve.duR = 0.0; + fitAve.uD = fitAve.duD = 0.0; + fitAve.p = fitAve.dp = 0.0; + average[0].flags |= ID_STAR_FIT_AVE; + fitStats->Nave ++; + } + + if (setRefColor) { + float colorMedian; + dsort (C_blue, NcBlue); + colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN; + average[0].refColorBlue = colorMedian; + dsort (C_red, NcRed); + colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN; + average[0].refColorRed = colorMedian; + } + + /* choose the result based on the chisq values */ + // XXXX for now, just use the mode as the result: + int result = mode; + + switch (result) { + case FIT_AVERAGE: + average[0].flags |= ID_STAR_USE_AVE; + fit = fitAve; + break; + case FIT_PM_ONLY: + average[0].flags |= ID_STAR_USE_PM; + fit = fitPM; + break; + case FIT_PM_AND_PAR: + average[0].flags |= ID_STAR_USE_PAR; + fit = fitPAR; + break; + } + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n", + average[0].R, + average[0].D, + fit.Ro, fit.Do, + 3600*(average[0].R - fit.Ro), + 3600*(average[0].D - fit.Do), + fit.uR, fit.uD, fit.p, fit.dp); + + // make sure that the fit succeeded + int status = TRUE; + status &= finite(fit.Ro); + status &= finite(fit.Do); + status &= finite(fit.dRo); + status &= finite(fit.dDo); + status &= finite(fit.uR); + status &= finite(fit.uD); + status &= finite(fit.duR); + status &= finite(fit.duD); + status &= finite(fit.p); + status &= finite(fit.dp); + if (!status) { + fitStats->Nskip ++; + return FALSE; // XXX ?? + } + + // what is the offset relative to the mean fit position? + coords.crval1 = average[0].R; + coords.crval2 = average[0].D; + + double dXoff, dYoff; + RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords); + float dPos = hypot (dXoff, dYoff); + if (dPos > MaxMeanOffset) { + if (fitStats->Noffset < 100) { + fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff); + } + fitStats->Noffset ++; + return FALSE; // XXX ?? + } + + average[0].R = fit.Ro; // RA in degrees + average[0].D = fit.Do; // DEC in degrees + average[0].dR = fit.dRo; // RA scatter in arcsec + average[0].dD = fit.dDo; // DEC scatter in arcsec + + average[0].uR = fit.uR; // RA proper motion in arcsec/year + average[0].uD = fit.uD; // DEC proper motion in arcsec/year + average[0].duR = fit.duR; // RA proper motion error in arcsec/year + average[0].duD = fit.duD; // DEC proper motion error in arcsec/year + + average[0].P = fit.p; // parallax in arcsec + average[0].dP = fit.dp; // parallax error in arcsec + + average[0].ChiSqAve = fitAve.chisq; + average[0].ChiSqPM = fitPM.chisq; + average[0].ChiSqPar = fitPAR.chisq; + + average[0].Tmean = (Tmean * 86400 * 365.25) + T2000; + average[0].Trange = (Trange * 86400 * 365.25); + average[0].Npos = fit.Nfit; + + // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) + average[0].flags &= ~ID_STAR_NO_ASTROM; + + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n", + average[0].R, + average[0].D, + fit.Ro, fit.Do, + 3600*(average[0].R - fit.Ro), + 3600*(average[0].D - fit.Do), + average[0].uR, + average[0].uD, + fitAve.chisq, fitPM.chisq, fitPAR.chisq); + return (TRUE); } + +// This function operates on both Measure and MeasureTiny. In the big stages, this should +// be called with just MeasureTiny set and Measure == NULL +int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) { + + off_t k; + + // set the default values + average[0].Rstk = NAN; // RA in degrees + average[0].Dstk = NAN; // DEC in degrees + average[0].dRstk = NAN; // RA scatter in arcsec + average[0].dDstk = NAN; // DEC scatter in arcsec + + /* calculate the average value of R,D for a single star */ + PMFit fitAve; + memset (&fitAve, 0, sizeof(fitAve)); + fitAve.chisq = NAN; + + if (average[0].Nmeasure == 0) return TRUE; + + int N = 0; + + int XVERB = FALSE; + XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC); + XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST); + + // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) + for (k = 0; k < average[0].Nmeasure; k++) { + + if (XVERB) { + char *date = ohana_sec_to_date (measure[k].t); + int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0; + fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n", k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig); + free (date); + } + + // SKIP everything except gpc1 stack data + if (!isGPC1stack(measure[k].photcode)) continue; + + // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst + int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); + if (!keepMeasure) { + continue; + } + + R[N] = getMeanR (&measure[k], average, secfilt); + D[N] = getMeanD (&measure[k], average, secfilt); + + // dX, dY : error in arcsec -- + dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); + dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); + + // allow a given photcode or measurement to be + // ignored if the error is NAN (for photcode, set astromErrSys to NaN) + if (isnan(dX[N])) continue; + if (isnan(dY[N])) continue; + + // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in + // arcseconds. dR, dD are the errors in those directions in degrees. IF we have + // non-circular errors (different values for X and Y), then dR and dD will be + // incorrect: they would need to be rotated to take out the position angle + dR[N] = dX[N] / 3600.0; + dD[N] = dY[N] / 3600.0; + + // XXX use a different flag for stack measurements? + // measure[k].dbFlags |= ID_MEAS_USED_OBJ; + // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } + + N++; + } // loop over measurements : average[0].Nmeasure + + // if we have too few good detections for the desired fit, or too limited a + // baseline, use a fit with fewer parameters. XXX if we have too few measurements + // for even the average position, consider including the lower-quality detections? + + // too few measurements for average position (require 2 values) + if (N < 1) return FALSE; // XXX ?? + + // find the mean position + StatType statsR, statsD; + liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) + liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) + + fitAve.Ro = statsR.mean; + fitAve.dRo = 3600.0*statsR.sigma; + + fitAve.Do = statsD.mean; + fitAve.dDo = 3600.0*statsD.sigma; + + fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq); + fitAve.Nfit = N; + + // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE; + fitStats->Nave ++; + + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", + average[0].R, + average[0].D, + fitAve.Ro, fitAve.Do, + 3600*(average[0].R - fitAve.Ro), + 3600*(average[0].D - fitAve.Do)); + + // make sure that the fit succeeded + int status = TRUE; + status &= finite(fitAve.Ro); + status &= finite(fitAve.Do); + status &= finite(fitAve.dRo); + status &= finite(fitAve.dDo); + if (!status) { + fitStats->Nskip ++; + return FALSE; + } + + // what is the offset relative to the mean fit position? + coords.crval1 = average[0].R; + coords.crval2 = average[0].D; + + double dXoff, dYoff; + RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords); + float dPos = hypot (dXoff, dYoff); + if (dPos > MaxMeanOffset) { + if (fitStats->Noffset < 100) { + fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff); + } + fitStats->Noffset ++; + return FALSE; + } + + // set the stack position values + average[0].Rstk = fitAve.Ro; // RA in degrees + average[0].Dstk = fitAve.Do; // DEC in degrees + average[0].dRstk = fitAve.dRo; // RA scatter in arcsec + average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec + + return (TRUE); +} + + /* fitting proper-motion and parallax: Index: /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/extra.c =================================================================== --- /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/extra.c (revision 37306) +++ /branches/eam_branches/ipp-20140813/Ohana/src/relastro/src/extra.c (revision 37306) @@ -0,0 +1,40 @@ +# include "relastro.h" + +// for now (20140710) I need to identify gpc1 chips explicitly. generalize in the future +int isGPC1chip (int photcode) { + + if ((photcode > 10000) && (photcode < 10077)) return TRUE; // g-band + if ((photcode > 10100) && (photcode < 10177)) return TRUE; // r-band + if ((photcode > 10200) && (photcode < 10277)) return TRUE; // i-band + if ((photcode > 10300) && (photcode < 10377)) return TRUE; // z-band + if ((photcode > 10400) && (photcode < 10477)) return TRUE; // y-band + if ((photcode > 10500) && (photcode < 10577)) return TRUE; // w-band + + return FALSE; +} + +// for now (20140710) I need to identify gpc1 stacks explicitly. generalize in the future +int isGPC1stack (int photcode) { + + if (photcode == 11000) return TRUE; // g-band + if (photcode == 11100) return TRUE; // r-band + if (photcode == 11200) return TRUE; // i-band + if (photcode == 11300) return TRUE; // z-band + if (photcode == 11400) return TRUE; // y-band + if (photcode == 11500) return TRUE; // w-band + + return FALSE; +} + +// for now (20140710) I need to identify gpc1 stacks explicitly. generalize in the future +int isGPC1warp (int photcode) { + + if (photcode == 12000) return TRUE; // g-band + if (photcode == 12100) return TRUE; // r-band + if (photcode == 12200) return TRUE; // i-band + if (photcode == 12300) return TRUE; // z-band + if (photcode == 12400) return TRUE; // y-band + if (photcode == 12500) return TRUE; // w-band + + return FALSE; +}