Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37523) @@ -53,4 +53,6 @@ $(SRC)/save_fake_stars.$(ARCH).o \ $(SRC)/match_fake_stars.$(ARCH).o \ +$(SRC)/fakeastro_2mass.$(ARCH).o \ +$(SRC)/make_2mass_measures.$(ARCH).o \ $(SRC)/remote_hosts.$(ARCH).o Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37523) @@ -6,5 +6,7 @@ # include -typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES} FakeastroOp; +# define RESETTIME { gettimeofday (&startTimer, (void *) NULL); } + +typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES, OP_2MASS} FakeastroOp; typedef struct { @@ -30,4 +32,5 @@ Average average; Measure measure; + StarPar starpar; int found; } Stars; @@ -94,4 +97,5 @@ float FAKEASTRO_RGAL; // parsec char FAKEASTRO_REF_EPOCH[80]; +char FAKEASTRO_2MASS_EPOCH[80]; float RADIUS; @@ -188,2 +192,4 @@ FakeAstro_Stars *make_fakeqsos (int Nstars, int ICRF); +int fakeastro_2mass (); +int make_2mass_measures (Catalog *catalog); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37523) @@ -57,4 +57,7 @@ strcpy (FAKEASTRO_REF_EPOCH, "2011/05/11,00:00:00"); // epoch of reference catalog } + if (!ScanConfig (config, "FAKEASTRO_2MASS_EPOCH", "%s", 0, FAKEASTRO_2MASS_EPOCH)) { + strcpy (FAKEASTRO_2MASS_EPOCH, "2000/01/01,00:00:00"); // epoch of reference catalog + } /* set the default search radius */ Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37523) @@ -11,13 +11,32 @@ case OP_GALAXY: /* the object analysis is a separate process iterating over catalogs */ - fakeastro_galaxy (); + fakeastro_galaxy (); // generate fake stars following a galaxy model + /* make_fakestars() -- generates FakeAstro_Stars with density,distance and uR,uD following galaxy model + make_fakeqsos() + sortStars() + make_subset() + save_fakestars() + fakestar_catalog() + insert_fakestar() -- converts suppied FakeAstro_Stars and populates average/secfilt/starpar structures + (fakestar_save_stars) [fakeastro_client / fakestar_load_stars] + */ exit (0); case OP_IMAGES: fakeastro_images (); + /* fakeastro_images_region() + load_fake_stars() + make_fake_images() + make_fake_stars() + fit_fake_stars() + save_fake_stars() + match_fake_stars() + */ exit (0); case OP_2MASS: fakeastro_2mass (); + /* make_2mass_measures() + */ exit (0); 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 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_2mass.c (revision 37523) @@ -8,5 +8,5 @@ SkyTableSetFilenames (skyTable, CATDIR, "cpt"); - SkyList *skyList = SkyListByPatch (skyTable, -1, &UserPatch); + SkyList *skylist = SkyListByPatch (skyTable, -1, &UserPatch); Catalog catalog; @@ -21,5 +21,5 @@ catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR; + catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR; catalog.Nsecfilt = GetPhotcodeNsecfilt (); if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) { @@ -43,6 +43,5 @@ dvo_catalog_free (&catalog); if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; - - + } exit (0); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c (revision 37523) @@ -27,6 +27,6 @@ gfits_define_bintable_column (&theader, "D", "RA", "", "degree", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "DEC", "", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "GLON", "", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "GLAT", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "GLON", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "GLAT", "", "degree", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "Ebv", "", "", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "dEbv", "", "", 1.0, 0.0); @@ -44,6 +44,6 @@ double *ra = NULL; ALLOCATE (ra , double , Nstars); for (i = 0; i < Nstars; i++) ra [i] = stars[i].R ; double *dec = NULL; ALLOCATE (dec , double , Nstars); for (i = 0; i < Nstars; i++) dec [i] = stars[i].D ; - double *glon = NULL; ALLOCATE (glon , double , Nstars); for (i = 0; i < Nstars; i++) glon [i] = stars[i].starpar.galLon ; - double *glat = NULL; ALLOCATE (glat , double , Nstars); for (i = 0; i < Nstars; i++) glat [i] = stars[i].starpar.galLat ; + float *glon = NULL; ALLOCATE (glon , float , Nstars); for (i = 0; i < Nstars; i++) glon [i] = stars[i].starpar.galLon ; + float *glat = NULL; ALLOCATE (glat , float , Nstars); for (i = 0; i < Nstars; i++) glat [i] = stars[i].starpar.galLat ; float *Ebv = NULL; ALLOCATE (Ebv , float , Nstars); for (i = 0; i < Nstars; i++) Ebv [i] = stars[i].starpar.Ebv ; float *dEbv = NULL; ALLOCATE (dEbv , float , Nstars); for (i = 0; i < Nstars; i++) dEbv [i] = stars[i].starpar.dEbv ; @@ -135,6 +135,6 @@ GET_COLUMN(ra , "RA" , double); GET_COLUMN(dec , "DEC" , double); - GET_COLUMN(glon , "GLON" , double); - GET_COLUMN(glat , "GLAT" , double); + GET_COLUMN(glon , "GLON" , float); + GET_COLUMN(glat , "GLAT" , float); GET_COLUMN(Ebv , "Ebv" , float); GET_COLUMN(dEbv , "dEbv" , float); @@ -159,4 +159,6 @@ stars[i].R = ra [i]; stars[i].D = dec [i]; + stars[i].starpar.R = ra [i]; + stars[i].starpar.D = dec [i]; stars[i].starpar.galLon = glon [i]; stars[i].starpar.galLat = glat [i]; Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c (revision 37523) @@ -12,4 +12,9 @@ if (Nstars > 225) order_use = 5; // 9 stars per cell if (Nstars > 360) order_use = 6; // 10 stars per cell + if (order_use < 1) { + fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nstars); + image[0].flags |= ID_IMAGE_ASTROM_FEW; + return FALSE; + } // I have Nstars with "reference" positions stars[].Rref,Dref and @@ -19,5 +24,5 @@ int CHIPMAP = (image[0].coords.Npolyterms == -1); - image[0].coords.Npolyterms = CHIPMAP ? -1 : order_use; + image[0].coords.Npolyterms = CHIPMAP ? 1 : order_use; CoordFit *fit = fit_init (image[0].coords.Npolyterms); 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 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_2mass_measures.c (revision 37523) @@ -3,4 +3,6 @@ // region corresponds to the catalog int make_2mass_measures (Catalog *catalog) { + + int Nsecfilt = GetPhotcodeNsecfilt (); Average *average = catalog->average; @@ -9,11 +11,40 @@ Measure *measure = catalog->measure; + time_t timeRef = ohana_date_to_sec(FAKEASTRO_REF_EPOCH); + time_t tzero_2mass = ohana_date_to_sec(FAKEASTRO_2MASS_EPOCH); + + // tzero, timeRef are in UNIX seconds, Toffset should be in years + double Toffset = (tzero_2mass - timeRef) / 365.25 / 86400.0; + + // use photcode to get zero point + PhotCode *codeJ = GetPhotcodebyName ("2MASS_J"); + // PhotCode *codeH = GetPhotcodebyName ("2MASS_H"); + // PhotCode *codeK = GetPhotcodebyName ("2MASS_K"); + + float ZP = 21.0; // this is chosen to give 16.0 mag stars a flux of 100.0 counts (crude, yes) + float SkyCts = 300.0; // this is chosen to make SN @ 16.0 = 5.0 ( + + int Nmeasure = catalog->Nmeasure; + int NMEASURE = catalog->Nmeasure; + + int Nsec = 4; // y-band + int i; for (i = 0; i < catalog->Naverage; i++) { - int nStar = average[i].starparOffset; + if (Nmeasure >= NMEASURE) { + NMEASURE = Nmeasure + 1000; + REALLOCATE (catalog[0].measure, Measure, NMEASURE); + measure = catalog[0].measure; + } - // which filter? - double y_PS1 = secfilt[i*Nsecfilt + 4].M - ZP; + int Nstarpar = average[i].starparOffset; + + // make a crude JHK color for now: + 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 + + double Minst = J_PS1 - ZP; + double Counts = pow(10.0, -0.4*Minst); @@ -21,17 +52,13 @@ // true position from src catalog - // XXX need to save the true position somewhere... - double Rtru = average[i].R; - double Dtru = average[i].D; + double Rtru = starpar[Nstarpar].R; + double Dtru = starpar[Nstarpar].D; // observed position is scattered from true position by: // * proper motion // * gaussian scatter (~ seeing) - double uR = starpar[nStar].uRA; // starpar are in arcsec / year - double uD = starpar[nStar].uDEC; + double uR = starpar[Nstarpar].uRA; // starpar are in arcsec / year + double uD = starpar[Nstarpar].uDEC; - // tzero, timeRef are in UNIX seconds, Toffset should be in years - double Toffset = (tzero_2mass - timeRef) / 365.25 / 86400.0; - // uR,uD in linear (arcsec / yr) double dRpm = uR*Toffset; @@ -48,25 +75,37 @@ double Dobs = Dtru + dDoff; - measure[nMeasure].M = Minst + ZP; - measure[nMeasure].dM = 1.0 / SN; + dvo_measure_init (&measure[Nmeasure]); - // stars[Nstars].measure.dMcal = ps1data[i].dMcal; - measure[nMeasure].Sky = sky[Nsec]; - measure[nMeasure].dSky = sqrt(sky[Nsec]); - measure[nMeasure].photFlags = 0; - measure[nMeasure].photFlags2 = 0; - measure[nMeasure].R = Robs; - measure[nMeasure].D = Dobs; - measure[nMeasure].airmass = airmass (image->secz, stars[Nstars].average.R, stars[Nstars].average.D, image->sidtime, image->latitude); - measure[nMeasure].az = azimuth (15.0*image->sidtime - stars[Nstars].average.R, stars[Nstars].average.D, image->latitude); - measure[nMeasure].Mcal = image->Mcal; - measure[nMeasure].t = image->tzero + 1e-4*stars[Nstars].measure.Yccd*image->trate; // trate is in 0.1 msec / row - measure[nMeasure].dt = Mtime; + measure[Nmeasure].R = Robs; + measure[Nmeasure].D = Dobs; - measure[nMeasure].imageID = 0; + measure[Nmeasure].M = J_PS1; + measure[Nmeasure].dM = 1.0 / SN; + + measure[Nmeasure].Sky = SkyCts; + measure[Nmeasure].dSky = sqrt(SkyCts); + measure[Nmeasure].photFlags = 0; + measure[Nmeasure].photFlags2 = 0; + measure[Nmeasure].airmass = 1.0; + measure[Nmeasure].az = 0.0; // irrelevant + measure[Nmeasure].Mcal = 0.0; + measure[Nmeasure].t = tzero_2mass; + measure[Nmeasure].dt = 0.0; + measure[Nmeasure].photcode = codeJ->code; + + measure[Nmeasure].averef = i; + measure[Nmeasure].objID = average[i].objID; + measure[Nmeasure].catID = average[i].catID; + + measure[Nmeasure].imageID = 0; // This is may optionally be replaced by the internal sequence (see FilterStars.c) - measure[nMeasure].detID = Nstars + 1; - Nstars ++; + measure[Nmeasure].detID = Nmeasure; + Nmeasure ++; } + + if (VERBOSE) fprintf (stderr, "added %d 2mass entries to %d stars\n", (int) Nmeasure, (int) catalog->Naverage); + + catalog->Nmeasure = Nmeasure; + return TRUE; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 37523) @@ -71,4 +71,6 @@ double Dtru = average[i].D; + stars[Nstars].starpar = starpar[nStar]; // keep this so we can save it on the generated stars + // XXX what is the epoch for the src catalog? J2000? stars[Nstars].Rref = Rtru; // add in any deviations we want here @@ -113,5 +115,5 @@ // stars[Nstars].measure.pltscale = ps1data[i].pltscale; - stars[Nstars].measure.M = Minst + ZP; + stars[Nstars].measure.M = Minst + ZP; // XX I need to compensate for the internal zero point of 25.0 stars[Nstars].measure.dM = 1.0 / SN; 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 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakeqsos.c (revision 37523) @@ -36,4 +36,6 @@ stars[i].found = FALSE; + stars[i].starpar.R = R; + stars[i].starpar.D = D; stars[i].starpar.galLon = NAN; stars[i].starpar.galLat = NAN; 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 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37523) @@ -96,4 +96,6 @@ stars[i].found = FALSE; + stars[i].starpar.R = R; + stars[i].starpar.D = D; stars[i].starpar.galLon = L; stars[i].starpar.galLat = B; 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 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c (revision 37523) @@ -8,5 +8,5 @@ double dX, dY, dR; off_t *N1, *N2; - off_t Nave, NAVE, Nmeas, NMEAS, Nmatch; + off_t Nave, NAVE, Nmeas, NMEAS, Nstarpar, NSTARPAR, Nmatch; unsigned int objID, catID; Coords tcoords; @@ -35,4 +35,5 @@ Nmatch = 0; NMEAS = Nmeas = catalog[0].Nmeasure; + NSTARPAR = Nstarpar = catalog[0].Nstarpar; // current max obj ID for this catalog @@ -221,4 +222,9 @@ REALLOCATE (catalog[0].measure, Measure, NMEAS); } + /* make sure there is space for next entry */ + if (Nstarpar >= NSTARPAR - 1) { + NSTARPAR = Nstarpar + 1000; + REALLOCATE (catalog[0].starpar, StarPar, NSTARPAR); + } if (Nave >= NAVE) { NAVE = Nave + 1000; @@ -231,4 +237,6 @@ catalog[0].average[Nave].D = stars[i].average.D; + catalog[0].average[Nave].Nstarpar = 1; + catalog[0].average[Nave].starparOffset = Nstarpar; catalog[0].average[Nave].Nmeasure = 1; catalog[0].average[Nave].measureOffset = Nmeas; @@ -260,4 +268,9 @@ } + // supply the starpar values from this detection + dvo_starpar_init (&catalog[0].starpar[Nstarpar]); + catalog[0].starpar[Nstarpar] = stars[i].starpar; + Nstarpar ++; + /* next[Nmeas] should always be -1 in this context (it is always the only measurement for the star) */ @@ -269,4 +282,5 @@ REALLOCATE (catalog[0].average, Average, Nave); REALLOCATE (catalog[0].measure, Measure, Nmeas); + REALLOCATE (catalog[0].starpar, StarPar, Nstarpar); catalog[0].sorted = FALSE; @@ -283,4 +297,5 @@ catalog[0].Naverage = Nave; catalog[0].Nmeasure = Nmeas; + catalog[0].Nstarpar = Nstarpar; catalog[0].Nsecf_mem = Nave*Nsecfilt; if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nmatch); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fake_stars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fake_stars.c (revision 37522) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fake_stars.c (revision 37523) @@ -1,5 +1,3 @@ # include "fakeastro.h" - -# define RESETTIME { gettimeofday (&startTimer, (void *) NULL); } int save_fake_stars (SkyTable *sky, Image *image, Stars *stars, int Nstars) { @@ -25,5 +23,5 @@ catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - catalog.catflags = LOAD_AVES | LOAD_SECF; + catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR; catalog.Nsecfilt = GetPhotcodeNsecfilt (); @@ -41,5 +39,4 @@ SetProtect (TRUE); - catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; dvo_catalog_update (&catalog, VERBOSE); SetProtect (FALSE);