Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39171) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39225) @@ -254,4 +254,8 @@ int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored double MIN_ERROR; + +double MIN_DISTANCE_MOD; +double MAX_DISTANCE_MOD; +double MAX_DISTANCE_MOD_ERR; int IMFIT_TOO_FEW; // need more than this number of stars to fit an image @@ -323,4 +327,7 @@ PhotCode *DCR_RED_PHOTCODE_POS, *DCR_RED_PHOTCODE_NEG; int DCR_RED_NSEC_POS, DCR_RED_NSEC_NEG; + +float *LoopWeight2MASS; +float *LoopWeightTycho; int ImagSelect; Index: trunk/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 39171) +++ trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 39225) @@ -159,6 +159,8 @@ GET_COLUMN(uR, "U_RA", float); GET_COLUMN(uD, "U_DEC", float); - GET_COLUMN(duR, "V_RA_ERR", float); - GET_COLUMN(duD, "V_DEC_ERR", float); + GET_COLUMN(duR, "U_RA_ERR", float); + GET_COLUMN(duD, "U_DEC_ERR", float); + GET_COLUMN(uRgal, "U_RA_GAL", float); + GET_COLUMN(uDgal, "U_DEC_GAL", float); GET_COLUMN(P, "PAR", float); GET_COLUMN(dP, "PAR_ERR", float); @@ -189,4 +191,5 @@ ALLOCATE (average, Average, Nrow); for (i = 0; i < Nrow; i++) { + dvo_average_init (&average[i]); average[i].R = R[i] ; average[i].D = D[i] ; @@ -197,4 +200,6 @@ average[i].duR = duR[i] ; average[i].duD = duD[i] ; + average[i].uRgal = uRgal[i] ; + average[i].uDgal = uDgal[i] ; average[i].P = P[i] ; average[i].dP = dP[i] ; @@ -230,4 +235,6 @@ free (duR); free (duD); + free (uRgal); + free (uDgal); free (P); free (dP); @@ -489,6 +496,8 @@ gfits_define_bintable_column (&theader, "E", "U_RA", "RA*cos(D) proper-motion", "", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "U_DEC", "DEC proper-motion", "", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "V_RA_ERR", "RA*cos(D) p-m error", "", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "V_DEC_ERR", "DEC p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_RA_ERR", "RA*cos(D) p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_DEC_ERR", "DEC p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_RA_GAL", "RA*cos(D) p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_DEC_GAL", "DEC p-m error", "", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "PAR", "parallax", "", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "PAR_ERR", "parallax error", "", 1.0, 0.0); @@ -526,4 +535,6 @@ float *duR ; ALLOCATE (duR , float , catalog->Naverage); float *duD ; ALLOCATE (duD , float , catalog->Naverage); + float *uRgal ; ALLOCATE (uRgal , float , catalog->Naverage); + float *uDgal ; ALLOCATE (uDgal , float , catalog->Naverage); float *P ; ALLOCATE (P , float , catalog->Naverage); float *dP ; ALLOCATE (dP , float , catalog->Naverage); @@ -560,4 +571,6 @@ duR[i] = average[i].duR ; duD[i] = average[i].duD ; + uRgal[i] = average[i].uRgal ; + uDgal[i] = average[i].uDgal ; P[i] = average[i].P ; dP[i] = average[i].dP ; @@ -591,6 +604,8 @@ gfits_set_bintable_column (&theader, &ftable, "U_RA", uR, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "U_DEC", uD, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "V_RA_ERR", duR, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "V_DEC_ERR", duD, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_RA_ERR", duR, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_DEC_ERR", duD, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_RA_GAL", uRgal, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_DEC_GAL", uDgal, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "PAR", P, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "PAR_ERR", dP, catalog->Naverage); @@ -604,5 +619,5 @@ gfits_set_bintable_column (&theader, &ftable, "NMEASURE", Nmeasure, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "NMISSING", Nmissing, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "NGALPHOT", Ngalphot, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NGALPHOT", Ngalphot, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset, catalog->Naverage); @@ -624,4 +639,6 @@ free (duR); free (duD); + free (uRgal); + free (uDgal); free (P); free (dP); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39171) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39225) @@ -35,4 +35,8 @@ if (!ScanConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX)) DPOS_MAX = 6.0; if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; + + if (!ScanConfig (config, "RELASTRO_MIN_DISTANCE_MOD", "%lf", 0, &MIN_DISTANCE_MOD)) MIN_DISTANCE_MOD = 7.5; + if (!ScanConfig (config, "RELASTRO_MAX_DISTANCE_MOD", "%lf", 0, &MAX_DISTANCE_MOD)) MAX_DISTANCE_MOD = 15.0; + if (!ScanConfig (config, "RELASTRO_MAX_DISTANCE_MOD_ERR", "%lf", 0, &MAX_DISTANCE_MOD_ERR)) MAX_DISTANCE_MOD_ERR = 1.0; if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) USE_FIXED_PIXCOORDS = FALSE; Index: trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- trunk/Ohana/src/relastro/src/GetAstromError.c (revision 39171) +++ trunk/Ohana/src/relastro/src/GetAstromError.c (revision 39225) @@ -73,74 +73,16 @@ int has2MASS = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS); - // serious hack: if the object has 2MASS, we set this internal bit and adjust the + // modest hack: if the object has 2MASS or Tycho, we set this internal bit and adjust the // weight to ensure the image is tied down to the 2mass frame - if (has2MASS) { - switch (Nloop) { - case 0: - dPtotal = dPtotal / 1000.0; - break; - case 1: - dPtotal = dPtotal / 300.0; - break; - case 2: - dPtotal = dPtotal / 100.0; - break; - case 3: - dPtotal = dPtotal / 30.0; - break; - case 4: - dPtotal = dPtotal / 10.0; - break; - case 5: - dPtotal = dPtotal / 10.0; - break; - default: - break; - } + if (has2MASS && LoopWeight2MASS) { + dPtotal = dPtotal / LoopWeight2MASS[Nloop]; } - if (is2MASS) { - switch (Nloop) { - case 0: - dPtotal = dPtotal / 1000.0; - break; - case 1: - dPtotal = dPtotal / 300.0; - break; - case 2: - dPtotal = dPtotal / 100.0; - break; - case 3: - dPtotal = dPtotal / 30.0; - break; - case 4: - dPtotal = dPtotal / 10.0; - break; - case 5: - dPtotal = dPtotal / 10.0; - break; - default: - break; - } + if (is2MASS && LoopWeight2MASS) { + dPtotal = dPtotal / LoopWeight2MASS[Nloop]; } - if (isTycho) { - switch (Nloop) { - case 0: - case 1: - dPtotal = dPtotal / 200.0; - break; - case 2: - case 3: - dPtotal = dPtotal / 100.0; - break; - case 4: - case 5: - dPtotal = dPtotal / 50.0; - break; - default: - break; - } + if (isTycho && LoopWeightTycho) { + dPtotal = dPtotal / LoopWeightTycho[Nloop]; } - return (dPtotal); } @@ -197,8 +139,29 @@ dPtotal = MAX (dPtotal, MIN_ERROR); + + // early on, we want 2MASS and Tycho to have a very high weight. This will force images + // to match the 2MASS / Tycho / ICRS reference frame. As Nloop gets higher, the weight + // needs to drop to allow the ps1 measurements to drive the solution + int is2MASS = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2011) && (measure[0].photcode <= 2013); + int isTycho = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2020) && (measure[0].photcode <= 2021); + int has2MASS = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS); + + // modest hack: if the object has 2MASS or Tycho, we set this internal bit and adjust the + // weight to ensure the image is tied down to the 2mass frame + + if (has2MASS && LoopWeight2MASS) { + dPtotal = dPtotal / LoopWeight2MASS[Nloop]; + } + if (is2MASS && LoopWeight2MASS) { + dPtotal = dPtotal / LoopWeight2MASS[Nloop]; + } + if (isTycho && LoopWeightTycho) { + dPtotal = dPtotal / LoopWeightTycho[Nloop]; + } + return (dPtotal); } -/* for a long time, psphot was either not reported position errors, or was reporting +/* for a long time, psphot was either not reporting position errors, or was reporting * completely wrong astrometry errors. This function lets us handle, in the * configuration, different strategies to generating a position error Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 39171) +++ trunk/Ohana/src/relastro/src/args.c (revision 39225) @@ -4,4 +4,5 @@ void usage_merge_source (void); void usage_merge_source_id (char *name); +float *ParseLoopWeights (char *rawlist); int args (int argc, char **argv) { @@ -529,4 +530,19 @@ remove_argument (N, &argc, argv); NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + // e.g., -loop-weights-2mass 1000,300,300,200,200,100 + // NOTE: this must come after -nloop above + LoopWeight2MASS = NULL; + if ((N = get_argument (argc, argv, "-loop-weights-2mass"))) { + remove_argument (N, &argc, argv); + LoopWeight2MASS = ParseLoopWeights (argv[N]); + remove_argument (N, &argc, argv); + } + LoopWeightTycho = NULL; + if ((N = get_argument (argc, argv, "-loop-weights-tycho"))) { + remove_argument (N, &argc, argv); + LoopWeightTycho = ParseLoopWeights (argv[N]); remove_argument (N, &argc, argv); } @@ -999,2 +1015,33 @@ } +float *ParseLoopWeights (char *rawlist) { + + float *weights = NULL; + ALLOCATE (weights, float, NLOOP); + + int Nloop = 0; + + /* parse the comma-separated list of photcodes */ + char *myList = strcreate(rawlist); + char *list = myList; + char *entry = NULL; + char *ptr = NULL; + while ((Nloop < NLOOP) && ((entry = strtok_r (list, ",", &ptr)) != NULL)) { + list = NULL; // pass NULL on successive strtok_r calls + + weights[Nloop] = atof(entry); + Nloop ++; + } + free (myList); + + if (Nloop == 0) { + fprintf (stderr, "syntax error parsing weights: %s\n", rawlist); + exit (3); + } + + while (Nloop < NLOOP) { + weights[Nloop] = weights[Nloop - 1]; + Nloop ++; + } + return weights; +} Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 39171) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 39225) @@ -117,4 +117,13 @@ if (!isfinite(catalog[0].average[i].uRgal)) continue; if (!isfinite(catalog[0].average[i].uDgal)) continue; + + // filter stars which are not well-handled by the galactic motion model + int mStarpar = catalog[0].average[i].starparOffset; + myAssert (mStarpar < catalog[0].Nstarpar, "oops"); + + StarPar *starpar = &catalog[0].starpar[mStarpar]; + if (starpar->DistMag < MIN_DISTANCE_MOD) continue; + if (starpar->DistMag > MAX_DISTANCE_MOD) continue; + if (starpar->dDistMag > MAX_DISTANCE_MOD_ERR) continue; } Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 39171) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 39225) @@ -134,4 +134,21 @@ if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); + if (LoopWeight2MASS) { + char *weightline = NULL; + for (i = 0; i < NLOOP; i++) { + strextend (&weightline, "%f,", LoopWeight2MASS[i]); + } + strextend (&command, "-loop-weights-2mass %s", weightline); + free (weightline); + } + if (LoopWeightTycho) { + char *weightline = NULL; + for (i = 0; i < NLOOP; i++) { + strextend (&weightline, "%f,", LoopWeightTycho[i]); + } + strextend (&command, "-loop-weights-tycho %s", weightline); + free (weightline); + } + strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 39171) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 39225) @@ -44,4 +44,8 @@ pcatalog->catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT; // don't need to load all data at this point pcatalog->Nsecfilt = GetPhotcodeNsecfilt (); + + if (USE_GALAXY_MODEL) { + pcatalog->catflags = pcatalog->catflags | DVO_LOAD_STARPAR; + } // loads Average, Measure, SecFilt @@ -171,4 +175,7 @@ strextend (&command, " -minerror %f", MIN_ERROR); strextend (&command, " -D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, " -D RELASTRO_MIN_DISTANCE_MOD %f", MIN_DISTANCE_MOD); + strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD %f", MAX_DISTANCE_MOD); + strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD_ERR %f", MAX_DISTANCE_MOD_ERR); if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm");