Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/include/relastro.h =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/include/relastro.h (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/include/relastro.h (revision 36808) @@ -83,7 +83,8 @@ double L, M; /* Focal Plane - pixels */ double X, Y; /* Chip Coords - pixels */ - double Mag; - double dMag; - double dPos; + float Mag; + float Color; + float dMag; + float dPos; int mask; int Nmeas; @@ -216,4 +217,8 @@ char WHERE_A[10000], WHERE_B[10000]; SkyRegionSelection SELECTION; + +char *DCR_COLOR_POS, *DCR_COLOR_NEG; +PhotCode *DCR_PHOTCODE_POS, *DCR_PHOTCODE_NEG; +int DCR_NSEC_POS, DCR_NSEC_NEG; int ImagSelect; @@ -526,2 +531,4 @@ int ImageTableSave (char *filename, Image *images, off_t Nimages); int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); + +float getColor (off_t meas, int cat); Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/BrightCatalog.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/BrightCatalog.c (revision 36808) @@ -167,5 +167,5 @@ GET_COLUMN(measureOffset, "OFF_MEASURE", int); GET_COLUMN(missingOffset, "OFF_MISSING", int); - GET_COLUMN(extendOffset, "OFF_EXTEND", int); + GET_COLUMN(refColor, "REF_COLOR", float); GET_COLUMN(flags, "FLAGS", int); GET_COLUMN(photFlagsUpper, "PHOTFLAGS_U", int); @@ -202,5 +202,5 @@ average[i].measureOffset = measureOffset[i] ; average[i].missingOffset = missingOffset[i] ; - average[i].extendOffset = extendOffset[i] ; + average[i].refColor = refColor[i] ; average[i].flags = flags[i] ; average[i].photFlagsUpper = photFlagsUpper[i] ; @@ -234,5 +234,5 @@ free (measureOffset); free (missingOffset); - free (extendOffset); + free (refColor); free (flags); free (photFlagsUpper); @@ -510,5 +510,5 @@ int *measureOffset ; ALLOCATE (measureOffset , int , catalog->Naverage); int *missingOffset ; ALLOCATE (missingOffset , int , catalog->Naverage); - int *extendOffset ; ALLOCATE (extendOffset , int , catalog->Naverage); + float *refColor ; ALLOCATE (refColor , float , catalog->Naverage); int *flags ; ALLOCATE (flags , int , catalog->Naverage); int *photFlagsUpper; ALLOCATE (photFlagsUpper, int , catalog->Naverage); @@ -543,5 +543,5 @@ measureOffset[i] = average[i].measureOffset ; missingOffset[i] = average[i].missingOffset ; - extendOffset[i] = average[i].extendOffset ; + refColor[i] = average[i].refColor ; flags[i] = average[i].flags ; photFlagsUpper[i] = average[i].photFlagsUpper ; @@ -575,5 +575,5 @@ gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "OFF_EXTEND", extendOffset, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "REF_COLOR", refColor, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_U", photFlagsUpper, catalog->Naverage); @@ -605,5 +605,5 @@ free (measureOffset); free (missingOffset); - free (extendOffset); + free (refColor); free (flags); free (photFlagsUpper); Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/FitChip.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/FitChip.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/FitChip.c (revision 36808) @@ -4,5 +4,5 @@ int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) { - int i, NstatFull, Nstat, Niter; + int i, NstatFull, Nstat, Niter, Nkeep; float dLsig, dMsig, dRsig; float dLsigFull, dMsigFull, dRsigFull; @@ -24,4 +24,5 @@ // measure the scatter distribution (use only the bright end detections) + Nkeep = 0; for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; @@ -42,5 +43,9 @@ continue; } - } + Nkeep ++; + } + + // I'm rejecting some points from the fit above; I count the remainders and actually + // use that count to set order_use below // figures to assess the fitting process: @@ -50,29 +55,15 @@ relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch); - // set the maximum order for the polynomial (based on number of stars) + // set the maximum order for the polynomial (based on number of stars kept above) int order_use = 0; - if (Nmatch > 5) order_use = 1; // 4 stars per polynomial term (per dimension) - if (Nmatch > 30) order_use = 2; // 5 stars per polynomial term (per dimension) - if (Nmatch > 60) order_use = 3; // 6 stars per polynomial term (per dimension) + if (Nkeep > 5) order_use = 1; // 4 stars per polynomial term (per dimension) + if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension) + if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension) if (order_use < 1) { - if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch); + if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep); image[0].flags |= ID_IMAGE_ASTROM_FEW; return FALSE; } image[0].coords.Npolyterms = order_use; - -# if (0) - if (CHIPORDER > 0) { - int Ndof_min = 3; - int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3); - int order_use = MIN (CHIPORDER, order_max); - if (order_use < 1) { - if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch); - image[0].flags |= ID_IMAGE_ASTROM_FEW; - return FALSE; - } - image[0].coords.Npolyterms = order_use; - } -# endif if (fit) fit_free (fit); @@ -84,30 +75,4 @@ fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); } - -# if (0) - // check if the fit has enough data points for the polynomial order - skip = FALSE; - switch (image[0].coords.Npolyterms) { - case 0: - case 1: - skip = (fit[0].Npts < 8); - break; - case 2: - skip = (fit[0].Npts < 11); - break; - case 3: - skip = (fit[0].Npts < 15); - break; - default: - fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms); - skip = TRUE; - } - if (skip) { - if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms); - fit_free (fit); - image[0].flags |= ID_IMAGE_ASTROM_FEW; - return FALSE; - } -# endif // measure the fit, update the coords & object coordinates @@ -147,4 +112,20 @@ } + int Ncolor = 0; + float *colorList = NULL; + ALLOCATE (colorList, float, Nmatch); + + for (i = 0; i < Nmatch; i++) { + if (!raw[i].mask) continue; + if (!isnan(ref[i].Color)) { + colorList[Ncolor] = ref[i].Color; + Ncolor ++; + } + } + fsort (colorList, Ncolor); + float colorMedian = (Ncolor > 0) ? colorList[(int)(0.5*Ncolor)] : NAN; + image[0].refColor = colorMedian; + free (colorList); + GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM); GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM); @@ -164,4 +145,6 @@ image[0].dYpixSys = dMsig; image[0].nFitAstrom = fit[0].Npts; + + // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom); if (fit) fit_free (fit); Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/ImageOps.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/ImageOps.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/ImageOps.c (revision 36808) @@ -333,4 +333,13 @@ */ +float getColor (off_t meas, int cat) { + + off_t i; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (NAN); + return (image[i].refColor); +} + void plot_images () { @@ -738,4 +747,6 @@ StarData *ref; + int Nsecfilt = GetPhotcodeNsecfilt(); + mosaic = NULL; moscoords = NULL; @@ -765,4 +776,8 @@ ref[i].dMag = measure[0].dM; ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); + + if ((DCR_NSEC_POS >= 0) && (DCR_NSEC_NEG >= -1)) { + ref[i].Color = catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_POS].M - catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_NEG].M; + } ref[i].mask = FALSE; Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/UpdateObjects.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/UpdateObjects.c (revision 36808) @@ -10,4 +10,5 @@ static double *T; static double *dT; +static double *C; void initObjectData (Catalog *catalog, int Ncatalog) { @@ -36,4 +37,6 @@ ALLOCATE (pX, double, MAX (1, Nmax)); ALLOCATE (pY, double, MAX (1, Nmax)); + + ALLOCATE (C, double, MAX (1, Nmax)); } @@ -54,4 +57,6 @@ free (pX); free (pY); + + free (C); } @@ -61,5 +66,5 @@ off_t j, k, m; - int i, N, Nsecfilt, mode, result, status, XVERB; + int i, N, Nc, Nsecfilt, mode, result, status, XVERB; StatType statsR, statsD; Coords coords; @@ -116,4 +121,5 @@ } + Nc = 0; N = 0; m = catalog[i].average[j].measureOffset; @@ -222,4 +228,10 @@ dD[N] = dY[N] / 3600.0; + float color = getColor (m+k, i); + if (!isnan(color)) { + C[Nc] = color; + Nc++; + } + N++; } @@ -341,4 +353,8 @@ Nave ++; } + + dsort (C, Nc); + float colorMedian = (Nc > 0) ? C[(int)(0.5*Nc)] : NAN; + catalog[i].average[j].refColor = colorMedian; /* choose the result based on the chisq values */ Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/args.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/args.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/args.c (revision 36808) @@ -292,4 +292,14 @@ remove_argument (N, &argc, argv); TimeSelect = TRUE; + } + + DCR_COLOR_POS = NULL; + DCR_COLOR_NEG = NULL; + if ((N = get_argument (argc, argv, "-dcr-color"))) { + remove_argument (N, &argc, argv); + DCR_COLOR_POS = strcreate(argv[N]); + remove_argument (N, &argc, argv); + DCR_COLOR_NEG = strcreate(argv[N]); + remove_argument (N, &argc, argv); } Index: branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/initialize.c =================================================================== --- branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/initialize.c (revision 36807) +++ branches/eam_branches/ipp-20140423/Ohana/src/relastro/src/initialize.c (revision 36808) @@ -7,4 +7,7 @@ if (RELASTRO_OP == OP_MERGE_SOURCE) return; + + if (DCR_COLOR_POS) fprintf (stderr, "DCR_COLOR_POS: %s\n", DCR_COLOR_POS); + if (DCR_COLOR_NEG) fprintf (stderr, "DCR_COLOR_NEG: %s\n", DCR_COLOR_NEG); if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); @@ -19,4 +22,22 @@ photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); + + DCR_NSEC_POS = DCR_NSEC_NEG = -1; + if (DCR_COLOR_POS) { + DCR_PHOTCODE_POS = GetPhotcodebyName (DCR_COLOR_POS); + if (!DCR_PHOTCODE_POS) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_POS); + exit (1); + } + DCR_NSEC_POS = GetPhotcodeNsec (DCR_PHOTCODE_POS[0].code); + } + if (DCR_COLOR_NEG) { + DCR_PHOTCODE_NEG = GetPhotcodebyName (DCR_COLOR_NEG); + if (!DCR_PHOTCODE_NEG) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_NEG); + exit (1); + } + DCR_NSEC_NEG = GetPhotcodeNsec (DCR_PHOTCODE_NEG[0].code); + } initstats (STATMODE);