Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 7080) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 8390) @@ -34,4 +34,11 @@ 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; + } + GetConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 7080) +++ trunk/Ohana/src/relastro/src/args.c (revision 8390) @@ -1,7 +1,9 @@ # include "relphot.h" +void usage (void); int args (int argc, char **argv) { int N; + double trange; /* define time */ @@ -14,10 +16,39 @@ } remove_argument (N, &argc, argv); - if (!str_to_time (argv[N], &TSTOP)) { - fprintf (stderr, "ERROR: syntax error\n"); - return (FALSE); + if (!str_to_dtime (argv[N], &trange)) { + if (!str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } } remove_argument (N, &argc, argv); TimeSelect = TRUE; + } + + /* specify portion of the sky */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + UserPatchSelect = FALSE; + 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); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatchSelect = TRUE; } @@ -142,25 +173,29 @@ } - if (argc != 3) { - fprintf (stderr, "ERROR: USAGE: nrphot (region) (photcode)\n"); - fprintf (stderr, " options: \n"); - fprintf (stderr, " -time (start) (stop)\n"); - fprintf (stderr, " -v\n"); - fprintf (stderr, " -plot\n"); - fprintf (stderr, " -plotdelay (seconds)\n"); - fprintf (stderr, " -statmode (mode)\n"); - fprintf (stderr, " -n (nloop)\n"); - fprintf (stderr, " -reset\n"); - fprintf (stderr, " -update\n"); - fprintf (stderr, " -params\n"); - fprintf (stderr, " -mosaic (mosaic)\n"); - fprintf (stderr, " -imfreeze\n"); - fprintf (stderr, " -grid\n"); - fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); - fprintf (stderr, " -instmag min max\n"); - fprintf (stderr, " \n"); - exit (2); - } - return (TRUE); + if ( UserPatchSelect && (argc != 2)) usage (); + if (!UserPatchSelect && (argc != 3)) usage (); + + return TRUE; } +void usage () { + fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n"); + fprintf (stderr, " or: relphot (photcode) -region RA RA DEC DEC\n"); + fprintf (stderr, " options: \n"); + fprintf (stderr, " -time (start) (stop)\n"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " -plot\n"); + fprintf (stderr, " -plotdelay (seconds)\n"); + fprintf (stderr, " -statmode (mode)\n"); + fprintf (stderr, " -n (nloop)\n"); + fprintf (stderr, " -reset\n"); + fprintf (stderr, " -update\n"); + fprintf (stderr, " -params\n"); + fprintf (stderr, " -mosaic (mosaic)\n"); + fprintf (stderr, " -imfreeze\n"); + fprintf (stderr, " -grid\n"); + fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); + fprintf (stderr, " -instmag min max\n"); + fprintf (stderr, " \n"); + exit (2); +} Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 7080) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 8390) @@ -1,5 +1,5 @@ # include "relphot.h" -int bcatalog (Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion) { +int bcatalog (Catalog *subcatalog, Catalog *catalog) { int i, j, offset, ecode; @@ -18,8 +18,11 @@ for (i = 0; i < catalog[0].Naverage; i++) { if (catalog[0].average[i].Nm < 2) continue; + + /* XXX this limitation is absurd if (catalog[0].average[i].R < fullregion[0].RA[0]) continue; if (catalog[0].average[i].R > fullregion[0].RA[1]) continue; if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue; if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue; + */ /* start with all stars good */ Index: trunk/Ohana/src/relastro/src/check_permissions.c =================================================================== --- trunk/Ohana/src/relastro/src/check_permissions.c (revision 7080) +++ (revision ) @@ -1,67 +1,0 @@ -# include "relphot.h" - -void check_permissions (char *basefile) { - - char *c, dir[256], filename[256]; - struct stat filestat; - uid_t uid; - gid_t gid; - int status, cmode; - - uid = getuid(); - gid = getgid(); - - /* check permission to write to directory */ - sprintf (filename, "%s", basefile); - c = strrchr (filename, '/'); - if (c == (char *) NULL) { - strcpy (dir, "."); - } else { - *c = 0; - strcpy (dir, filename); - } - status = stat (dir, &filestat); - if (status == -1) { - fprintf (stderr, "directory %s does not exist, creating...\n", dir); - cmode = S_IRWXU | S_IRWXG | S_IRWXO; - status = mkdir (dir, cmode); - if (status == -1) { - fprintf (stderr, "ERROR: can't create %s\n", dir); - exit (1); - } - } - status = stat (dir, &filestat); - if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRWXU)) || - ((gid == filestat.st_gid) && (filestat.st_mode & S_IRWXG)) || - (filestat.st_mode & S_IRWXO)) { - } else { - fprintf (stderr, "ERROR: can't write to %s\n", dir); - exit (1); - } - - /* check permission to write to file */ - sprintf (filename, "%s", basefile); - status = stat (filename, &filestat); - if (status == 0) { /* file exists, are permissions OK? */ - if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || - ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || - ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { - } else { - fprintf (stderr, "ERROR: can't write to %s\n", filename); - exit (1); - } - } - - /* check permission to write to backup file */ - sprintf (filename, "%s~", basefile); - status = stat (filename, &filestat); - if (status == 0) { /* file exists, are permissions OK? */ - if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || - ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || - ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { - } else { - fprintf (stderr, "ERROR: can't write to %s\n", filename); - exit (1); - } - } -} Index: trunk/Ohana/src/relastro/src/find_images.c =================================================================== --- trunk/Ohana/src/relastro/src/find_images.c (revision 7080) +++ (revision ) @@ -1,207 +1,0 @@ -# include "relphot.h" - -/* this function returns a list of all images which overlap the given - set of region files. All images in the image catalog are tested - once, so there is no check that an image already has been included. - LineNum stores the locations in the Image database of the list of - images */ - -Image *find_images (FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum) { - - Image *timage, *image; - int i, j, k, m, found, nimage, Ntimage, NIMAGE; - int InRange, ecode; - double Xc[5], Yc[5], Xi[5], Yi[5], r, d, dx, dy; - int *line_number; - Coords tcoords; - - if (VERBOSE) fprintf (stderr, "finding images\n"); - - /* we make positional comparisons in the projection of catalog */ - tcoords.crval1 = 0.5*(region[0].RA[0] + region[0].RA[1]); - tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]); - tcoords.crpix1 = tcoords.crpix2 = 0.0; - tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - strcpy (tcoords.ctype, "RA---TAN"); - - timage = gfits_table_get_Image (&db[0].ftable, &Ntimage, &db[0].swapped); - - nimage = 0; - NIMAGE = 100; - ALLOCATE (image, Image, NIMAGE); - ALLOCATE (line_number, int, NIMAGE); - - for (i = 0; i < Ntimage; i++) { - - /* select images by photcode */ - ecode = GetPhotcodeEquivCodebyCode (timage[i].source); - if (ecode != photcode[0].code) continue; - - /* select images by time */ - if (TimeSelect) { - if (timage[i].tzero < TSTART) continue; - if (timage[i].tzero > TSTOP) continue; - } - - /* ignore image set to SKIP - if (timage[i].code & ID_IMAGE_SKIP) continue; */ - - /* define image corners */ - Xi[0] = 0; Yi[0] = 0; - Xi[1] = timage[i].NX; Yi[1] = 0; - Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; - Xi[3] = 0; Yi[3] = timage[i].NY; - Xi[4] = 0; Yi[4] = 0; - found = FALSE; - /* transform to tcoords */ - for (j = 0; j < 5; j++) { - XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords); - InRange = RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords); - if (!InRange) goto imskip; - } - /* compare with each region file */ - for (m = 0; (m < Nregion) && !found; m++) { - /* define catalog corners */ - Xc[0] = region[m].RA[0]; Yc[0] = region[m].DEC[0]; - Xc[1] = region[m].RA[1]; Yc[1] = region[m].DEC[0]; - Xc[2] = region[m].RA[1]; Yc[2] = region[m].DEC[1]; - Xc[3] = region[m].RA[0]; Yc[3] = region[m].DEC[1]; - Xc[4] = region[m].RA[0]; Yc[4] = region[m].DEC[0]; - for (j = 0; j < 5; j++) { - r = Xc[j]; d = Yc[j]; - RD_to_XY (&Xc[j], &Yc[j], r, d, &tcoords); - } - dx = 0.02*(Xc[2] - Xc[0]); - dy = 0.02*(Yc[2] - Yc[0]); - Xc[0] -= dx; Yc[0] -= dy; - Xc[1] += dx; Yc[1] -= dy; - Xc[2] += dx; Yc[2] += dy; - Xc[3] -= dx; Yc[3] += dy; - Xc[4] -= dx; Yc[4] -= dy; - - /* check if image corner inside catalog */ - for (j = 0; (j < 4) && !found; j++) { - found |= corner_check (&Xi[j], &Yi[j], &Xc[0], &Yc[0]); - } - /* check if catalog corner inside image */ - for (j = 0; (j < 4) && !found; j++) { - found |= corner_check (&Xc[j], &Yc[j], &Xi[0], &Yi[0]); - } - /* check if edges cross */ - for (j = 0; (j < 4) && !found; j++) { - for (k = 0; (k < 4) && !found; k++) { - found |= edge_check (&Xi[j], &Yi[j], &Xc[k], &Yc[k]); - } - } - if (found) { - image[nimage] = timage[i]; - /* always allow a 'few' image to succeed, if possible */ - if (image[nimage].code & ID_IMAGE_FEW) { - image[nimage].code &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR); - } - if (RESET) { - assignMcal (&image[nimage], (double *) NULL, -1); - image[nimage].dMcal_PS = NO_MAG; - image[nimage].code &= ~ID_IMAGE_POOR; - } - line_number[nimage] = i; - nimage ++; - if (nimage == NIMAGE) { - NIMAGE += 100; - REALLOCATE (image, Image, NIMAGE); - REALLOCATE (line_number, int, NIMAGE); - } - } - } - imskip: - continue; - } - - if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); - - REALLOCATE (image, Image, MAX (nimage, 1)); - REALLOCATE (line_number, int, MAX (nimage, 1)); - - *Nimage = nimage; - *LineNum = line_number; - return (image); -} - -/* check if line between points 0 and 1 of x1 - crosses line between points 0 and 1 of x2 */ -int edge_check (double *x1, double *y1, double *x2, double *y2) { - - double theta1, theta2; - double Theta1, Theta2; - - theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); - theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); - - if (theta1*theta2 < 0.0) { - return (FALSE); - } - - if (fabs(theta1) < fabs(theta2)) { - return (FALSE); - } - - Theta1 = theta1; - Theta2 = theta2; - theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); - theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); - - - if (theta1*theta2 < 0.0) { - return (FALSE); - } - - if (fabs(theta1) < fabs(theta2)) { - return (FALSE); - } - - return (TRUE); - -} - -/* check if point x1,y1 is in box formed by x2[0-4] */ -int corner_check (double *x1, double *y1, double *x2, double *y2) { - - int i; - double theta; - - theta = 0; - - for (i = 0; i < 4; i++) { - theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); - } - if (fabs(theta) > 6) { - return (TRUE); - } else { - return (FALSE); - } -} - -/* returns the opening angle between the three points (2 is in middle) - in range -pi to pi */ - -double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) { - - double dx1, dy1, dx2, dy2, ct, st, theta; - - dx1 = x1 - x2; - dy1 = y1 - y2; - - dx2 = x3 - x2; - dy2 = y3 - y2; - - ct = (dx1*dx2 + dy1*dy2); - st = (dx1*dy2 - dx2*dy1); - - theta = atan2 (st, ct); - - return (theta); - -} - Index: trunk/Ohana/src/relastro/src/find_regions.c =================================================================== --- trunk/Ohana/src/relastro/src/find_regions.c (revision 7080) +++ (revision ) @@ -1,102 +1,0 @@ -# include "relphot.h" - -/* given a list of images, find all region files which overlap these images. - if a region file does not exist, we should recognize this fact somehow */ - -GSCRegion *find_regions (Image *image, int Nimage, int *Nregions, GSCRegion *fullregion) { - - int i, nregions; - double x, y; - GSCRegion *regions; - double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec; - - if (VERBOSE) fprintf (stderr, "finding regions: "); - - /* find RA and DEC range of images */ - MinRA = MinDEC = +361; - MaxRA = MaxDEC = -361; - for (i = 0; i < Nimage; i++) { - for (x = 0.0; x <= 1.0; x+=0.25) { - for (y = 0.0; y <= 1.0; y+=0.25) { - XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); - MinRA = MIN (ra, MinRA); - MaxRA = MAX (ra, MaxRA); - MinDEC = MIN (dec, MinDEC); - MaxDEC = MAX (dec, MaxDEC); - } - } - } - - /* if images overlap 0,360 boundary, then use -180 to +180 range */ - if (MaxRA - MinRA > 180) { - MinRA = +361; - MaxRA = -361; - for (i = 0; i < Nimage; i++) { - for (x = 0.0; x <= 1.0; x+=0.25) { - for (y = 0.0; y <= 1.0; y+=0.25) { - XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); - if (ra > 180) ra -= 360; - MinRA = MIN (ra, MinRA); - MaxRA = MAX (ra, MaxRA); - } - } - } - } - - fullregion[0].RA[0] = MinRA; - fullregion[0].RA[1] = MaxRA; - fullregion[0].DEC[0] = MinDEC; - fullregion[0].DEC[1] = MaxDEC; - - /* find all region files containing specified ra, dec range */ - regions = get_regions (MinRA, MaxRA, MinDEC, MaxDEC, &nregions); - - if (VERBOSE) fprintf (stderr, "%d found\n", nregions); - - *Nregions = nregions; - return (regions); -} - - -void getfullregion (Image *image, int Nimage, GSCRegion *fullregion) { - - int i; - double x, y; - double MinRA, MaxRA, MinDEC, MaxDEC, ra, dec; - - /* find RA and DEC range of images */ - MinRA = MinDEC = +361; - MaxRA = MaxDEC = -361; - for (i = 0; i < Nimage; i++) { - for (x = 0.0; x <= 1.0; x+=0.25) { - for (y = 0.0; y <= 1.0; y+=0.25) { - XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); - MinRA = MIN (ra, MinRA); - MaxRA = MAX (ra, MaxRA); - MinDEC = MIN (dec, MinDEC); - MaxDEC = MAX (dec, MaxDEC); - } - } - } - - /* if images overlap 0,360 boundary, then use -180 to +180 range */ - if (MaxRA - MinRA > 180) { - MinRA = +361; - MaxRA = -361; - for (i = 0; i < Nimage; i++) { - for (x = 0.0; x <= 1.0; x+=0.25) { - for (y = 0.0; y <= 1.0; y+=0.25) { - XY_to_RD (&ra, &dec, image[i].NX*(1.1*x - 0.05), image[i].NY*(1.1*y - 0.05), &image[i].coords); - if (ra > 180) ra -= 360; - MinRA = MIN (ra, MinRA); - MaxRA = MAX (ra, MaxRA); - } - } - } - } - - fullregion[0].RA[0] = MinRA; - fullregion[0].RA[1] = MaxRA; - fullregion[0].DEC[0] = MinDEC; - fullregion[0].DEC[1] = MaxDEC; -} Index: trunk/Ohana/src/relastro/src/free_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/free_catalogs.c (revision 7080) +++ (revision ) @@ -1,13 +1,0 @@ -# include "relphot.h" - -void free_catalogs (Catalog *catalog, int Ncatalog) { - - int i; - - for (i = 0; i < Ncatalog; i++) { - free (catalog[i].secfilt); - free (catalog[i].average); - free (catalog[i].measure); - } - -} Index: trunk/Ohana/src/relastro/src/gcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/gcatalog.c (revision 7080) +++ (revision ) @@ -1,26 +1,0 @@ -# include "relphot.h" - -int gcatalog (Catalog *catalog, int FINAL) { - - if (FINAL) { - catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; - } else { - catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; - } - - /* CATMODE and CATFORMAT determined from catalog */ - - if (!load_catalog (catalog, VERBOSE)) { - fprintf (stderr, "ERROR: failure loading catalog %s\n", catalog[0].filename); - exit (1); - } - - if (VERBOSE) { - fprintf (stderr, "read %d stars from catalog file %s (%d measurements)\n", - catalog[0].Naverage, catalog[0].filename, catalog[0].Nmeasure); - } - return (TRUE); -} - - - Index: trunk/Ohana/src/relastro/src/get_regions.c =================================================================== --- trunk/Ohana/src/relastro/src/get_regions.c (revision 7080) +++ (revision ) @@ -1,167 +1,0 @@ -# include "relphot.h" -double BigDecBounds[] = {0.0, 7.5, 15.0, 22.5, 30.0, 37.5, 45.0, - 52.5, 60.0, 67.5, 75.0, 82.5, 90.0, - 0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, - -52.5, -60.0, -67.5, -75.0, -82.5, -90.0}; - -double Areas[] = {2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0, - 2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0}; - -char *DecSections[] = {"N0000", "N0730", "N1500", "N2230", "N3000", "N3730", "N4500", - "N5230", "N6000", "N6730", "N7500", "N8230", "weirdness", - "S0000", "S0730", "S1500", "S2230", "S3000", "S3730", "S4500", - "S5230", "S6000", "S6730", "S7500", "S8230", "weirdness"}; - -char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", - "n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", - "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", - "s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"}; - -char *disk[] = {"disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", - "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "weirdness", - "disk 1", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", - "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "weirdness"}; - -int NBigRASections [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3, 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3}; - -int NDecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, - 0, 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48}; - -/* returns a list of region files within the desired RA, DEC region */ -GSCRegion *get_regions (double minRa, double maxRa, double minDec, double maxDec, int *Nregions) { - - char buffer[28800], temp[50]; - GSCRegion *regions; - FILE *f; - double RA0, RA1, DEC0, DEC1; - int i, j, NBigDec; - int NLINES, done, NREGIONS, nregions; - - f = fopen (GSCFILE, "r"); - if (f == NULL) { - fprintf (stderr, "ERROR: can't find regions file %s\n", GSCFILE); - exit (0); - } - - NREGIONS = 50; - ALLOCATE (regions, GSCRegion, NREGIONS); - nregions = 0; - - if ((minDec <= -90) || (maxDec >= 90)) { - minRa = 0; - maxRa = 360; - } - - /* use the pole regions, if near pole */ - if (maxDec > 86.25) { - sprintf (regions[nregions].filename, "n8230/pole.cpt"); - regions[nregions].RA[0] = 0; - regions[nregions].RA[1] = 360; - regions[nregions].DEC[0] = 86.25; - regions[nregions].DEC[1] = 90.0; - nregions ++; - if (nregions == NREGIONS) { - NREGIONS += 50; - REALLOCATE (regions, GSCRegion, NREGIONS); - } - } - - if (minDec > 86.25) { - *Nregions = nregions; - return (regions); - } - - if ((minDec < 0) && (maxDec > 0)) { - /* Search Both Sides */ - NBigDec = 0; - } else { - /* find large DEC region (directory) */ - NBigDec = -1; - for (i = 0; i < 12; i++) { - if ((minDec >= BigDecBounds[i]) && (minDec < BigDecBounds[i+1])) { - NBigDec = i; - break; - } - } - if (NBigDec < 0) { - for (i = 13; i < 24; i++) { - if ((maxDec < BigDecBounds[i]) && (maxDec >= BigDecBounds[i+1])) { - NBigDec = i; - break; - } - } - } - } - if (NBigDec < 0) { - fprintf (stderr, "ERROR: Dec out of range: %f\n", minDec); - *Nregions = 0; - return ((GSCRegion *) NULL); - } - - /* count lines before section */ - NLINES = 0; - for (i = 0; i < NBigDec; i++) { - NLINES += NDecLines[i]; - } - fseek (f, 5*2880 + 48*NLINES, SEEK_SET); - - /* should be in this section. if not, there is a problem counting... */ - /* careful with the 0,360.0 boundary **/ - done = FALSE; - for (j = 0; !done && (NBigDec + j < 25); j++) { - fread (buffer, 48*NDecLines[NBigDec + j], 1, f); - for (i = 0; (i < NDecLines[NBigDec + j]); i++) { - strncpy (temp, &buffer[i*48], 48); - temp[49] = 0; - hstgsc_hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); - if (RA1 < RA0) RA1 += 360.0; - if ((DEC1 > 0) && (minDec < DEC1) && (maxDec > DEC0) && (minRa < RA1) && (maxRa > RA0)) { - temp[5] = 0; - sprintf (regions[nregions].filename, "%s/%s.cpt", Dec2Sections[NBigDec + j], &temp[1]); - regions[nregions].RA[0] = RA0; - regions[nregions].RA[1] = RA1; - regions[nregions].DEC[0] = DEC0; - regions[nregions].DEC[1] = DEC1; - nregions ++; - if (nregions == NREGIONS) { - NREGIONS += 50; - REALLOCATE (regions, GSCRegion, NREGIONS); - } - } - if ((DEC1 < 0) && (minDec < DEC0) && (maxDec > DEC1) && (minRa < RA1) && (maxRa > RA0)) { - temp[5] = 0; - sprintf (regions[nregions].filename, "%s/%s.cpt", Dec2Sections[NBigDec + j], &temp[1]); - regions[nregions].RA[0] = RA0; - regions[nregions].RA[1] = RA1; - regions[nregions].DEC[0] = DEC0; - regions[nregions].DEC[1] = DEC1; - nregions ++; - if (nregions == NREGIONS) { - NREGIONS += 50; - REALLOCATE (regions, GSCRegion, NREGIONS); - } - } - if (((DEC1 > 0) && (maxDec <= DEC1)) || ((DEC1 < 0) && (minDec >= DEC1))) { - done = TRUE; - } - } - if (done && (minDec < 0) && (maxDec > 0) && (BigDecBounds[NBigDec + j + 1] > 0)) { - /* skip remaining north sections, try south sections */ - /* count lines before section */ - NLINES = 0; - for (i = 0; i < 13; i++) { - NLINES += NDecLines[i]; - } - fseek (f, 5*2880 + 48*NLINES, SEEK_SET); - done = FALSE; - j = 12; - } - } - - REALLOCATE (regions, GSCRegion, MAX (1,nregions)); - *Nregions = nregions; - - fclose (f); - return (regions); - -} Index: trunk/Ohana/src/relastro/src/image-db.c =================================================================== --- trunk/Ohana/src/relastro/src/image-db.c (revision 7080) +++ (revision ) @@ -1,2 +1,0 @@ -# include "relphot.h" - Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 7080) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 8390) @@ -3,14 +3,16 @@ void initialize (int argc, char **argv) { - /* are these set correctly? */ + int N; + ConfigInit (&argc, argv); args (argc, argv); - if ((photcode = GetPhotcodebyName (argv[2])) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[2]); + N = UserPatchSelect ? 1 : 2; + if ((photcode = GetPhotcodebyName (argv[N])) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); exit (1); } if ((photcode[0].type != PHOT_PRI) && (photcode[0].type != PHOT_SEC)) { - fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[2]); + fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[N]); exit (1); } Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 7080) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 8390) @@ -1,39 +1,37 @@ # include "relphot.h" -Catalog *load_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion) { +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) { int i, Nstar; - char filename[256]; Catalog *catalog, tcatalog; if (VERBOSE) fprintf (stderr, "loading catalog data\n"); - ALLOCATE (catalog, Catalog, Nregion); + ALLOCATE (catalog, Catalog, skylist[0].Nregions); - /* load data from each region file, only use bright stars */ - for (i = 0; i < Nregion; i++) { - sprintf (filename, "%s/%s", CATDIR, region[i].filename); - tcatalog.filename = filename; /* don't free region before catalog! */ - switch (lock_catalog (&tcatalog, LCK_SOFT)) { - case 0: - fprintf (stderr, "ERROR: can't lock file %s\n", region[i].filename); - exit (1); - case 1: - gcatalog (&tcatalog, FALSE); /* load from disk */ - break; - case 2: - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); - tcatalog.Naverage = 0; - tcatalog.Nmeasure = 0; - break; - default: - fprintf (stderr, "weird lock_catalog exit state\n"); + // load data from each region file, only use bright stars + for (i = 0; i < skylist[0].Nregions; i++) { + + // set up the basic catalog info + tcatalog.filename = skylist[0].filename[i]; + tcatalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + tcatalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + tcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + + if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename); exit (1); } - bcatalog (&catalog[i], &tcatalog, fullregion); + if (VERBOSE && !tcatalog.Nave_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); + + // select only the brighter stars + bcatalog (&catalog[i], &tcatalog); + dvo_catalog_unlock (&tcatalog); + dvo_catalog_free (&tcatalog); } + // XXX keep this test? Nstar = 0; - for (i = 0; i < Nregion; i++) { + for (i = 0; i < skylist[0].Nregions; i++) { Nstar += catalog[i].Naverage; } @@ -42,4 +40,13 @@ exit (0); } + + // XXX consider only returning the populated catalogs + *Ncatalog = skylist[0].Nregions; return (catalog); } + +/* this function loads all relevant catalog files for the first pass. it currently loads the data + read only (SOFT lock) since it assumes the image table has been locked. if we go to the new + addstar locking paradigm, in which the images and catalogs are updated independently, then we may + need to use an XCLD lock here. +*/ Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 7080) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 8390) @@ -1,35 +1,37 @@ # include "relphot.h" -GSCRegion *load_images (FITS_DB *db, char *seed, int *nregion, GSCRegion *fullregion) { +SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) { - GSCRegion *region; - int Nregion; - Image *image; - int *Line; - int Nimage; + Image *image, *subset; + int Nimage, Nsubset, Nchar; + int *LineNumber; - /* read entire db table */ - if (!dvo_image_load (db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db[0].filename); + SkyTable *sky = NULL; + SkyList *skylist = NULL; - region = name_region (seed, &Nregion); + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area + if (RegionSelect) { + skylist = SkyListByPatch (sky, -1, region); + } else { + Nchar = strlen(regionName); + if (!strcmp (®ionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0; + skylist = SkyListByName (sky, regionName); + } - /* find images which supplied all measurements */ - image = find_images (db, region, Nregion, &Nimage, &Line); - free (Line); + // convert database table to internal structure + image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); - /* find all region files which touch images */ - region = find_regions (image, Nimage, &Nregion, fullregion); - free (image); + // select the images which overlap the selected sky regions + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); - /* find images which supplied all measurements for all regions */ - image = find_images (db, region, Nregion, &Nimage, &Line); - getfullregion (image, Nimage, fullregion); + gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); - gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, Line, Nimage); - - initImages (image, Nimage); - initMosaics (image, Nimage); + initImages (subset, Nsubset); + initMosaics (subset, Nsubset); - *nregion = Nregion; - return (region); + return (skylist); } Index: trunk/Ohana/src/relastro/src/misc.c =================================================================== --- trunk/Ohana/src/relastro/src/misc.c (revision 7080) +++ (revision ) @@ -1,182 +1,0 @@ -# include "relphot.h" - -void sortA (double *X, int N) { - - int l,j,ir,i; - double tX; - - if (N < 2) return; - - l = N >> 1; - ir = N - 1; - for (;;) { - if (l > 0) { - l--; - tX = X[l]; - } - else { - tX = X[ir]; - X[ir] = X[0]; - if (--ir == 0) { - X[0] = tX; - return; - } - } - i = l; - j = (l << 1) + 1; - while (j <= ir) { - if (j < ir && X[j] < X[j+1]) j++; - if (tX < X[j]) { - X[i] = X[j]; - j += (i=j) + 1; - } - else j = ir + 1; - } - X[i] = tX; - } -} - -void sortB (double *X, double *Y, int N) { - - int l,j,ir,i; - double tX, tY; - - if (N < 2) return; - - l = N >> 1; - ir = N - 1; - for (;;) { - if (l > 0) { - l--; - tX = X[l]; - tY = Y[l]; - } - else { - tX = X[ir]; - X[ir] = X[0]; - tY = Y[ir]; - Y[ir] = Y[0]; - if (--ir == 0) { - X[0] = tX; - Y[0] = tY; - return; - } - } - i = l; - j = (l << 1) + 1; - while (j <= ir) { - if (j < ir && X[j] < X[j+1]) j++; - if (tX < X[j]) { - X[i] = X[j]; - Y[i] = Y[j]; - j += (i=j) + 1; - } - else j = ir + 1; - } - X[i] = tX; - Y[i] = tY; - } -} - -void sortC (double *X, double *Y, double *F1, double *F2, int N) { - - int l,j,ir,i; - double tX, tY; - double t1, t2; - - if (N < 2) return; - - l = N >> 1; - ir = N - 1; - for (;;) { - if (l > 0) { - l--; - tX = X[l]; - tY = Y[l]; - t1 = F1[l]; - t2 = F2[l]; - } - else { - tX = X[ir]; - X[ir] = X[0]; - tY = Y[ir]; - Y[ir] = Y[0]; - t1 = F1[ir]; - F1[ir] = F1[0]; - t2 = F2[ir]; - F2[ir] = F2[0]; - if (--ir == 0) { - X[0] = tX; - Y[0] = tY; - F1[0] = t1; - F2[0] = t2; - return; - } - } - i = l; - j = (l << 1) + 1; - while (j <= ir) { - if (j < ir && X[j] < X[j+1]) j++; - if (tX < X[j]) { - X[i] = X[j]; - Y[i] = Y[j]; - F1[i] = F1[j]; - F2[i] = F2[j]; - j += (i=j) + 1; - } - else j = ir + 1; - } - X[i] = tX; - Y[i] = tY; - F1[i] = t1; - F2[i] = t2; - } -} - -void sortD (double *X, double *Y, double *Z, int N) { - - int l,j,ir,i; - double tX, tY, tZ; - - if (N < 2) return; - - l = N >> 1; - ir = N - 1; - for (;;) { - if (l > 0) { - l--; - tX = X[l]; - tY = Y[l]; - tZ = Z[l]; - } - else { - tX = X[ir]; - X[ir] = X[0]; - tY = Y[ir]; - Y[ir] = Y[0]; - tZ = Z[ir]; - Z[ir] = Z[0]; - if (--ir == 0) { - X[0] = tX; - Y[0] = tY; - Z[0] = tZ; - return; - } - } - i = l; - j = (l << 1) + 1; - while (j <= ir) { - if (j < ir && X[j] < X[j+1]) j++; - if (tX < X[j]) { - X[i] = X[j]; - Y[i] = Y[j]; - Z[i] = Z[j]; - j += (i=j) + 1; - } - else j = ir + 1; - } - X[i] = tX; - Y[i] = tY; - Z[i] = tZ; - } -} Index: trunk/Ohana/src/relastro/src/name_region.c =================================================================== --- trunk/Ohana/src/relastro/src/name_region.c (revision 7080) +++ (revision ) @@ -1,24 +1,0 @@ -# include "relphot.h" - -GSCRegion *name_region (char *name, int *Nregions) { - - GSCRegion *region; - Header header; - - ALLOCATE (region, GSCRegion, 1); - - sprintf (region[0].filename, "%s/%s", CATDIR, name); - if (!gfits_read_header (region[0].filename, &header)) { - if (VERBOSE) fprintf (stderr, "catalog file not yet available: %s\n", region[0].filename); - exit (0); - } - gfits_scan (&header, "RA0", "%lf", 1, ®ion[0].RA[0]); - gfits_scan (&header, "RA1", "%lf", 1, ®ion[0].RA[1]); - gfits_scan (&header, "DEC0", "%lf", 1, ®ion[0].DEC[0]); - gfits_scan (&header, "DEC1", "%lf", 1, ®ion[0].DEC[1]); - gfits_free_header (&header); - - *Nregions = 1; - return (region); -} - Index: trunk/Ohana/src/relastro/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/reload_catalogs.c (revision 7080) +++ trunk/Ohana/src/relastro/src/reload_catalogs.c (revision 8390) @@ -1,8 +1,9 @@ # include "relphot.h" -void reload_catalogs (GSCRegion *region, int Nregion, GSCRegion *fullregion) { +void reload_catalogs (SkyList *skylist) { int i; - char filename[256]; + int status; + struct stat filestat; Catalog catalog; @@ -10,23 +11,29 @@ /* load data from each region file */ - for (i = 0; i < Nregion; i++) { - sprintf (filename, "%s/%s", CATDIR, region[i].filename); - catalog.filename = filename; /* don't free region before catalog! */ - switch (lock_catalog (&catalog, LCK_XCLD)) { - case 0: - fprintf (stderr, "ERROR: can't lock file %s\n", region[i].filename); - exit (1); - case 1: - gcatalog (&catalog, TRUE); /* load from disk */ - break; - case 2: - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); - unlock_catalog (&catalog); + for (i = 0; i < skylist[0].Nregions; i++) { + catalog.filename = skylist[0].filename[i]; + + // only update existing db tables + status = stat (catalog.filename, &filestat); + if ((status == -1) && (errno == ENOENT)) { + if (VERBOSE) fprintf (stderr, "no file %s, skipping\n", catalog.filename); continue; - break; - default: - fprintf (stderr, "weird lock_catalog exit state\n"); + } + + 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_MEAS | LOAD_MISS | LOAD_SECF; + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); exit (1); } + if (VERBOSE && (catalog.Nave_disk == 0)) { + fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + initImageBins (&catalog, 1); initMosaicBins (&catalog, 1); @@ -37,6 +44,7 @@ setMrelFinal (&catalog); - wcatalog (&catalog); - unlock_catalog (&catalog); + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); freeImageBins (1); Index: trunk/Ohana/src/relastro/src/relphot.c =================================================================== --- trunk/Ohana/src/relastro/src/relphot.c (revision 7080) +++ trunk/Ohana/src/relastro/src/relphot.c (revision 8390) @@ -1,4 +1,3 @@ # include "relphot.h" -/** I'm not currently setting the lockfiles. this should be fixed! **/ int main (int argc, char **argv) { @@ -6,22 +5,31 @@ int i, status, Ncatalog; Catalog *catalog; - GSCRegion fullregion, *region; FITS_DB db; + + SkyList *skylist = NULL; /* get configuration info, args */ initialize (argc, argv); + /* register database handle with shutdown procedure */ set_db (&db); + + /* lock and load the image db table */ status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); - /* load regions and images based on seed region */ - // XXX EAM : load images and regions based on bounds, not name - region = load_images (&db, argv[1], &Ncatalog, &fullregion); + /* load regions and images based on specified sky patch */ + // XXX need to mimic old-style load by passing patch name + // XXX need to reduce number of global variables in use. + // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true + skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect); + + /* unlock, if we can (else, unlocked below) */ if (!UPDATE) dvo_image_unlock (&db); /* load catalog data from region files */ - catalog = load_catalogs (region, Ncatalog, &fullregion); + catalog = load_catalogs (skylist, &Ncatalog); /* match measurements with images, mosaics */ @@ -80,5 +88,5 @@ /* at this point, we have correct cal coeffs in the image/mosaic structures */ - free_catalogs (catalog, Ncatalog); + for (i = 0; i < Ncatalog; i++) dvo_catalog_free (&catalog[i]); freeImageBins (Ncatalog); freeMosaicBins (Ncatalog); @@ -86,5 +94,5 @@ /* load catalog data from region files, update Mrel include all data */ - reload_catalogs (region, Ncatalog, &fullregion); + reload_catalogs (skylist); setMcalFinal (); dvo_image_update (&db, VERBOSE); Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 8390) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 8390) @@ -0,0 +1,236 @@ +# include "relphot.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). All images in the image catalog are tested once, so there is no check that an + image already has been included. LineNum stores the locations in the Image database of the + list of images */ + +typedef struct { + double Xc[5]; + double Yc[5]; + double Rc; + double Dc; +} SkyRegionCoords; + +Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) { + + Image *image; + int i, j, k, m, found, nimage, NIMAGE; + int InRange, ecode; + double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + int *line_number; + Coords tcoords; + SkyRegionCoords *skycoords; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + // the comparison is made in the catalog local projection. below we set crval1,2 + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + skycoords[i].Yc[4] = skycoords[i].Yc[0]; + + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); + dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + } + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + nimage = 0; + NIMAGE = 100; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, int, NIMAGE); + + // go through the complete list of images, selecting ones which overlap any region + for (i = 0; i < Ntimage; i++) { + + /* exclude images by photcode */ + ecode = GetPhotcodeEquivCodebyCode (timage[i].source); + if (ecode != photcode[0].code) continue; + + /* exclude images by time */ + if (TimeSelect) { + if (timage[i].tzero < TSTART) continue; + if (timage[i].tzero > TSTOP) continue; + } + + /* define image corners */ + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; + found = FALSE; + + /* transform to ra,dec */ + for (j = 0; j < 5; j++) { + XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + } + + /* compare with each region file */ + for (m = 0; (m < skylist[0].Nregions) && !found; m++) { + + /* we make positional comparisons in the projection of catalog */ + tcoords.crval1 = skycoords[m].Rc; + tcoords.crval2 = skycoords[m].Dc; + + /* transform to ra,dec */ + InRange = TRUE; + for (j = 0; (j < 5) && InRange; j++) { + InRange = RD_to_XY (&Xi[j], &Yi[j], Ri[j], Di[j], &tcoords); + } + if (!InRange) continue; + + /* check if image corner inside catalog */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]); + } + /* check if catalog corner inside image */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]); + } + /* check if edges cross */ + for (j = 0; (j < 4) && !found; j++) { + for (k = 0; (k < 4) && !found; k++) { + found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]); + } + } + if (!found) continue; + + image[nimage] = timage[i]; + /* always allow 'few' images to succeed, if possible */ + if (image[nimage].code & ID_IMAGE_FEW) { + image[nimage].code &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR); + } + if (RESET) { + assignMcal (&image[nimage], (double *) NULL, -1); + image[nimage].dMcal_PS = NO_MAG; + image[nimage].code &= ~ID_IMAGE_POOR; + } + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, int, NIMAGE); + } + } + } + + if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, int, MAX (nimage, 1)); + free (skycoords); + + *Nimage = nimage; + *LineNumber = line_number; + return (image); +} + +/* check if line between points 0 and 1 of x1 + crosses line between points 0 and 1 of x2 */ +int edge_check (double *x1, double *y1, double *x2, double *y2) { + + double theta1, theta2; + double Theta1, Theta2; + + theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); + theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + Theta1 = theta1; + Theta2 = theta2; + theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); + theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); + + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + return (TRUE); + +} + +/* check if point x1,y1 is in box formed by x2[0-4] */ +int corner_check (double *x1, double *y1, double *x2, double *y2) { + + int i; + double theta; + + theta = 0; + + for (i = 0; i < 4; i++) { + theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); + } + if (fabs(theta) > 6) { + return (TRUE); + } else { + return (FALSE); + } +} + +/* returns the opening angle between the three points (2 is in middle) + in range -pi to pi */ + +double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) { + + double dx1, dy1, dx2, dy2, ct, st, theta; + + dx1 = x1 - x2; + dy1 = y1 - y2; + + dx2 = x3 - x2; + dy2 = y3 - y2; + + ct = (dx1*dx2 + dy1*dy2); + st = (dx1*dy2 - dx2*dy1); + + theta = atan2 (st, ct); + + return (theta); + +} + Index: trunk/Ohana/src/relastro/src/sort.c =================================================================== --- trunk/Ohana/src/relastro/src/sort.c (revision 8390) +++ trunk/Ohana/src/relastro/src/sort.c (revision 8390) @@ -0,0 +1,182 @@ +# include "relphot.h" + +void sortA (double *X, int N) { + + int l,j,ir,i; + double tX; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + if (--ir == 0) { + X[0] = tX; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + } +} + +void sortB (double *X, double *Y, int N) { + + int l,j,ir,i; + double tX, tY; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + } +} + +void sortC (double *X, double *Y, double *F1, double *F2, int N) { + + int l,j,ir,i; + double tX, tY; + double t1, t2; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + t1 = F1[l]; + t2 = F2[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + t1 = F1[ir]; + F1[ir] = F1[0]; + t2 = F2[ir]; + F2[ir] = F2[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + F1[0] = t1; + F2[0] = t2; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + F1[i] = F1[j]; + F2[i] = F2[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + F1[i] = t1; + F2[i] = t2; + } +} + +void sortD (double *X, double *Y, double *Z, int N) { + + int l,j,ir,i; + double tX, tY, tZ; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tZ = Z[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tZ = Z[ir]; + Z[ir] = Z[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + Z[0] = tZ; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + Z[i] = Z[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + Z[i] = tZ; + } +} Index: trunk/Ohana/src/relastro/src/wcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/wcatalog.c (revision 7080) +++ (revision ) @@ -1,29 +1,0 @@ -# include "relphot.h" - -void wcatalog (Catalog *catalog) { - - /* CATMODE and CATFORMAT determined from catalog */ - - if (!save_catalog (catalog, VERBOSE)) { - fprintf (stderr, "ERROR: failure saving catalog\n"); - exit (1); - } - - if (catalog[0].Naverage) { - free (catalog[0].average); - catalog[0].Naverage = 0; - } - if (catalog[0].Nmeasure) { - free (catalog[0].measure); - catalog[0].Nmeasure = 0; - } - if (catalog[0].Nmissing) { - free (catalog[0].missing); - catalog[0].Nmissing = 0; - } - if (catalog[0].Nsecfilt) { - free (catalog[0].secfilt); - catalog[0].Nsecfilt = 0; - } - gfits_free_header (&catalog[0].header); -} Index: trunk/Ohana/src/relastro/src/wimages.c =================================================================== --- trunk/Ohana/src/relastro/src/wimages.c (revision 7080) +++ (revision ) @@ -1,45 +1,0 @@ -# include "relphot.h" - -/* add 1 image to image db */ -void wimages (FITS_DB *db, Image *image, int Nimage) { - - int status, Nimages; - - /* set permissions somewhere ... - mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; - fchmod (ImageCat, mode); */ - - /* adjust header */ - Nimages = 0; - gfits_scan (&db[0].header, "NIMAGES", "%d", 1, &Nimages); - Nimages ++; - gfits_modify (&db[0].header, "NIMAGES", "%d", 1, Nimages); - - gfits_table_to_vtable (&db[0].ftable, &db[0].vtable, 0, 0); - for (i = 0; i < Nimage; i++) { - gfits_vadd_rows (&db[0].vtable, (char *) image, 1, sizeof(Image)); - } - - /* check that primary header and table header agree */ - if (Nimages != db[0].header.Naxis[1]) { - fprintf (stderr, "header / table length mismatch!\n"); - } - - SetProtect (TRUE); - status = dvo_image_update (db, VERBOSE); - SetProtect (FALSE); - - if (!status) { - fprintf (stderr, "ERROR: failed writing data to image catalog\n"); - exit (1); - } -} - -/* add the new image to the image table: - - use the table header (from ftable) - - create an empty vtable - - add the new image to the end of the vtable - - - write out the vtable, either as TEXT or FITS - * dvo_image_update performs the needed byte-swap on the image data -*/