Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37463) @@ -30,4 +30,5 @@ $(SRC)/fakeastro.$(ARCH).o \ $(SRC)/fakeastro_galaxy.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ $(SRC)/make_fakestars.$(ARCH).o \ $(SRC)/make_subset.$(ARCH).o \ @@ -37,4 +38,12 @@ $(SRC)/insert_fakestar.$(ARCH).o \ $(SRC)/gaussian.$(ARCH).o \ +$(SRC)/fakeastro_images.$(ARCH).o \ +$(SRC)/load_template_images.$(ARCH).o \ +$(SRC)/make_fake_images.$(ARCH).o \ +$(SRC)/get_image_patch.$(ARCH).o \ +$(SRC)/make_fake_stars.$(ARCH).o \ +$(SRC)/make_fake_stars_catalog.$(ARCH).o \ +$(SRC)/save_fake_stars.$(ARCH).o \ +$(SRC)/match_fake_stars.$(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 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37463) @@ -14,4 +14,11 @@ int found; // assigned to an object? } FakeAstro_Stars; + +typedef struct { + Average average; + Measure measure; + Lensing *lensing; // optionally carry out the lensing measurements + int found; +} Stars; /* used in find_matches, find_matches_refstars */ @@ -41,4 +48,8 @@ char *INPUT; +char *IMAGES_INPUT; +char *CATDIR_INPUT; +char *CATDIR_OUTPUT; + int PARALLEL; int PARALLEL_MANUAL; @@ -50,4 +61,6 @@ int FORCE; +float RADIUS; + SkyRegion UserPatch; @@ -58,6 +71,6 @@ void ConfigInit PROTO((int *argc, char **argv)); void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); -int args PROTO((int argc, char **argv)); -int args_client PROTO((int argc, char **argv)); +int args PROTO((int *argc, char **argv)); +int args_client PROTO((int *argc, char **argv)); void set_db PROTO((FITS_DB *in)); @@ -70,5 +83,5 @@ void initialize_client (int argc, char **argv); -int fakeastro_galaxy (SkyList *skylistInput, HostTable *hosts); +int fakeastro_galaxy (); FakeAstro_Stars *make_fakestars (int Nstars); int sortStars (FakeAstro_Stars *stars, int Nstars); @@ -96,2 +109,14 @@ int strextend (char *input, char *format,...); + +int fakeastro_images (); +Image *load_template_images (int *nimage); +Image *make_fake_images (Image *image, int *nfakeImage); + +SkyRegion *get_image_patch (Image *image); +Stars *make_fake_stars (SkyTable *sky, SkyRegion *patch, Image *image, int *nstars); +Stars *make_fake_stars_catalog (Stars *stars, int *nstars, SkyRegion *patch, Catalog *catalog, Image *image); +int save_fake_stars (SkyTable *sky, SkyRegion *patch, Image *image, Stars *stars, int Nstars); +int match_fake_stars (SkyRegion *region, Stars *stars, unsigned int NstarsIn, Catalog *catalog, Image *image); + +int InitStar (Stars *star); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37463) @@ -19,9 +19,13 @@ // if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; - // force CATDIR to be absolute (so parallel mode will work) - GetConfig (config, "CATDIR", "%s", 0, CATDIR); - char *tmpcatdir = abspath (CATDIR, DVO_MAX_PATH); - strcpy (CATDIR, tmpcatdir); - free (tmpcatdir); + if (FAKEASTRO_OP == OP_GALAXY) { + // force CATDIR to be absolute (so parallel mode will work) + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + char *tmpcatdir = abspath (CATDIR, DVO_MAX_PATH); + strcpy (CATDIR, tmpcatdir); + free (tmpcatdir); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + } GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); @@ -30,14 +34,31 @@ GetConfig (config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); - sprintf (ImageCat, "%s/Images.dat", CATDIR); - if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) SKY_DEPTH = 2; if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) SKY_TABLE[0] = 0; + + /* set the default search radius */ + if (!ScanConfig (config, "ADDSTAR_RADIUS", "%f", 0, &RADIUS)) { + GetConfig (config, "RADIUS", "%f", 0, &RADIUS); + } + if (RADIUS < 0.0001) { + fprintf (stderr, "absurd match radius %f\n", RADIUS); + exit (2); + } if (*CATMODE == 0) strcpy (CATMODE, "RAW"); if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + char *CATDIR_CHECK = (FAKEASTRO_OP == OP_GALAXY) ? CATDIR : CATDIR_OUTPUT; + + // 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); + } + /* update master photcode table if not defined */ - sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR_CHECK); if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile, TRUE)) { fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/StarOps.c (revision 37463) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/StarOps.c (revision 37463) @@ -0,0 +1,12 @@ +# include "fakeastro.h" + +int InitStar (Stars *star) { + + + dvo_measure_init (&star[0].measure); + dvo_average_init (&star[0].average); + star[0].found = -1; // found == -1 -> not yet found (use enums?) + + star[0].lensing = NULL; // we only populate this if needed + return TRUE; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37463) @@ -3,5 +3,5 @@ void usage_client (void); -int args (int argc, char **argv) { +int args (int *argc, char **argv) { int N; @@ -10,12 +10,12 @@ FAKEASTRO_OP = OP_NONE; - if ((N = get_argument (argc, argv, "-images"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-images"))) { + remove_argument (N, argc, argv); FAKEASTRO_OP = OP_IMAGES; } - if ((N = get_argument (argc, argv, "-galaxy"))) { + if ((N = get_argument (*argc, argv, "-galaxy"))) { if (FAKEASTRO_OP != OP_NONE) usage(); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); FAKEASTRO_OP = OP_GALAXY; } @@ -28,23 +28,23 @@ UserPatch.Dmin = -90; UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-region"))) { + remove_argument (N, argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, argc, argv); UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, argc, argv); UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } - if ((N = get_argument (argc, argv, "-catalog"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-catalog"))) { + remove_argument (N, argc, argv); UserPatch.Rmin = atof (argv[N]); UserPatch.Rmax = UserPatch.Rmin + 0.001; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); UserPatch.Dmin = atof (argv[N]); UserPatch.Dmax = UserPatch.Dmin + 0.001; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } @@ -52,20 +52,20 @@ // based on the state of the SkyTable PARALLEL = FALSE; - if ((N = get_argument (argc, argv, "-parallel"))) { + if ((N = get_argument (*argc, argv, "-parallel"))) { PARALLEL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // this is a test mode : rather than launching the remote jobs and waiting for completion, // fakeastro will simply list the remote command and wait for the user to signal completion PARALLEL_MANUAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-manual"))) { + if ((N = get_argument (*argc, argv, "-parallel-manual"))) { PARALLEL = TRUE; // -parallel-manual implies -parallel PARALLEL_MANUAL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // this is a test mode : rather than launching the fakeastro_client jobs remotely, they are // run in serial via 'system' PARALLEL_SERIAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if ((N = get_argument (*argc, argv, "-parallel-serial"))) { if (PARALLEL_MANUAL) { fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); @@ -74,53 +74,53 @@ PARALLEL = TRUE; // -parallel-serial implies -parallel PARALLEL_SERIAL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // MaxDensityUse = FALSE; - // if ((N = get_argument (argc, argv, "-max-density"))) { - // remove_argument (N, &argc, argv); + // if ((N = get_argument (*argc, argv, "-max-density"))) { + // remove_argument (N, argc, argv); // MaxDensityValue = atof(argv[N]); - // remove_argument (N, &argc, argv); + // remove_argument (N, argc, argv); // MaxDensityUse = TRUE; // } FORCE = FALSE; - if ((N = get_argument (argc, argv, "-force"))) { + if ((N = get_argument (*argc, argv, "-force"))) { FORCE = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { + if ((N = get_argument (*argc, argv, "-v"))) { VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-vv"))) { VERBOSE = VERBOSE2 = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } if (FAKEASTRO_OP == OP_IMAGES) { // mandatory arguments to fakeastro -images -input catdir -output -catdir -images images.fits - if ((N = get_argument (argc, argv, "-input"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-input"))) { + remove_argument (N, argc, argv); CATDIR_INPUT = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } else { fprintf (stderr, "missing -input (catdir)\n"); exit (2); } - if ((N = get_argument (argc, argv, "-output"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-output"))) { + remove_argument (N, argc, argv); CATDIR_OUTPUT = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } else { fprintf (stderr, "missing -output (catdir)\n"); exit (2); } - if ((N = get_argument (argc, argv, "-images"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-input-images"))) { + remove_argument (N, argc, argv); IMAGES_INPUT = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } else { fprintf (stderr, "missing -images (images)\n"); @@ -129,9 +129,8 @@ } - if (argc != 1) usage (); return TRUE; } -int args_client (int argc, char **argv) { +int args_client (int *argc, char **argv) { int N; @@ -146,42 +145,42 @@ HOST_ID = 0; - if ((N = get_argument (argc, argv, "-hostID"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-hostID"))) { + remove_argument (N, argc, argv); HOST_ID = atoi (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } if (!HOST_ID) usage_client(); HOSTDIR = NULL; - if ((N = get_argument (argc, argv, "-hostdir"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-hostdir"))) { + remove_argument (N, argc, argv); HOSTDIR = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } if (!HOSTDIR) usage_client(); CPT_FILE = NULL; - if ((N = get_argument (argc, argv, "-cpt"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-cpt"))) { + remove_argument (N, argc, argv); CPT_FILE = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } if (!CPT_FILE) usage_client(); INPUT = NULL; - if ((N = get_argument (argc, argv, "-input"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-input"))) { + remove_argument (N, argc, argv); INPUT = strcreate (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } if (!INPUT) usage_client(); - if ((N = get_argument (argc, argv, "-images"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-images"))) { + remove_argument (N, argc, argv); FAKEASTRO_OP = OP_IMAGES; } - if ((N = get_argument (argc, argv, "-galaxy"))) { + if ((N = get_argument (*argc, argv, "-galaxy"))) { if (FAKEASTRO_OP != OP_NONE) usage(); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); FAKEASTRO_OP = OP_GALAXY; } @@ -193,44 +192,43 @@ UserPatch.Dmin = -90; UserPatch.Dmax = +90; - if ((N = get_argument (argc, argv, "-region"))) { - remove_argument (N, &argc, argv); - UserPatch.Rmin = atof (argv[N]); - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-region"))) { + remove_argument (N, argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, argc, argv); UserPatch.Rmax = atof (argv[N]); - remove_argument (N, &argc, argv); - UserPatch.Dmin = atof (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, argc, argv); UserPatch.Dmax = atof (argv[N]); - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } - if ((N = get_argument (argc, argv, "-catalog"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (*argc, argv, "-catalog"))) { + remove_argument (N, argc, argv); UserPatch.Rmin = atof (argv[N]); UserPatch.Rmax = UserPatch.Rmin + 0.001; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); UserPatch.Dmin = atof (argv[N]); UserPatch.Dmax = UserPatch.Dmin + 0.001; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // MaxDensityUse = FALSE; - // if ((N = get_argument (argc, argv, "-max-density"))) { - // remove_argument (N, &argc, argv); + // if ((N = get_argument (*argc, argv, "-max-density"))) { + // remove_argument (N, argc, argv); // MaxDensityValue = atof(argv[N]); - // remove_argument (N, &argc, argv); + // remove_argument (N, argc, argv); // MaxDensityUse = TRUE; // } VERBOSE = VERBOSE2 = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { + if ((N = get_argument (*argc, argv, "-v"))) { VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-vv"))) { + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-vv"))) { VERBOSE = VERBOSE2 = TRUE; - remove_argument (N, &argc, argv); - } - - if (argc != 1) usage_client (); + remove_argument (N, argc, argv); + } + return TRUE; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37463) @@ -2,4 +2,6 @@ int main (int argc, char **argv) { + + gauss_init (50000); /* get configuration info, args */ @@ -9,9 +11,9 @@ case OP_GALAXY: /* the object analysis is a separate process iterating over catalogs */ - fakeastro_galaxy (skylist, hosts); + fakeastro_galaxy (); exit (0); case OP_IMAGES: - // fakeastro_images (skylist); + fakeastro_images (); exit (0); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c (revision 37463) @@ -3,7 +3,5 @@ int main (int argc, char **argv) { - // need to construct these options with args_loadstarpar... - ConfigInit (&argc, argv); - args_client (argc, argv); + initialize_client (argc, argv); // client is called with a pointer to the file to be loaded Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c (revision 37463) @@ -1,9 +1,11 @@ # include "fakeastro.h" -int fakeastro_galaxy (SkyList *skylistInput, HostTable *hosts) { +int fakeastro_galaxy () { + + int n, i; SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); - SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + SkyList *skylistInput = SkyListByPatch (sky, -1, &UserPatch); // load the list of hosts @@ -17,5 +19,4 @@ // ensure that the paths are absolute path names - int i; for (i = 0; i < hosts->Nhosts; i++) { char *tmppath = abspath (hosts->hosts[i].pathname, DVO_MAX_PATH); @@ -27,17 +28,4 @@ init_remote_hosts (); } - - int n, i; - - // XXX check for an existing CATDIR and exit? - - /* outline - - * generate N stars (random draw from galaxy mode) - * assign them to catalogs and distribute - - */ - - gauss_init (50000); int Nloop = 1; @@ -93,2 +81,10 @@ return TRUE; } + + /* outline + + * generate N stars (random draw from galaxy mode) + * assign them to catalogs and distribute + + */ + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_images.c (revision 37463) @@ -1,30 +1,88 @@ # include "fakeastro.h" -int fakeastro_images (HostTable *hosts) { +int fakeastro_images () { + + int i, j; + + FITS_DB db; + + /*** update the image table ***/ + snprintf (ImageCat, DVO_MAX_PATH, "%s/Images.dat", CATDIR_OUTPUT); + + /* setup image table format and lock */ + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + int status = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + /* load or create the image table */ + if (db.dbstate == LCK_EMPTY) { + if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); + dvo_image_create (&db, GetZeroPoint()); + } else { + if (!dvo_image_load (&db, VERBOSE, FALSE)) { + Shutdown ("can't read image catalog %s", db.filename); + } + } SkyTable *skyTableInput = SkyTableLoadOptimal (CATDIR_INPUT, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (skyTableInput, CATDIR_INPUT, "cpt"); - // SkyList *skyListInput = SkyListByPatch (skyInpt, -1, &UserPatch); SkyTable *skyTableOutput = SkyTableLoadOptimal (CATDIR_OUTPUT, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (skyTableOutput, CATDIR_OUTPUT, "cpt"); - // SkyList *skyListOutput = SkyListByPatch (skyInpt, -1, &UserPatch); - int Nimages; - Images *images = load_template_images (&Nimages); + int Nimage = 0; + int NIMAGE = 1000; + Image *image = NULL; - for (i = 0; i < Nimages; i++) { + ALLOCATE (image, Image, NIMAGE); + + int Nrefimage; + Image *refimage = load_template_images (&Nrefimage); + + for (i = 0; i < Nrefimage; i++) { // we only want to make fake images for the exposures - if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + if (strcmp(&refimage[i].coords.ctype[4], "-DIS")) continue; - Images *fakeImages = make_fake_images (&images[i], &NfakeImages); + int NfakeImage; + Image *fakeImage = make_fake_images (&refimage[i], &NfakeImage); - for (j = 0; j < NfakeImages; j++) { - Stars *fakeStars = make_fake_stars (skyTableInput, skyTableOutput, &fakeImages[i], &NfakeStars); + for (j = 0; j < NfakeImage; j++) { + + // we only want to make fake stars for the fake chips + if (strcmp(&fakeImage[j].coords.ctype[4], "-WRP")) continue; + + SkyRegion *patch = get_image_patch (&fakeImage[j]); + + int NfakeStars; + Stars *fakeStars = make_fake_stars (skyTableInput, patch, &fakeImage[j], &NfakeStars); - // save fake stars... + save_fake_stars (skyTableOutput, patch, &fakeImage[j], fakeStars, NfakeStars); + + free (fakeStars); + free (patch); } + + if (Nimage + NfakeImage >= NIMAGE) { + NIMAGE += 1000; + REALLOCATE (image, Image, NIMAGE); + } + + for (j = 0; j < NfakeImage; j++) { + memcpy (&image[Nimage], &fakeImage[j], sizeof(Image)); + } + + free (fakeImage); } + + /* add the new image and save */ + dvo_image_addrows (&db, image, Nimage); + SetProtect (TRUE); + dvo_image_update (&db, VERBOSE); + SetProtect (FALSE); + dvo_image_unlock (&db); /* unlock? */ + + exit (0); } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/get_image_patch.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/get_image_patch.c (revision 37463) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/get_image_patch.c (revision 37463) @@ -0,0 +1,36 @@ +# include "fakeastro.h" + +SkyRegion *get_image_patch (Image *image) { + + SkyRegion *region; + + // find the R,D coords of the 4 corners and 4 edge midpoints + + static double Xpt[] = {0.0, 0.5, 1.0, 0.0, 1.0, 0.0, 0.5, 1.0}; + static double Ypt[] = {0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0}; + + double Rmin = +480.0; + double Rmax = -360.0; + double Dmin = +90.0; + double Dmax = -90.0; + + int i; + for (i = 0; i < 8; i++) { + double R, D; + XY_to_RD (&R, &D, Xpt[i]*image->NX, Ypt[i]*image->NY, &image->coords); + + Rmin = MIN(R,Rmin); + Rmax = MAX(R,Rmax); + Dmin = MIN(D,Dmin); + Dmax = MAX(D,Dmax); + } + + ALLOCATE (region, SkyRegion, 1); + + region->Rmin = Rmin; + region->Rmax = Rmax; + region->Dmin = Dmin; + region->Dmax = Dmax; + + return region; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c (revision 37463) @@ -12,14 +12,7 @@ if (get_argument (argc, argv, "--help")) usage(); + args (&argc, argv); ConfigInit (&argc, argv); - args (argc, argv); - - // check for existence of CATDIR - struct stat filestat; - int status = stat (CATDIR, &filestat); - if (!FORCE && (status == 0)) { - fprintf (stderr, "directory %s exist, refusing to contaminate\n", CATDIR); - exit (1); - } + if (argc != 1) usage (); } @@ -27,6 +20,7 @@ void initialize_client (int argc, char **argv) { + args_client (&argc, argv); ConfigInit (&argc, argv); - args_client (argc, argv); + if (argc != 1) usage_client (); } Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/load_template_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/load_template_images.c (revision 37463) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/load_template_images.c (revision 37463) @@ -0,0 +1,34 @@ +# include "fakeastro.h" + +Image *load_template_images (int *nimage) { + + FITS_DB db; + + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + int status = dvo_image_lock (&db, IMAGES_INPUT, 60.0, LCK_SOFT); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + // load the image table + if (db.dbstate == LCK_EMPTY) { + dvo_image_unlock (&db); // unlock input + *nimage = 0; + return NULL; + } + + // this operation reads the PHU header (db.header) + if (!dvo_image_load (&db, VERBOSE, TRUE)) { + Shutdown ("can't read input image catalog %s", db.filename); + } + dvo_image_unlock (&db); // unlock input + + off_t Nimage; + Image *image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + return (NULL); + } + + *nimage = Nimage; + return image; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_images.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_images.c (revision 37463) @@ -21,4 +21,9 @@ for (iy = 0; iy < 8; iy++) { + if ((ix == 0) && (iy == 0)) continue; + if ((ix == 7) && (iy == 0)) continue; + if ((ix == 0) && (iy == 7)) continue; + if ((ix == 7) && (iy == 7)) continue; + fakeImage[N] = fakeImage[0]; @@ -30,11 +35,11 @@ fakeImage[N].coords.crpix1 = (ix - 3)*4900; fakeImage[N].coords.crpix2 = (iy - 3)*4900; - fakeImage[N].coords.pc1_1 = 1.0 - fakeImage[N].coords.pc2_2 = 1.0 + fakeImage[N].coords.pc1_1 = 1.0; + fakeImage[N].coords.pc2_2 = 1.0; } else { fakeImage[N].coords.crpix1 = (3 - ix)*4900; fakeImage[N].coords.crpix2 = (4 - iy)*4900; - fakeImage[N].coords.pc1_1 = -1.0 - fakeImage[N].coords.pc2_2 = -1.0 + fakeImage[N].coords.pc1_1 = -1.0; + fakeImage[N].coords.pc2_2 = -1.0; } @@ -53,5 +58,6 @@ fakeImage[N].NY = 4850; - fakeImage[N].photcode = 10000 + iy*10 + ix; + // XXX need a way to choose this better... + fakeImage[N].photcode = 10100 + iy*10 + ix; fakeImage[N].exptime = 45.0; N++; Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars.c (revision 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars.c (revision 37463) @@ -1,60 +1,37 @@ # include "fakeastro.h" -Stars *make_fake_stars (SkyTable *skyTableInput, SkyTable *skyTableOutput, Image *image, int *nfakeStars) { +Stars *make_fake_stars (SkyTable *sky, SkyRegion *patch, Image *image, int *nstars) { - // find the R,D coords of the 4 corners and 4 edge midpoints - - static double Xpt[] = {0, 2425, 4850, 0, 4850, 0, 2425, 4850}; - static double Ypt[] = {0, 0, 0, 2425, 2425, 4850, 4850, 4850}; - - double Rmin = +480.0; - double Rmax = -360.0; - double Dmin = +90.0; - double Dmax = -90.0; - - for (i = 0; i < 8; i++) { - double R, D; - XY_to_RD (&R, &D, Xpt[i], Ypt[i], image->coords); - - Rmin = MIN(R,Rmin); - Rmax = MAX(R,Rmax); - Dmin = MIN(D,Dmin); - Dmax = MAX(D,Dmax); - - } - - SkyList *skyInput = SkyListByBounds (skyTableInput, -1, Rmin, Rmax, Dmin, Dmax); + SkyList *skylist = SkyListByPatch (sky, -1, patch); int Nstars = 0; + Stars *stars = NULL; + + Catalog catalog; // load stars from database in these regions - for (i = 0; i < skyInput->Nregions; i++) { + int i; + for (i = 0; i < skylist->Nregions; i++) { // set the parameters which guide catalog open/load/create - catInput.filename = skyInput[0].filename[i]; - catInput.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data - catInput.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - catInput.catflags = LOAD_AVES | LOAD_SECF; - catInput.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catInput, skylist[0].regions[i], VERBOSE, "r")) { + catalog.filename = skylist[0].filename[i]; + 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_STARPAR; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) { continue; } // generate fake measurements for this image - stars = make_fake_stars_catalog (catInput, stars, &Nstars); + stars = make_fake_stars_catalog (stars, &Nstars, patch, &catalog, image); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); } - // will these match, or should I mangle the input cpt names? - SkyList *skyOutput = SkyListByBounds (skyTableOutput, -1, Rmin, Rmax, Dmin, Dmax); + SkyListFree (skylist); - // load stars from database in these regions - for (i = 0; i < skyOutput->Nregions; i++) { - if (!dvo_catalog_open (&catOutput, skyOutput[0].regions[i], VERBOSE, "w")) { - fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); - exit (2); - } - - save_fake_stars (stars, Nstars, catOutput); - } + *nstars = Nstars; + return stars; } 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 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c (revision 37463) @@ -1,3 +1,7 @@ # include "fakeastro.h" +# define SCALE 0.001 + +// sky[Nsecfilt], crude and hard-wired for now +float sky[] = {100.0, 150.0, 200.0, 300.0, 400.0}; // things to figure out: @@ -6,24 +10,53 @@ // * what about QSOs? -Stars *make_fake_stars_catalog (Catalog *catalog, Image *image, int *nfakeStars) { +Stars *make_fake_stars_catalog (Stars *stars, int *nstars, SkyRegion *region, Catalog *catalog, Image *image) { + + int Nstars = *nstars; + *nstars += catalog->Naverage; + + time_t timeRef = ohana_date_to_sec("2011/05/11,00:00:00"); + + if (!stars) { + ALLOCATE (stars, Stars, *nstars); + } else { + REALLOCATE (stars, Stars, *nstars); + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + // use photcode to get zero point + PhotCode *code = GetPhotcodebyCode (image->photcode); + int Nsec = GetPhotcodeNsec (code->equiv); + myAssert (Nsec >= 0, "undefined Nsec?"); float Mtime = 2.5*log10(image->exptime); - float ZP = image->Mcal + zeropt; + + // XXX fix this!! + double plateScale = 0.257; + + // XXX put in airmass? + float ZP = SCALE*code->C - image->Mcal + Mtime; + // float ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; // generate a set of measurements for each star entry Average *average = catalog->average; + SecFilt *secfilt = catalog->secfilt; StarPar *starpar = catalog->starpar; + int i; for (i = 0; i < catalog->Naverage; i++) { - int nStar = average[i].starparOffset + // true position from src catalog + if (!IN_REGION(average[i].R, average[i].D)) continue; - InitStar (&stars[i]); + int nStar = average[i].starparOffset; + + InitStar (&stars[Nstars]); // which filter? - double Minst = secfilt[i*Nsecfilt + Ns].M - Mtime - ZP; + double Minst = secfilt[i*Nsecfilt + Nsec].M - ZP; double Counts = pow(10.0, -0.4*Minst); - double SkyCts = Something; + double SkyCts = sky[Nsec]; double SN = Counts / sqrt(SkyCts + Counts); @@ -36,12 +69,13 @@ // * proper motion // * gaussian scatter (~ seeing) - double uR = starpar[nStar].uR; - double uD = starpar[nStar].uD; + double uR = starpar[nStar].uRA / 1000.0; // starpar are (currently) in mas / year + double uD = starpar[nStar].uDEC / 1000.0; - double Toffset = average - image.tzero; + // tzero, timeRef are in UNIX seconds, Toffset should be in years + double Toffset = (image->tzero - timeRef) / 365.25 / 86400.0; // uR,uD in linear (arcsec / yr) - double dRoff = uR*Toffset; - double dDoff = uD*Toffset; + double dRpm = uR*Toffset; + double dDpm = uD*Toffset; // uR,uD in linear arcsec @@ -49,35 +83,41 @@ double dDsee = rnd_gauss (0.0, 1.0 / SN); + double dRoff = (dRpm + dRsee) / 3600.0; + double dDoff = (dDpm + dDsee) / 3600.0; + double Robs = Rtru + dRoff / cos(Dtru*DEG_RAD); double Dobs = Dtru + dDoff; double X, Y; - RD_to_XY (&X, &Y, Robs, Dobs, image->coords); + RD_to_XY (&X, &Y, Robs, Dobs, &image->coords); + if (X < 0) continue; + if (Y < 0) continue; + if (X > image->NX) continue; + if (Y > image->NY) continue; - stars[i].measure.Xccd = X; - stars[i].measure.Yccd = Y; - stars[i].measure.dXccd = 1.0 / SN / plateScale; - stars[i].measure.dYccd = 1.0 / SN / plateScale; + stars[Nstars].measure.Xccd = X; + stars[Nstars].measure.Yccd = Y; + stars[Nstars].measure.dXccd = 1.0 / SN / plateScale; + stars[Nstars].measure.dYccd = 1.0 / SN / plateScale; - // stars[i].measure.posangle = ToShortDegrees(ps1data[i].posangle); - // stars[i].measure.pltscale = ps1data[i].pltscale; + // stars[Nstars].measure.posangle = ToShortDegrees(ps1data[i].posangle); + // stars[Nstars].measure.pltscale = ps1data[i].pltscale; - if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = Minst + ZeroPt; - } - stars[i].measure.dM = 1.0 / SN; + stars[Nstars].measure.M = Minst + ZP; + stars[Nstars].measure.dM = 1.0 / SN; - // stars[i].measure.dMcal = ps1data[i].dMcal; - stars[i].measure.Sky = X?; - stars[i].measure.dSky = X?; + // stars[Nstars].measure.dMcal = ps1data[i].dMcal; + stars[Nstars].measure.Sky = sky[Nsec]; + stars[Nstars].measure.dSky = sqrt(sky[Nsec]); - stars[i].measure.photFlags = ps1data[i].flags; - stars[i].measure.photFlags2 = ps1data[i].flags2; + stars[Nstars].measure.photFlags = 0; + stars[Nstars].measure.photFlags2 = 0; // this is may optionally be replaced by the internal sequence (see FilterStars.c) - stars[i].measure.detID = ps1data[i].detID; + stars[Nstars].measure.detID = Nstars + 1; + + Nstars ++; } + return stars; } 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 37462) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37463) @@ -35,7 +35,24 @@ */ - double z = rnd_gauss (0.0, Z_GAL); - double r = sqrt(drand48()) * R_GAL; - double Lrad = drand48() * 2 * M_PI; + if (i % 100000 == 0) fprintf (stderr, "."); + double z,r,L,B,R,D,Lrad,Brad; + + int inPatch = FALSE; + while (!inPatch) { + z = rnd_gauss (0.0, Z_GAL); + r = sqrt(drand48()) * R_GAL; + Lrad = drand48() * 2 * M_PI; + Brad = atan2(z,r); + + 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; + } // double x = r*cos(L); @@ -43,5 +60,4 @@ double distance = sqrt (SQ(r) + SQ(z)); - double Brad = atan2(z,r); double uL_gal = (A_oort * cos(2.0*Lrad) + B_oort) * cos(Brad) * iFkap; @@ -59,10 +75,4 @@ double Mr = rnd_gauss (11.25, 1.0); - double L = Lrad*DEG_RAD; - double B = Brad*DEG_RAD; - - double R, D; - ApplyTransform (&R, &D, L, B, transform); - // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf double Rrad = R*RAD_DEG; @@ -104,4 +114,5 @@ stars[i].starpar.uDEC = uD; } + fprintf (stderr, "\n"); return stars; } 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 37463) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c (revision 37463) @@ -0,0 +1,296 @@ +# include "fakeastro.h" + +int match_fake_stars (SkyRegion *region, Stars *stars, unsigned int NstarsIn, Catalog *catalog, Image *image) { + + off_t i, j, n, N, J, Jmin, status, Nstars; + double RADIUS2, Rmin; + double *X1, *Y1, *X2, *Y2; + double dX, dY, dR; + off_t *N1, *N2; + off_t Nave, NAVE, Nmeas, NMEAS, Nmatch; + unsigned int objID, catID; + Coords tcoords; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + // which equiv photcode are we? + PhotCode *code = GetPhotcodebyCode (image->photcode); + int Nsec = GetPhotcodeNsec (code->equiv); + myAssert (Nsec >= 0, "undefined Nsec?"); + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, double, NstarsIn); + ALLOCATE (Y1, double, NstarsIn); + ALLOCATE (N1, off_t, NstarsIn); + + /** allocate local arrays (catalog) **/ + NAVE = Nave = catalog[0].Naverage; + ALLOCATE (X2, double, NAVE); + ALLOCATE (Y2, double, NAVE); + ALLOCATE (N2, off_t, NAVE); + ALLOCATE (catalog[0].found_t, off_t, NAVE); + /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */ + + /* internal counters */ + Nmatch = 0; + NMEAS = Nmeas = catalog[0].Nmeasure; + + // current max obj ID for this catalog + objID = catalog[0].objID; + catID = catalog[0].catID; + + /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has + * the advantage that every point in R,D has a mapping to a unique X,Y. However, note + * that not all possible X,Y points map back to R,D and the local plate scale changes + * far from the projection pole. We use the center of the region (catalog) for crval1,2. + */ + InitCoords (&tcoords, "DEC--ARC"); + tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); + if (region[0].Dmax < 90) { + tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax); + } else { + tcoords.crval2 = 90.0; + } + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + + /* build spatial index (RA sort) referencing input array sequence */ + Nstars = 0; + for (i = 0; i < NstarsIn; i++) { + status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].average.R, stars[i].average.D, &tcoords); + if (!status) continue; + N1[Nstars] = i; + Nstars ++; + } + if (Nstars < 1) { + if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename); + free (X1); + free (Y1); + free (N1); + free (X2); + free (Y2); + free (N2); + return (0); + } + if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars); + + /* build spatial index (RA sort) */ + for (i = 0; i < Nave; i++) { + RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + catalog[0].found_t[N2[i]] = -1; + } + if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave); + + /* RADIUS is global from ConfigInit.c */ + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + if (!finite(X1[i]) || !finite(Y1[i])) { + i++; + continue; + } + if (!finite(X2[j]) || !finite(Y2[j])) { + j++; + continue; + } + + /* negative dX: j is too large */ + dX = X1[i] - X2[j]; + if (dX <= -1.02*RADIUS) { + i++; + continue; + } + /* positive dX, i is too large */ + if (dX >= 1.02*RADIUS) { + j++; + continue; + } + + if (stars[N1[i]].found != -1) { + /* this star has already been assigned to an object in this or another catalog */ + i++; + continue; + } + + /* this block will match a given detection to the closest object within range of that detection. + XXX note that this matches ALL detections within range of the single object to that same object + this is bad, but I cannot just go in linear order (ie, mark off each object as they are + used). I should make a list of all Nobj * Ndet pairs in range and choose the matches + based on their separations. UGH + */ + + /* within match range; look for matches */ + Jmin = -1; + Rmin = RADIUS2; + for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) { + /* find closest match for this detection */ + dX = X1[i] - X2[J]; + dY = Y1[i] - Y2[J]; + dR = dX*dX + dY*dY; + if (dR > RADIUS2) continue; + if (dR < Rmin) { + Rmin = dR; + Jmin = J; + } + } + + /* no match, try next detection */ + if (Jmin == -1) { + i++; + continue; + } + + /*** a match is found, add to average, measure ***/ + Nmatch ++; + n = N2[Jmin]; + N = N1[i]; + + /* make sure there is space for next entry */ + if (Nmeas >= NMEAS) { + NMEAS = Nmeas + 1000; + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + + /** add measurements for this star **/ + + // set the new measurements + catalog[0].measure[Nmeas] = stars[N].measure; + + // measure now carries R,D (not dR,dD) + // note that ReadStarsFITS does not set measure.R,D and average.R,D + catalog[0].measure[Nmeas].R = stars[N].average.R; + catalog[0].measure[Nmeas].D = stars[N].average.D; + catalog[0].measure[Nmeas].dbFlags = 0; + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].objID = catalog[0].average[n].objID; + catalog[0].measure[Nmeas].catID = catalog[0].catID; + + float dRoff = dvoOffsetR(&catalog[0].measure[Nmeas], &catalog[0].average[n]); + + // rationalize dR + if (dRoff > +180.0*3600.0) { + // average on high end of boundary, move star up + catalog[0].measure[Nmeas].R += 360.0; + dRoff -= 360.0*3600.0; + } + if (dRoff < -180.0*3600.0) { + // average on low end of boundary, move star down + catalog[0].measure[Nmeas].R -= 360.0; + dRoff += 360.0*3600.0; + } + if (fabs(dRoff) > 10*RADIUS) { + // take declination into account and check again. + double cosD = cos(RAD_DEG*catalog[0].average[n].D); + if (fabs(dRoff*cosD) > 10*RADIUS) { + fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", + catalog[0].average[n].R, catalog[0].average[n].D, + stars[N].average.R, stars[N].average.D, + X1[i], X2[Jmin], + Y1[i], Y2[Jmin]); + // XXX abort on this? -- this is a bad failure... + } + } + + /* set the average magnitude if not already set and the photcode.equiv is not 0 */ + /* in UPDATE mode, this value is not saved; use relphot to recalculate */ + if (Nsec > -1) { + if (isnan(catalog[0].secfilt[n*Nsecfilt+Nsec].M)) { + catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF); + } + } + + /* Nm is updated, but not written out in -update mode (for existing entries) + Nm is recalculated in build_meas_links if loaded table is not sorted */ + stars[N].found = Nmeas; + catalog[0].found_t[n] = Nmeas; + catalog[0].average[n].Nmeasure ++; + Nmeas ++; + i++; + } + + /** incorporate unmatched image stars, if this star is in field of this catalog **/ + /* these new entries are all written out in UPDATE mode */ + for (i = 0; i < Nstars; i ++) { + if (stars[i].found != -1) continue; + if (!IN_REGION (stars[i].average.R, stars[i].average.D)) continue; + + /* make sure there is space for next entry */ + if (Nmeas >= NMEAS - 1) { + NMEAS = Nmeas + 1000; + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + if (Nave >= NAVE) { + NAVE = Nave + 1000; + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); + } + + dvo_average_init (&catalog[0].average[Nave]); + catalog[0].average[Nave].R = stars[i].average.R; + catalog[0].average[Nave].D = stars[i].average.D; + + catalog[0].average[Nave].Nmeasure = 1; + catalog[0].average[Nave].measureOffset = Nmeas; + catalog[0].average[Nave].objID = objID; + catalog[0].average[Nave].catID = catID; + + objID ++; + + for (j = 0; j < Nsecfilt; j++) { + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); + } + + // supply the measurments from this detection + dvo_measure_init (&catalog[0].measure[Nmeas]); + catalog[0].measure[Nmeas] = stars[i].measure; + + // the following measure elements cannot be set until here: + catalog[0].measure[Nmeas].R = stars[i].average.R; + catalog[0].measure[Nmeas].D = stars[i].average.D; + catalog[0].measure[Nmeas].dbFlags = 0; + catalog[0].measure[Nmeas].averef = Nave; + catalog[0].measure[Nmeas].objID = catalog[0].average[Nave].objID; + catalog[0].measure[Nmeas].catID = catalog[0].catID; + + /* set the average magnitude if not already set and the photcode.equiv is not 0 */ + /* XXX NOTE : this value is not saved; use relphot to recalculate */ + if (Nsec > -1) { + catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas], MAG_CLASS_PSF); + } + + /* next[Nmeas] should always be -1 in this context (it is always the only + measurement for the star) */ + stars[i].found = Nmeas; + Nmeas ++; + Nave ++; + } + + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + catalog[0].sorted = FALSE; + + /* note stars which have been found in this catalog */ + for (i = 0; i < NstarsIn; i++) { + if (stars[i].found > -1) { + stars[i].found = -2; + } + } + + /* check if the catalog has changed? if no change, no need to write */ + catalog[0].objID = objID; // new max value, save on catalog close + catalog[0].Naverage = Nave; + catalog[0].Nmeasure = Nmeas; + 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); + + free (X1); + free (Y1); + free (N1); + free (N2); + free (X2); + free (Y2); + + return (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 37463) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fake_stars.c (revision 37463) @@ -0,0 +1,51 @@ +# include "fakeastro.h" + +# define RESETTIME { gettimeofday (&startTimer, (void *) NULL); } + +int save_fake_stars (SkyTable *sky, SkyRegion *patch, Image *image, Stars *stars, int Nstars) { + + Catalog catalog; + + // will these match, or should I mangle the input cpt names? + SkyList *skylist = SkyListByPatch (sky, -1, patch); + + int Naverage = 0; + int Nmeasure = 0; + + INITTIME; + + // load stars from database in these regions + int i; + for (i = 0; i < skylist->Nregions; i++) { + + // set the parameters which guide catalog open/load/create + catalog.filename = skylist[0].filename[i]; + 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.Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + match_fake_stars (skylist[0].regions[i], stars, Nstars, &catalog, image); + + /* report total updated values */ + Naverage += catalog.Naverage; + Nmeasure += catalog.Nmeasure; + + SetProtect (TRUE); + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING; + dvo_catalog_update (&catalog, VERBOSE); + SetProtect (FALSE); + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; + } + + SkyListFree (skylist); + return TRUE; +}