Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37564) @@ -87,15 +87,17 @@ int VERBOSE; int VERBOSE2; +int TESTING; int FORCE; +int UNIFORM_RADEC; -int FAKEASTRO_NLOOP; -int FAKEASTRO_NSTARS; -int FAKEASTRO_NQSO_ICRF; -int FAKEASTRO_NQSO_ZERO; -float FAKEASTRO_ZGAL; // parsec -float FAKEASTRO_RGAL; // parsec -char FAKEASTRO_REF_EPOCH[80]; -char FAKEASTRO_2MASS_EPOCH[80]; +int FAKEASTRO_NLOOP; +int FAKEASTRO_NSTARS; +int FAKEASTRO_NQSO_ICRF; +int FAKEASTRO_NQSO_ZERO; +float FAKEASTRO_ZGAL; // parsec +float FAKEASTRO_RGAL; // parsec +char FAKEASTRO_REF_EPOCH[80]; +char FAKEASTRO_2MASS_EPOCH[80]; float RADIUS; Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37564) @@ -19,5 +19,5 @@ // if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; - if (FAKEASTRO_OP == OP_GALAXY) { + if ((FAKEASTRO_OP == OP_GALAXY) || (FAKEASTRO_OP == OP_2MASS)) { // force CATDIR to be absolute (so parallel mode will work) GetConfig (config, "CATDIR", "%s", 0, CATDIR); @@ -73,12 +73,15 @@ if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); - char *CATDIR_CHECK = (FAKEASTRO_OP == OP_GALAXY) ? CATDIR : CATDIR_OUTPUT; + char *CATDIR_CHECK = (FAKEASTRO_OP == OP_IMAGES) ? CATDIR_OUTPUT : CATDIR; - // check for existence of CATDIR - struct stat filestat; - int status = stat (CATDIR_CHECK, &filestat); - if (!FORCE && (status == 0)) { - fprintf (stderr, "directory %s exists, refusing to contaminate\n", CATDIR); - exit (1); + // OP_2MASS is adding detections to an existing db, the others require and empty db + if (FAKEASTRO_OP != OP_2MASS) { + // check for existence of CATDIR + struct stat filestat; + int status = stat (CATDIR_CHECK, &filestat); + if (!FORCE && (status == 0)) { + fprintf (stderr, "directory %s exists, refusing to contaminate\n", CATDIR); + exit (1); + } } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37564) @@ -19,4 +19,15 @@ remove_argument (N, argc, argv); FAKEASTRO_OP = OP_GALAXY; + } + + if ((N = get_argument (*argc, argv, "-2mass"))) { + remove_argument (N, argc, argv); + FAKEASTRO_OP = OP_2MASS; + } + + TESTING = FALSE; + if ((N = get_argument (*argc, argv, "-testing"))) { + remove_argument (N, argc, argv); + TESTING = TRUE; } @@ -84,4 +95,10 @@ // MaxDensityUse = TRUE; // } + + UNIFORM_RADEC = FALSE; + if ((N = get_argument (*argc, argv, "-uniform-radec"))) { + UNIFORM_RADEC = TRUE; + remove_argument (N, argc, argv); + } FORCE = FALSE; Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_2mass.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_2mass.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_2mass.c (revision 37564) @@ -23,5 +23,5 @@ catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR; catalog.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) { + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); exit (1); @@ -37,8 +37,14 @@ SetProtect (TRUE); - dvo_catalog_update (&catalog, VERBOSE); + if (!dvo_catalog_update (&catalog, VERBOSE)) { + fprintf (stderr, "ERROR: failure to update %s\n", catalog.filename); + exit (3); + } SetProtect (FALSE); - dvo_catalog_unlock (&catalog); + if (!dvo_catalog_unlock (&catalog)) { + fprintf (stderr, "ERROR: failure to unlock %s\n", catalog.filename); + exit (2); + } dvo_catalog_free (&catalog); if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_2mass_measures.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_2mass_measures.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_2mass_measures.c (revision 37564) @@ -44,4 +44,16 @@ double J_PS1 = secfilt[i*Nsecfilt + Nsec].M - 1.0; // make all stars have y-J = 1.0 if (J_PS1 > 16.0) continue; // only generate 2MASS detections for objects with J_PS1 < 16.0 + + if (isnan(secfilt[i*Nsecfilt + Nsec].M)) { + // look for a non-NAN secfilt mag and just use that (it is not super important) + int ns; + for (ns = 0; ns < Nsecfilt; ns++) { + if (!isnan(secfilt[i*Nsecfilt + ns].M)) break; + } + if (ns == Nsecfilt) continue; // no non-nan + + J_PS1 = secfilt[i*Nsecfilt + ns].M; // pretend secfilt.M = J + if (J_PS1 > 16.0) continue; // only generate 2MASS detections for objects with J_PS1 < 16.0 + } double Minst = J_PS1 - ZP; @@ -101,5 +113,5 @@ // This is may optionally be replaced by the internal sequence (see FilterStars.c) - measure[Nmeasure].detID = Nmeasure; + measure[Nmeasure].detID = 0; Nmeasure ++; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakeqsos.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakeqsos.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakeqsos.c (revision 37564) @@ -12,4 +12,6 @@ float dPmax = sin(RAD_DEG*UserPatch.Dmax); float dP = dPmax - dPmin; + + CoordTransform *transback = InitTransform (COORD_CELESTIAL, COORD_GALACTIC); int i; @@ -32,7 +34,9 @@ break; } + double L, B; + ApplyTransform (&L, &B, R, D, transback); // Mr will be interpretted as m_r - double Mr = rnd_gauss (17.0, 1.0); + double Mr = rnd_gauss (18.0, 1.5); stars[i].R = R; @@ -43,6 +47,6 @@ stars[i].starpar.R = R; stars[i].starpar.D = D; - stars[i].starpar.galLon = NAN; - stars[i].starpar.galLat = NAN; + stars[i].starpar.galLon = L; + stars[i].starpar.galLat = B; // how shall I distinguish ICRF and ZERO quasars Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37564) @@ -19,4 +19,9 @@ CoordTransform *transback = InitTransform (COORD_CELESTIAL, COORD_GALACTIC); + float dR = UserPatch.Rmax - UserPatch.Rmin; + float dPmin = sin(RAD_DEG*UserPatch.Dmin); + float dPmax = sin(RAD_DEG*UserPatch.Dmax); + float dP = dPmax - dPmin; + int i; for (i = 0; i < Nstars; i++) { @@ -35,20 +40,35 @@ double z,r,L,B,R,D,Lrad,Brad; - int inPatch = FALSE; - while (!inPatch) { - z = rnd_gauss (0.0, FAKEASTRO_ZGAL); - r = sqrt(drand48()) * FAKEASTRO_RGAL; - Lrad = drand48() * 2 * M_PI; - Brad = atan2(z,r); + if (UNIFORM_RADEC) { + // we can generate a distribution which is uniform on the sky, in which case + // we can limit to the selected patch analyically + r = pow(drand48(), 0.33333) * FAKEASTRO_RGAL; + z = 0.0; + double phi = dP * drand48() + dPmin; + D = DEG_RAD * asin(phi); // random in degrees + R = drand48() * dR + UserPatch.Rmin; // random in degrees + R = ohana_normalize_angle (R); + + ApplyTransform (&L, &B, R, D, transback); + Lrad = L*RAD_DEG; + Brad = B*RAD_DEG; + } else { + int inPatch = FALSE; + while (!inPatch) { + z = rnd_gauss (0.0, FAKEASTRO_ZGAL); + r = sqrt(drand48()) * FAKEASTRO_RGAL; + Lrad = drand48() * 2 * M_PI; + Brad = atan2(z,r); - L = Lrad*DEG_RAD; - B = Brad*DEG_RAD; + L = Lrad*DEG_RAD; + B = Brad*DEG_RAD; - ApplyTransform (&R, &D, L, B, transform); - if (R < UserPatch.Rmin) continue; - if (R > UserPatch.Rmax) continue; - if (D < UserPatch.Dmin) continue; - if (D > UserPatch.Dmax) continue; - break; + ApplyTransform (&R, &D, L, B, transform); + if (R < UserPatch.Rmin) continue; + if (R > UserPatch.Rmax) continue; + if (D < UserPatch.Dmin) continue; + if (D > UserPatch.Dmax) continue; + break; + } } @@ -66,9 +86,25 @@ double uL = uL_gal + uL_sol; double uB = uB_gal + uB_sol; - + + // XXX: amplify motion to make tests easier: + if (TESTING) { + uL *= 100.0; + uB *= 100.0; + } + // crude Mr distribution from Bochanski et al 2010 // http://iopscience.iop.org/1538-3881/139/6/2679/pdf/aj_139_6_2679.pdf - - double Mr = rnd_gauss (11.25, 1.0); + + // two gaussian distributions: + // 75% in narrow gauss, with Mr = 11.25, sigma = 1.0 + // 25% is wide gauss, with Mr = 10.0, sigma = 3.0 + // first choose which gauss: + int bigPeak = (drand48() > 0.60); + double Mr; + if (bigPeak) { + Mr = rnd_gauss (11.25, 1.0); + } else { + Mr = rnd_gauss (10.00, 2.5); + } // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c (revision 37563) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c (revision 37564) @@ -252,4 +252,7 @@ dvo_starpar_init (&catalog[0].starpar[Nstarpar]); catalog[0].starpar[Nstarpar] = stars[i].starpar; + catalog[0].starpar[Nstarpar].averef = Nave; + catalog[0].starpar[Nstarpar].objID = catalog[0].average[Nave].objID; + catalog[0].starpar[Nstarpar].catID = catalog[0].catID; Nstarpar ++;