Index: trunk/Ohana/src/addstar/Makefile =================================================================== --- trunk/Ohana/src/addstar/Makefile (revision 34259) +++ trunk/Ohana/src/addstar/Makefile (revision 34260) @@ -17,19 +17,21 @@ FULL_LDFLAGS = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS) -addstar : $(BIN)/addstar.$(ARCH) -addstard : $(BIN)/addstard.$(ARCH) -addstart : $(BIN)/addstart.$(ARCH) -addstarc : $(BIN)/addstarc.$(ARCH) -mkacc-2mass : $(BIN)/mkacc-2mass.$(ARCH) -sedstar : $(BIN)/sedstar.$(ARCH) -load2mass : $(BIN)/load2mass.$(ARCH) -loadwise : $(BIN)/loadwise.$(ARCH) -loadsupercos: $(BIN)/loadsupercos.$(ARCH) -gztest : $(BIN)/gztest.$(ARCH) -mkcmf : $(BIN)/mkcmf.$(ARCH) +addstar : $(BIN)/addstar.$(ARCH) +addstard : $(BIN)/addstard.$(ARCH) +addstart : $(BIN)/addstart.$(ARCH) +addstarc : $(BIN)/addstarc.$(ARCH) +mkacc-2mass : $(BIN)/mkacc-2mass.$(ARCH) +sedstar : $(BIN)/sedstar.$(ARCH) +load2mass : $(BIN)/load2mass.$(ARCH) +loadwise : $(BIN)/loadwise.$(ARCH) +dumpskycells : $(BIN)/dumpskycells.$(ARCH) +findskycell : $(BIN)/findskycell.$(ARCH) +loadsupercos : $(BIN)/loadsupercos.$(ARCH) +gztest : $(BIN)/gztest.$(ARCH) +mkcmf : $(BIN)/mkcmf.$(ARCH) all: addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells -INSTALL = addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells +INSTALL = addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells findskycell # I need to fix the client/server version of addstar now that I have dropped Stars @@ -299,4 +301,8 @@ $(SRC)/SetSignals.$(ARCH).o \ +FINDSKYCELL = \ +$(SRC)/findskycell.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o + $(ADDSTARC) : $(INC)/addstar.h $(ADDSTARD) : $(INC)/addstar.h @@ -306,4 +312,5 @@ $(SKYCELLS) : $(INC)/addstar.h $(DUMPSKYCELLS) : $(INC)/addstar.h +$(FINDSKYCELL) : $(INC)/addstar.h $(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h $(LOAD-WISE) : $(INC)/addstar.h $(INC)/WISE.h @@ -322,4 +329,5 @@ $(BIN)/skycells.$(ARCH) : $(SKYCELLS) $(BIN)/dumpskycells.$(ARCH) : $(DUMPSKYCELLS) +$(BIN)/findskycell.$(ARCH) : $(FINDSKYCELL) $(BIN)/mkcmf.$(ARCH) : $(MKCMF) Index: trunk/Ohana/src/addstar/doc/tamas_rings.c =================================================================== --- trunk/Ohana/src/addstar/doc/tamas_rings.c (revision 34260) +++ trunk/Ohana/src/addstar/doc/tamas_rings.c (revision 34260) @@ -0,0 +1,59 @@ + +/* + * Prototype implementation of the optimal rings + * by Tamas Budavari + * 2012-07-23 + * + * NOTES: + * - Does only one hemisphere (re-run or mirror) + * - Will leave a hole on top to be covered later + * - Based on PS1 internal draft and the C# code + * - Centers seem to match previous output + * - No further testing has beed done + */ + +#include +#include +#include + +int main(int argc, char* argv[]) +{ + double d2r = M_PI / 180; + + // parameter 'a' is the cell size in degrees + double adeg = 3.955; + + // half of 'a' in radians and its atan + double halfa = adeg / 2 * d2r; + double halftheta = atan(halfa); + + // loop init + double d = 0; // starting Decl. - could change this... + int ring = 0; // ring ID + + while (d < M_PI / 2 - halftheta) + { + double dm = d - halftheta; // eq.5 + if (d == 0) dm = 0; // initial + + int m = (int)ceil(M_PI * cos(dm) / halftheta); // eq.6 + double ip = 2 * M_PI / m; // eq.7 + double dp = atan(tan(d + halftheta) * cos(ip / 2)); // eq.9 + + // printf("%d %d\n", ring, m); // ring & # of cells in that ring + int i; // dump centers of the cells + for (i=0; i -enum {SQUARES, TRIANGLES, LOCAL, RINGS}; +enum {SQUARES, TRIANGLES, LOCAL, RINGS, TAMAS}; enum {TETRAHEDRON, CUBE, OCTOHEDRON, DODECAHEDRON, ICOSAHEDRON}; @@ -95,4 +95,5 @@ int sky_tessellation_squares PROTO((FITS_DB *db, int level, int Nmax)); int sky_tessellation_rings PROTO((FITS_DB *db, int level, int Nmax)); +int sky_tessellation_tamas PROTO((FITS_DB *db, int level, int Nmax)); int sky_triangle_to_image PROTO((Image *image, SkyTriangle *triangle)); @@ -104,4 +105,5 @@ SkyRectangle *sky_rectangle_ring PROTO((float dec, float dDEC, int *nring, char *format)); +SkyRectangle *sky_rectangle_tamas PROTO((double *Dec, double dm, double halfa, double halftheta, int *nring, char *format)); SkyTriangle *sky_divide_triangles PROTO((SkyTriangle *in, int *ntriangles)); Index: trunk/Ohana/src/addstar/src/BoundaryTreeIO.c =================================================================== --- trunk/Ohana/src/addstar/src/BoundaryTreeIO.c (revision 34260) +++ trunk/Ohana/src/addstar/src/BoundaryTreeIO.c (revision 34260) @@ -0,0 +1,291 @@ +# include "addstar.h" + +# define GET_COLUMN_NEW(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +# define GET_COLUMN_RAW(OUT,NAME,TYPE) \ + OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +BoundaryTree *BoundaryTreeLoad(char *filename) { + + int i, j, nz, nb, Ncol; + off_t Nrow; + char type[16]; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + BoundaryTree *tree = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + goto escape; + } + + ALLOCATE (tree, BoundaryTree, 1); + + gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin); + gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset); + + ftable.header = &theader; + + /*** zone information table ***/ + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN_RAW(tree->Nband, "NBAND", int); + GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); + tree->Nzone = Nrow; + + // allocate the storage arrays + ALLOCATE (tree->ra, double *, tree->Nzone); + ALLOCATE (tree->dec, double *, tree->Nzone); + ALLOCATE (tree->cell, int *, tree->Nzone); + ALLOCATE (tree->name, char **, tree->Nzone); + for (i = 0; i < tree->Nzone; i++) { + ALLOCATE (tree->ra[i], double, tree->Nband[i]); + ALLOCATE (tree->dec[i], double, tree->Nband[i]); + ALLOCATE (tree->cell[i], int, tree->Nband[i]); + ALLOCATE (tree->name[i], char *, tree->Nband[i]); + for (j = 0; j < tree->Nband[i]; j++) { + ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH); + } + } + } + + /*** cell information table ***/ + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN_NEW(R, "RA", double); + GET_COLUMN_NEW(D, "DEC", double); + GET_COLUMN_NEW(zone, "ZONE", int); + GET_COLUMN_NEW(band, "BAND", int); + GET_COLUMN_NEW(index, "INDEX", int); + GET_COLUMN_NEW(name, "NAME", char); // XXX how is this done? + gfits_free_header (&theader); + gfits_free_table (&ftable); + fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow); + + // assign the storage arrays + for (i = 0; i < Nrow; i++) { + nz = zone[i]; + nb = band[i]; + tree->ra[nz][nb] = R[i]; + tree->dec[nz][nb] = D[i]; + tree->cell[nz][nb] = i; // XXX ? + memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); + } + + free (R ); + free (D ); + free (zone ); + free (band ); + free (index ); + free (name ); + } + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return tree; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (tree) free (tree); + + fclose (f); + return NULL; +} + +// we are passed a BoundaryTree structure, write it to a FITS table (3 ext) +int BoundaryTreeSave(char *filename, BoundaryTree *tree) { + + int i, nz, nb; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open boundary tree file for output %s\n", filename); + return FALSE; + } + + // we need some information in the header to define the layout + gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin); + gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset); + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + /*** zone information table ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "ZONE_DATA"); + + gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NBAND", "number of cells in each zone", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA_ORIGIN", "origin of ra cell sequence", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA_OFFSET", "offset per cell of ra cell sequence", "degree/cell", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + int *zone = NULL; ALLOCATE (zone, int, tree->Nzone); + + // assign the storage arrays + for (i = 0; i < tree->Nzone; i++) { + zone[i] = i; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "NBAND", tree->Nband, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN", tree->RA_origin, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET", tree->RA_offset, tree->Nzone); + free (zone); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** cell information table ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "CELL_DATA"); + + char fmt[16]; + snprintf (fmt, 16, "%dA", BOUNDARY_TREE_NAME_LENGTH); + gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + int Ncell = 0; + for (i = 0; i < tree->Nzone; i++) { + Ncell += tree->Nband[i]; + } + + // create intermediate storage arrays + // NOTE: we have to unroll the 2D arrays in tree into 1D arrays + double *R ; ALLOCATE (R, double, Ncell); + double *D ; ALLOCATE (D, double, Ncell); + int *zone ; ALLOCATE (zone, int, Ncell); + int *band ; ALLOCATE (band, int, Ncell); + int *index ; ALLOCATE (index, int, Ncell); + char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TREE_NAME_LENGTH); + + // NOTE: a table column of characters must be fixed width, and is passed as a + // contiguous array of Nchar * Nrow values + + // assign the storage arrays + i = 0; + for (nz = 0; nz < tree->Nzone; nz++) { + for (nb = 0; nb < tree->Nband[nz]; nb++) { + R[i] = tree->ra[nz][nb]; + D[i] = tree->dec[nz][nb]; + zone[i] = nz; + band[i] = nb; + index[i] = i; // or tree->cells[nz][nb] ? + memcpy(&name[i*BOUNDARY_TREE_NAME_LENGTH], tree->name[nz][nb], BOUNDARY_TREE_NAME_LENGTH); + i++; + } + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Ncell); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Ncell); + gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, Ncell); + gfits_set_bintable_column (&theader, &ftable, "BAND", band, Ncell); + gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell); + gfits_set_bintable_column (&theader, &ftable, "NAME", name, Ncell); + + free (R ); + free (D ); + free (zone ); + free (band ); + free (index ); + free (name ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +// the boundary tree... +// given an (ra,dec) pair, find the containing projection cell. + +int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { + + // first, find the containing zone + + // if we know dDEC, we can get the bin instantly: + *zone = (dec - tree->DEC_origin) / tree->DEC_offset; + + if (*zone < 0) return FALSE; + if (*zone >= tree->Nzone) return FALSE; + + // now select the RA bin for that zone + *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone]; + + if (*band < 0) return FALSE; + if (*band >= tree->Nband[*zone]) *band = 0; + + return TRUE; +} + Index: trunk/Ohana/src/addstar/src/FilterStars.c =================================================================== --- trunk/Ohana/src/addstar/src/FilterStars.c (revision 34259) +++ trunk/Ohana/src/addstar/src/FilterStars.c (revision 34260) @@ -84,4 +84,19 @@ stars[N].measure.Map += MTIME - dMs; } + if (!isnan(stars[N].measure.Mkron)) { + stars[N].measure.Mkron += MTIME - dMs; + } + if (!isnan(stars[N].measure.FluxPSF)) { + stars[N].measure.FluxPSF /= image[0].exptime; + } + if (!isnan(stars[N].measure.dFluxPSF)) { + stars[N].measure.dFluxPSF /= image[0].exptime; + } + if (!isnan(stars[N].measure.FluxKron)) { + stars[N].measure.FluxKron /= image[0].exptime; + } + if (!isnan(stars[N].measure.dFluxKron)) { + stars[N].measure.dFluxKron /= image[0].exptime; + } // the external ID is supplied, but do we trust it? @@ -93,5 +108,11 @@ double mjd; mjd = ohana_sec_to_mjd (image[0].tzero); - stars[N].measure.extID = CreatePSPSDetectionID(mjd, image[0].ccdnum, stars[N].measure.detID); + int isStack = ((image[0].photcode >= 11000) && (image[0].photcode <= 11400)); + + if (isStack) { + stars[N].measure.extID = CreatePSPSStackDetectionID(image[0].sourceID, image[0].externID, stars[N].measure.detID); + } else { + stars[N].measure.extID = CreatePSPSDetectionID(mjd, image[0].ccdnum, stars[N].measure.detID); + } } else { stars[N].measure.extID = 0; Index: trunk/Ohana/src/addstar/src/ReadStarsFITS.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 34259) +++ trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 34260) @@ -10,5 +10,5 @@ Header theader; FTable table; - Stars *stars; + Stars *stars; // Stars contains Average and Measure if (in_theader == NULL) { @@ -89,23 +89,51 @@ InitStar (&stars[i]); - stars[i].measure.Xccd = smpdata[i].X; - stars[i].measure.Yccd = smpdata[i].Y; + stars[i].measure.Xccd = smpdata[i].X; + stars[i].measure.Yccd = smpdata[i].Y; + stars[i].measure.dXccd = NAN_S_SHORT; // not provided by SMPDATA: + stars[i].measure.dYccd = NAN_S_SHORT; // not provided by SMPDATA: + + stars[i].measure.posangle = NAN_S_SHORT; // not provided by SMPDATA: + stars[i].measure.pltscale = NAN; // not provided by SMPDATA: if ((smpdata[i].M >= ZeroPt) || isnan(smpdata[i].M)) { - stars[i].measure.M = NAN; - stars[i].measure.Map = NAN; - } else { - stars[i].measure.M = smpdata[i].M; - stars[i].measure.Map = smpdata[i].M; - } - - stars[i].measure.dM = smpdata[i].dM*0.001; - + stars[i].measure.M = NAN; + stars[i].measure.Map = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = smpdata[i].M; + stars[i].measure.Map = smpdata[i].M; + stars[i].measure.FluxPSF = pow(10.0, -0.4*smpdata[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * smpdata[i].dM; + } + stars[i].measure.dM = smpdata[i].dM*0.001; + stars[i].measure.dMcal = NAN; // not provided by SMPDATA: + + stars[i].measure.Mkron = NAN; // not provided by SMPDATA: + stars[i].measure.dMkron = NAN; // not provided by SMPDATA: + stars[i].measure.FluxKron = NAN; // not provided by SMPDATA: + stars[i].measure.dFluxKron = NAN; // not provided by SMPDATA: + + stars[i].measure.Sky = NAN; // not provided by SMPDATA: + stars[i].measure.dSky = NAN; // not provided by SMPDATA: + + stars[i].measure.psfChisq = NAN; // not provided by SMPDATA: + stars[i].measure.psfQual = NAN; // not provided by SMPDATA: + stars[i].measure.psfNdof = NAN_S_INT; // not provided by SMPDATA: + stars[i].measure.psfNpix = NAN_S_INT; // not provided by SMPDATA: + stars[i].measure.crNsigma = NAN; // not provided by SMPDATA: + stars[i].measure.extNsigma = NAN; // not provided by SMPDATA: + + stars[i].measure.FWx = ToShortPixels (smpdata[i].fx); + stars[i].measure.FWy = ToShortPixels (smpdata[i].fy); + stars[i].measure.theta = ToShortDegrees (smpdata[i].df); + + stars[i].measure.Mxx = NAN_S_SHORT; // not provided by SMPDATA: + stars[i].measure.Mxy = NAN_S_SHORT; // not provided by SMPDATA: + stars[i].measure.Myy = NAN_S_SHORT; // not provided by SMPDATA: + // the dophot type information gets pushed into the upper 2 bytes of photFlags - stars[i].measure.photFlags = (smpdata[i].dophot << 16); - - stars[i].measure.FWx = ToShortPixels (smpdata[i].fx); - stars[i].measure.FWy = ToShortPixels (smpdata[i].fy); - stars[i].measure.theta = ToShortDegrees (smpdata[i].df); + stars[i].measure.photFlags = (smpdata[i].dophot << 16); } *nstars = Nstars; @@ -131,28 +159,51 @@ for (i = 0; i < Nstars; i++) { InitStar (&stars[i]); - stars[i].measure.Xccd = ps1data[i].X; - stars[i].measure.Yccd = ps1data[i].Y; - - stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); - stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); + stars[i].measure.Xccd = ps1data[i].X; + stars[i].measure.Yccd = ps1data[i].Y; + stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); + stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); + stars[i].measure.posangle = NAN_S_SHORT; // not provided by PS1_DEV_0: + stars[i].measure.pltscale = NAN; // not provided by PS1_DEV_0: + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; - } - stars[i].measure.Map = NAN; - stars[i].measure.dM = ps1data[i].dM; - stars[i].measure.Sky = ps1data[i].sky; - stars[i].measure.dSky = ps1data[i].dSky; - - stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); - stars[i].measure.FWy = ToShortPixels(ps1data[i].fy); - stars[i].measure.theta = ToShortDegrees(ps1data[i].df); - - stars[i].measure.psfChisq = ps1data[i].psfChisq; - stars[i].measure.psfQual = ps1data[i].psfQual; - - stars[i].measure.detID = ps1data[i].detID; + stars[i].measure.M = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.FluxPSF = pow(10.0, -0.4*ps1data[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * ps1data[i].dM; + } + stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = NAN; // not provided by PS1_DEV_0: + stars[i].measure.Map = NAN; // not provided by PS1_DEV_0: + + stars[i].measure.Mkron = NAN; // not provided by PS1_DEV_0: + stars[i].measure.dMkron = NAN; // not provided by PS1_DEV_0: + stars[i].measure.FluxKron = NAN; // not provided by PS1_DEV_0: + stars[i].measure.dFluxKron = NAN; // not provided by PS1_DEV_0: + + stars[i].measure.Sky = ps1data[i].sky; + stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQual = ps1data[i].psfQual; + stars[i].measure.psfNdof = NAN_S_INT; // not provided by PS1_DEV_0: + stars[i].measure.psfNpix = NAN_S_INT; // not provided by PS1_DEV_0: + stars[i].measure.crNsigma = NAN; // not provided by PS1_DEV_0: + stars[i].measure.extNsigma = NAN; // not provided by PS1_DEV_0: + + stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); + stars[i].measure.FWy = ToShortPixels(ps1data[i].fy); + stars[i].measure.theta = ToShortDegrees(ps1data[i].df); + + stars[i].measure.Mxx = NAN_S_SHORT; // not provided by PS1_DEV_0: + stars[i].measure.Mxy = NAN_S_SHORT; // not provided by PS1_DEV_0: + stars[i].measure.Myy = NAN_S_SHORT; // not provided by PS1_DEV_0: + + stars[i].measure.photFlags = 0; // not provided by PS1_DEV_0: + + stars[i].measure.detID = ps1data[i].detID; } *nstars = Nstars; @@ -182,17 +233,37 @@ stars[i].measure.Xccd = ps1data[i].X; stars[i].measure.Yccd = ps1data[i].Y; - stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); + stars[i].measure.posangle = NAN_S_SHORT; // not provided by PS1_DEV_1: + stars[i].measure.pltscale = NAN; // not provided by PS1_DEV_1: + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; - } - stars[i].measure.Map = NAN; + stars[i].measure.M = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.FluxPSF = pow(10.0, -0.4*ps1data[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * ps1data[i].dM; + } stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = NAN; // not provided by PS1_DEV_1: + stars[i].measure.Map = NAN; // not provided by PS1_DEV_1: + + stars[i].measure.Mkron = NAN; // not provided by PS1_DEV_1: + stars[i].measure.dMkron = NAN; // not provided by PS1_DEV_1: + stars[i].measure.FluxKron = NAN; // not provided by PS1_DEV_1: + stars[i].measure.dFluxKron = NAN; // not provided by PS1_DEV_1: + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQual = ps1data[i].psfQual; + stars[i].measure.psfNdof = NAN_S_INT; // not provided by PS1_DEV_1: + stars[i].measure.psfNpix = NAN_S_INT; // not provided by PS1_DEV_1: + stars[i].measure.crNsigma = ps1data[i].crNsigma; + stars[i].measure.extNsigma = ps1data[i].extNsigma; stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); @@ -200,11 +271,12 @@ stars[i].measure.theta = ToShortDegrees(ps1data[i].df); - stars[i].measure.psfChisq = ps1data[i].psfChisq; - stars[i].measure.psfQual = ps1data[i].psfQual; - stars[i].measure.crNsigma = ps1data[i].crNsigma; - stars[i].measure.extNsigma = ps1data[i].extNsigma; - - stars[i].measure.detID = ps1data[i].detID; + stars[i].measure.Mxx = NAN_S_SHORT; // not provided by PS1_DEV_1: + stars[i].measure.Mxy = NAN_S_SHORT; // not provided by PS1_DEV_1: + stars[i].measure.Myy = NAN_S_SHORT; // not provided by PS1_DEV_1: + stars[i].measure.photFlags = ps1data[i].flags; + + // this is may optionally be replaced by the internal sequence (see FilterStars.c) + stars[i].measure.detID = ps1data[i].detID; } *nstars = Nstars; @@ -229,6 +301,6 @@ if (table[0].header[0].Naxis[0] == 136) { - stars = Convert_PS1_V1_Alt (table, nstars); - return (stars); + stars = Convert_PS1_V1_Alt (table, nstars); + return (stars); } @@ -252,15 +324,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.FluxPSF = pow(10.0, -0.4*ps1data[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * ps1data[i].dM; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = NAN; // not provided by PS1_V1: + stars[i].measure.dMkron = NAN; // not provided by PS1_V1: + stars[i].measure.FluxKron = NAN; // not provided by PS1_V1: + stars[i].measure.dFluxKron = NAN; // not provided by PS1_V1: + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -277,5 +358,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; @@ -328,15 +409,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.FluxPSF = pow(10.0, -0.4*ps1data[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * ps1data[i].dM; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = NAN; // not provided by PS1_V1_Alt: + stars[i].measure.dMkron = NAN; // not provided by PS1_V1_Alt: + stars[i].measure.FluxKron = NAN; // not provided by PS1_V1_Alt: + stars[i].measure.dFluxKron = NAN; // not provided by PS1_V1_Alt: + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -353,5 +443,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; @@ -396,15 +486,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + stars[i].measure.FluxPSF = NAN; + stars[i].measure.dFluxPSF = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.FluxPSF = pow(10.0, -0.4*ps1data[i].M); + stars[i].measure.dFluxPSF = stars[i].measure.FluxPSF * ps1data[i].dM; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = NAN; // not provided by PS1_V2: + stars[i].measure.dMkron = NAN; // not provided by PS1_V2: + stars[i].measure.FluxKron = NAN; // not provided by PS1_V2: + stars[i].measure.dFluxKron = NAN; // not provided by PS1_V2: + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -421,5 +520,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; @@ -464,15 +563,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN; + stars[i].measure.dMkron = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN; + + // these fluxes are converted from counts to counts/sec in FilterStars.c + stars[i].measure.FluxPSF = ps1data[i].Flux; + stars[i].measure.dFluxPSF = ps1data[i].dFlux; + stars[i].measure.FluxKron = ps1data[i].kronFlux; + stars[i].measure.dFluxKron = ps1data[i].kronFluxErr; + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -489,5 +597,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; @@ -496,6 +604,9 @@ // the Average fields and the following Measure fields are set in FilterStars after - // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID, - // averef is set in find_matches, dbFlags is zero on ingest. + // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID. + + // averef is set in find_matches + + // dbFlags is zero on ingest. // the following fields are currently not being set anywhere: t_msec @@ -514,6 +625,6 @@ if (table[0].header[0].Naxis[0] == 196) { - stars = Convert_PS1_SV1_Alt (table, nstars); - return (stars); + stars = Convert_PS1_SV1_Alt (table, nstars); + return (stars); } @@ -537,15 +648,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN; + stars[i].measure.dMkron = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN; + + // these fluxes are converted from counts to counts/sec in FilterStars.c + stars[i].measure.FluxPSF = ps1data[i].Flux; + stars[i].measure.dFluxPSF = ps1data[i].dFlux; + stars[i].measure.FluxKron = ps1data[i].kronFlux; + stars[i].measure.dFluxKron = ps1data[i].kronFluxErr; + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -562,5 +682,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; @@ -607,15 +727,24 @@ if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { - stars[i].measure.M = NAN; - } else { - stars[i].measure.M = ps1data[i].M + ZeroPt; + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; } stars[i].measure.dM = ps1data[i].dM; stars[i].measure.dMcal = ps1data[i].dMcal; stars[i].measure.Map = ps1data[i].Map + ZeroPt; - + + stars[i].measure.Mkron = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN; + stars[i].measure.dMkron = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN; + + // these fluxes are converted from counts to counts/sec in FilterStars.c + stars[i].measure.FluxPSF = ps1data[i].Flux; + stars[i].measure.dFluxPSF = ps1data[i].dFlux; + stars[i].measure.FluxKron = ps1data[i].kronFlux; + stars[i].measure.dFluxKron = ps1data[i].kronFluxErr; + stars[i].measure.Sky = ps1data[i].sky; stars[i].measure.dSky = ps1data[i].dSky; - + stars[i].measure.psfChisq = ps1data[i].psfChisq; stars[i].measure.psfQual = ps1data[i].psfQual; @@ -632,5 +761,5 @@ stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); - + stars[i].measure.photFlags = ps1data[i].flags; Index: trunk/Ohana/src/addstar/src/StarOps.c =================================================================== --- trunk/Ohana/src/addstar/src/StarOps.c (revision 34259) +++ trunk/Ohana/src/addstar/src/StarOps.c (revision 34260) @@ -3,6 +3,7 @@ int InitStar (Stars *star) { - memset (&star[0].average, 0, sizeof(Average)); - memset (&star[0].measure, 0, sizeof(Measure)); + + dvo_measure_init (&star[0].measure); + dvo_average_init (&star[0].average); star[0].found = -1; // found == -1 -> not yet found (use enums?) Index: trunk/Ohana/src/addstar/src/args_skycells.c =================================================================== --- trunk/Ohana/src/addstar/src/args_skycells.c (revision 34259) +++ trunk/Ohana/src/addstar/src/args_skycells.c (revision 34260) @@ -37,4 +37,7 @@ if (!strcasecmp (argv[N], "rings")) { MODE = RINGS; + } + if (!strcasecmp (argv[N], "tamas")) { + MODE = TAMAS; } remove_argument (N, &argc, argv); @@ -179,4 +182,16 @@ } remove_argument (N, &argc, argv); + } + if (MODE == TAMAS) { + CELLSIZE = 3.955; + if ((N = get_argument (argc, argv, "-cellsize"))) { + remove_argument (N, &argc, argv); + CELLSIZE = strtod (argv[N], &ptr); + if ((*ptr != 0) || (CELLSIZE < 0.0)) { + fprintf (stderr, "-cellsize requires a floating-point argument\n"); + help (); + } + remove_argument (N, &argc, argv); + } } Index: trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches.c (revision 34259) +++ trunk/Ohana/src/addstar/src/find_matches.c (revision 34260) @@ -250,40 +250,15 @@ if (!IN_REGION (stars[i].average.R, stars[i].average.D)) continue; + 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].dR = 0; - catalog[0].average[Nave].dD = 0; catalog[0].average[Nave].Nmeasure = NSTAR_GROUP; - catalog[0].average[Nave].Nmissing = 0; - catalog[0].average[Nave].Nextend = 0; - catalog[0].average[Nave].measureOffset = Nmeas; - catalog[0].average[Nave].missingOffset = -1; - catalog[0].average[Nave].extendOffset = -1; - - catalog[0].average[Nave].uR = 0; - catalog[0].average[Nave].uD = 0; - catalog[0].average[Nave].duR = 0; - catalog[0].average[Nave].duD = 0; - catalog[0].average[Nave].P = 0; - catalog[0].average[Nave].dP = 0; - - catalog[0].average[Nave].Xp = 0; - catalog[0].average[Nave].ChiSqAve = 0.0; - catalog[0].average[Nave].ChiSqPM = 0.0; - catalog[0].average[Nave].ChiSqPar = 0.0; - catalog[0].average[Nave].Tmean = 0; - catalog[0].average[Nave].Trange = 0; - catalog[0].average[Nave].Npos = 0; - catalog[0].average[Nave].objID = objID; catalog[0].average[Nave].catID = catID; - catalog[0].average[Nave].flags = 0; + if (PSPS_ID) { - catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, - catalog[0].average[Nave].D); - } else { - catalog[0].average[Nave].extID = 0; + catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D); } @@ -291,15 +266,5 @@ for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Map = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Mstdev = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000; - catalog[0].secfilt[Nave*Nsecfilt+j].flags = 0; + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); } Index: trunk/Ohana/src/addstar/src/find_matches_closest.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 34259) +++ trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 34260) @@ -252,39 +252,15 @@ if (!IN_REGION (stars[i].average.R, stars[i].average.D)) continue; + 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].dR = 0; - catalog[0].average[Nave].dD = 0; catalog[0].average[Nave].Nmeasure = NSTAR_GROUP; - catalog[0].average[Nave].Nmissing = 0; - catalog[0].average[Nave].Nextend = 0; - catalog[0].average[Nave].measureOffset = Nmeas; - catalog[0].average[Nave].missingOffset = -1; - catalog[0].average[Nave].extendOffset = -1; - - catalog[0].average[Nave].uR = 0; - catalog[0].average[Nave].uD = 0; - catalog[0].average[Nave].duR = 0; - catalog[0].average[Nave].duD = 0; - catalog[0].average[Nave].P = 0; - catalog[0].average[Nave].dP = 0; - - catalog[0].average[Nave].Xp = 0; - catalog[0].average[Nave].ChiSqAve = 0.0; - catalog[0].average[Nave].ChiSqPM = 0.0; - catalog[0].average[Nave].ChiSqPar = 0.0; - catalog[0].average[Nave].Tmean = 0; - catalog[0].average[Nave].Trange = 0; - catalog[0].average[Nave].Npos = 0; - catalog[0].average[Nave].objID = objID; catalog[0].average[Nave].catID = catID; - catalog[0].average[Nave].flags = 0; + if (PSPS_ID) { catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D); - } else { - catalog[0].average[Nave].extID = 0; } @@ -292,15 +268,5 @@ for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Map = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Mstdev = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000; - catalog[0].secfilt[Nave*Nsecfilt+j].flags = 0; + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); } Index: trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c (revision 34259) +++ trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c (revision 34260) @@ -255,14 +255,16 @@ if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue; + dvo_average_init (&catalog[0].average[Nave]); catalog[0].average[Nave].R = stars[N][0].average.R; catalog[0].average[Nave].D = stars[N][0].average.D; catalog[0].average[Nave].Nmeasure = NREFSTAR_GROUP; - catalog[0].average[Nave].Nmissing = 0; - catalog[0].average[Nave].Nextend = 0; - catalog[0].average[Nave].measureOffset = Nmeas; - catalog[0].average[Nave].missingOffset = -1; - catalog[0].average[Nave].extendOffset = -1; + catalog[0].average[Nave].objID = objID; + catalog[0].average[Nave].catID = catID; + + if (PSPS_ID) { + catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D); + } if (ACCEPT_MOTION) { @@ -275,32 +277,4 @@ catalog[0].average[Nave].P = stars[N][0].average.P; catalog[0].average[Nave].dP = stars[N][0].average.dP; - } else { - catalog[0].average[Nave].dR = 0; - catalog[0].average[Nave].dD = 0; - catalog[0].average[Nave].uR = 0; - catalog[0].average[Nave].uD = 0; - catalog[0].average[Nave].duR = 0; - catalog[0].average[Nave].duD = 0; - catalog[0].average[Nave].P = 0; - catalog[0].average[Nave].dP = 0; - catalog[0].average[Nave].Xp = 0; - } - - catalog[0].average[Nave].Xp = 0; - catalog[0].average[Nave].ChiSqAve = 0.0; - catalog[0].average[Nave].ChiSqPM = 0.0; - catalog[0].average[Nave].ChiSqPar = 0.0; - catalog[0].average[Nave].Tmean = 0; - catalog[0].average[Nave].Trange = 0; - catalog[0].average[Nave].Npos = 0; - - catalog[0].average[Nave].objID = objID; - catalog[0].average[Nave].catID = catID; - catalog[0].average[Nave].flags = 0; - if (PSPS_ID) { - catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, - catalog[0].average[Nave].D); - } else { - catalog[0].average[Nave].extID = 0; } @@ -308,11 +282,5 @@ for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); } Index: trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 34259) +++ trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 34260) @@ -227,14 +227,16 @@ if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue; + dvo_average_init (&catalog[0].average[Nave]); catalog[0].average[Nave].R = stars[N][0].average.R; catalog[0].average[Nave].D = stars[N][0].average.D; catalog[0].average[Nave].Nmeasure = NREFSTAR_GROUP; - catalog[0].average[Nave].Nmissing = 0; - catalog[0].average[Nave].Nextend = 0; - catalog[0].average[Nave].measureOffset = Nmeas; - catalog[0].average[Nave].missingOffset = -1; - catalog[0].average[Nave].extendOffset = -1; + catalog[0].average[Nave].objID = objID; + catalog[0].average[Nave].catID = catID; + + if (PSPS_ID) { + catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D); + } if (ACCEPT_MOTION) { @@ -247,45 +249,10 @@ catalog[0].average[Nave].P = stars[N][0].average.P; catalog[0].average[Nave].dP = stars[N][0].average.dP; - } else { - catalog[0].average[Nave].dR = 0; - catalog[0].average[Nave].dD = 0; - catalog[0].average[Nave].uR = 0; - catalog[0].average[Nave].uD = 0; - catalog[0].average[Nave].duR = 0; - catalog[0].average[Nave].duD = 0; - catalog[0].average[Nave].P = 0; - catalog[0].average[Nave].dP = 0; - catalog[0].average[Nave].Xp = 0; - } - - catalog[0].average[Nave].Xp = 0; - catalog[0].average[Nave].ChiSqAve = 0.0; - catalog[0].average[Nave].ChiSqPM = 0.0; - catalog[0].average[Nave].ChiSqPar = 0.0; - catalog[0].average[Nave].Tmean = 0; - catalog[0].average[Nave].Trange = 0; - catalog[0].average[Nave].Npos = 0; - - catalog[0].average[Nave].objID = objID; - catalog[0].average[Nave].catID = catID; - catalog[0].average[Nave].flags = 0; - if (PSPS_ID) { - catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, - catalog[0].average[Nave].D); - } else { - catalog[0].average[Nave].extID = 0; - } - + } objID ++; for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); } Index: trunk/Ohana/src/addstar/src/findskycell.c =================================================================== --- trunk/Ohana/src/addstar/src/findskycell.c (revision 34260) +++ trunk/Ohana/src/addstar/src/findskycell.c (revision 34260) @@ -0,0 +1,259 @@ +# include "addstar.h" + +static double RA_offset_RINGS_V3[] = {360.000000, 40.000000, 24.000000, 17.142857, 13.333333, 10.909091, 9.230769, 8.000000, 7.200000, + 6.545455, 6.000000, 5.625000, 5.294118, 5.000000, 4.736842, 4.556962, 4.390244, + 4.285714, 4.186047, 4.090909, 4.044944, 4.044944, 4.000000, 4.000000, 4.044944, + 4.044944, 4.090909, 4.186047, 4.285714, 4.390244, 4.556962, 4.736842, 5.000000, + 5.294118, 5.625000, 6.000000, 6.545455, 7.200000, 8.000000, 9.230769, 10.909091, + 13.333333, 17.142857, 24.000000, 40.000000, 360.000000}; + +// in the general case, projection cell centers are arbitrary +// in a more specific case, DEC[i] = DEC_origin + DEC_offset*zone +// in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone] + +enum {TREE_NONE, TREE_MAKE, TREE_USE}; + +void usage (void) { + fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n"); + fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n"); + fprintf (stderr, " (datafile) should contain a list of RA,DEC pairs\n"); + exit (2); +} + +int mktree (char *treefile, char *catdir); +int apply_tree (char *treefile, char *datafile); + +int main (int argc, char **argv) { + + int N; + char *treefile = NULL; + + // what does this program do? + + // 1) load the image table for a tessellation and generate the boundary tree + + // 2) convert RA,DEC (or list?) to cell ID + + if (get_argument (argc, argv, "-help")) usage (); + if (get_argument (argc, argv, "-h")) usage (); + + /* extra error messages */ + int MODE = TREE_NONE; + if ((N = get_argument (argc, argv, "-mktree"))) { + MODE = TREE_MAKE; + remove_argument (N, &argc, argv); + treefile = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-tree"))) { + MODE = TREE_USE; + remove_argument (N, &argc, argv); + treefile = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + // generate the boundary tree + if (argc < 2) usage(); + if (!MODE) usage(); + + if (MODE == TREE_MAKE) { + mktree (treefile, argv[1]); + exit (0); + } + + apply_tree (treefile, argv[1]); + exit (0); +} + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +int mktree (char *treefile, char *catdir) { + + int i, j, zone, band, status; + FITS_DB db; + Image *image; + off_t Nimage; + double x, y, ra, dec; + + char imagefile[DVO_MAX_PATH]; + snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir); + + status = dvo_image_lock (&db, imagefile, 2.0, LCK_XCLD); + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + /* load or create the image table */ + if (db.dbstate == LCK_EMPTY) Shutdown ("can't read image catalog %s", db.filename); + + if (!dvo_image_load (&db, TRUE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + + // generate an empty BoundaryTree + BoundaryTree tree; + tree.FixedGridDEC = TRUE; + tree.FixedGridRA = TRUE; + + // for the moment, I'm going to hardwire the DEC bands to match RINGS.V3 + tree.DEC_origin = -92.0; + tree.DEC_offset = 4.0; + + tree.Nzone = 46; + + ALLOCATE (tree.Nband, int, tree.Nzone); + ALLOCATE (tree.NBAND, int, tree.Nzone); + + ALLOCATE (tree.RA_origin, double, tree.Nzone); + ALLOCATE (tree.RA_offset, double, tree.Nzone); + + ALLOCATE (tree.ra, double *, tree.Nzone); + ALLOCATE (tree.dec, double *, tree.Nzone); + ALLOCATE (tree.cell, int *, tree.Nzone); + ALLOCATE (tree.name, char **, tree.Nzone); + + // NOTE 1: the RA_origin, RA_offsets for RINGS.V3 are defined so that the first projection cell + // overlaps the RA = 0,360 boundary. This make the split a bit of a hack. We end up + // with Nbands, but the max boundary of the last band only goes to 360.0 - + // 0.5*RA_offset. To get the right band number for the boundary region, if the + // calculation for the band number lands beyond the Nbands (ie, band >= Nbands), then + // we need to loop back to the first band. + + // NOTE 2: when we generate the zone & bands initially, we do not know the number of + // bands in the end. we cannot use the test of band >= Nband unless we set an absurdly + // large default value. Thus the value of 1000000 below. + + // assign the bands for RINGS.V3 + for (zone = 0; zone < tree.Nzone; zone++) { + tree.Nband[zone] = 1000000; + tree.NBAND[zone] = 10; + tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone]; + tree.RA_offset[zone] = RA_offset_RINGS_V3[zone]; + ALLOCATE (tree.ra[zone], double, tree.NBAND[zone]); + ALLOCATE (tree.dec[zone], double, tree.NBAND[zone]); + ALLOCATE (tree.cell[zone], int, tree.NBAND[zone]); + ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]); + for (band = 0; band < tree.NBAND[zone]; band++) { + tree.ra[zone][band] = NAN; + tree.dec[zone][band] = NAN; + tree.cell[zone][band] = -1; + ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH); + } + } + + // find the RA,DEC of the image centers & assign to cells + for (i = 0; i < Nimage; i++) { + x = 0.5*image[i].NX; + y = 0.5*image[i].NY; + XY_to_RD (&ra, &dec, x, y, &image[i].coords); + + if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) { + fprintf (stderr, "mismatch!\n"); + continue; + } + // fprintf (stderr, "%d %f %f %f %f %f %d %d\n", i, x, y, tree.RA_offset[zone], ra, dec, zone, band); + + if (band >= tree.NBAND[zone]) { + int start = tree.NBAND[zone]; + tree.NBAND[zone] = band + 10; + REALLOCATE (tree.ra[zone], double, tree.NBAND[zone]); + REALLOCATE (tree.dec[zone], double, tree.NBAND[zone]); + REALLOCATE (tree.cell[zone], int, tree.NBAND[zone]); + REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]); + for (j = start; j < tree.NBAND[zone]; j++) { + tree.ra[zone][j] = NAN; + tree.dec[zone][j] = NAN; + tree.cell[zone][j] = -1; + ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH); + } + } + tree.ra[zone][band] = ra; + tree.dec[zone][band] = dec; + tree.cell[zone][band] = i; + memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH); + } + + // figure out the max band value for each zone? + for (zone = 0; zone < tree.Nzone; zone++) { + int found_last = FALSE; + int last_band = -1; + for (band = 0; band < tree.NBAND[zone]; band++) { + // all cells should be filled + if (tree.cell[zone][band] < 0) { + if (!found_last) { + found_last = TRUE; + last_band = band; + } + } else { + if (found_last) { + fprintf (stderr, "error: empty cell (%d,%d) after last band (%d)\n", zone, band, last_band); + } + } + } + if (last_band == -1) { + last_band = tree.NBAND[zone]; + } + tree.Nband[zone] = last_band; + // fprintf (stderr, "last_band : %d, Nband: %d\n", last_band, tree.Nband[zone]); + } + + struct timeval start, stop; + gettimeofday (&start, (void *) NULL); + + int Npts = 10000000; + + // test : find skycell for NN random points on the sky + long A = time(NULL); + long B = A + 10000; + srand48(B); + for (i = 0; i < Npts; i++) { + ra = 360.0 * drand48(); + dec = 180.0 * drand48() - 90.0; + if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) { + fprintf (stderr, "failure for %f,%f\n", ra, dec); + } + } + MARKTIME("-- test %d pts: %f sec\n", Npts, dtime); + + BoundaryTreeSave (treefile, &tree); + + return TRUE; +} + +int apply_tree (char *treefile, char *datafile) { + + BoundaryTree *tree = BoundaryTreeLoad (treefile); + if (!tree) { + fprintf (stderr, "error loading boundary tree file %s\n", treefile); + exit (2); + } + + FILE *f = fopen (datafile, "r"); + if (!f) { + fprintf (stderr, "error opening data file %s\n", datafile); + exit (3); + } + + double ra, dec; + int Nvalue = 0; + while ((Nvalue = fscanf (f, "%lf %lf", &ra, &dec)) != EOF) { + + int zone, band; + if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) { + fprintf (stderr, "error finding cell for %f,%f\n", ra, dec); + continue; + } + + fprintf (stdout, "%10.6f %10.6f %3d %3d %s\n", ra, dec, zone, band, tree->name[zone][band]); + } + + exit (0); +} Index: trunk/Ohana/src/addstar/src/mkcmf.c =================================================================== --- trunk/Ohana/src/addstar/src/mkcmf.c (revision 34259) +++ trunk/Ohana/src/addstar/src/mkcmf.c (revision 34260) @@ -14,4 +14,5 @@ void gauss_init (int Nbin); double rnd_gauss (double mean, double sigma); +void writeStars_PS1_V3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars); @@ -276,4 +277,8 @@ if (!strcmp(type, "PS1_V2")) { writeStars_PS1_V2 (&ftable, X, Y, M, Flag, Nstars); + found = TRUE; + } + if (!strcmp(type, "PS1_V3")) { + writeStars_PS1_V3 (&ftable, X, Y, M, Flag, Nstars); found = TRUE; } @@ -572,2 +577,60 @@ } +void writeStars_PS1_V3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { + + int i; + CMF_PS1_V3 *stars; + float flux, fSN; + + // XXX add gaussian-distributed noise based on counts + // this needs to make different output 'stars' entries depending on the desired type + ALLOCATE (stars, CMF_PS1_V3, Nstars); + gauss_init (2048); + for (i = 0; i < Nstars; i++) { + stars[i].detID = i; + + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + + stars[i].X = X[i]; + stars[i].Y = Y[i]; + stars[i].M = M[i]; + stars[i].Map = M[i] - 0.05; + + if (ADDNOISE) { + stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M += fSN*rnd_gauss(0.0, 1.0); + } + + // randomly give poor PSFQF values + if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) { + stars[i].psfQual = 0.25; + } else { + stars[i].psfQual = PSFQUAL; + } + + stars[i].dX = FX * fSN; + stars[i].dY = FY * fSN; + stars[i].dM = fSN; + + stars[i].Mpeak = M[i] + 1.0; + stars[i].sky = SKY; + stars[i].dSky = DSKY; + stars[i].psfChisq = PSFCHI; + stars[i].crNsigma = CRN; + stars[i].extNsigma = EXTN; + stars[i].fx = FX; + stars[i].fy = FY; + stars[i].df = DF; + stars[i].nFrames = 1; + stars[i].flags = Flag[i]; + + stars[i].kronFlux = flux * 1.25; + stars[i].kronFluxErr = fSN * flux * 1.25; + } + + gfits_table_set_CMF_PS1_V3 (ftable, stars, Nstars); + gfits_modify (ftable->header, "EXTTYPE", "%s", 1, "PS1_V3"); +} + Index: trunk/Ohana/src/addstar/src/psps_ids.c =================================================================== --- trunk/Ohana/src/addstar/src/psps_ids.c (revision 34259) +++ trunk/Ohana/src/addstar/src/psps_ids.c (revision 34260) @@ -19,4 +19,19 @@ uint64_t +CreatePSPSStackDetectionID(int sourceID, int imageID, int detID) +{ + // sourceID : ID of database + table that tracked the image (< 0x100 = 256) + // imageID : external ID of the image which provided the detections (< 0x1000.0000 ~ 2.7e8) + // detID : detection sequence in image (< 0x1000.0000 ~ 2.7e8) + + assert (detID < 0x10000000); + assert (imageID < 0x10000000); + assert (sourceID < 0x100); + + uint64_t detectid = ((uint64_t)sourceID << 56) + ((uint64_t)imageID << 28) + (uint64_t)detID; + return detectid; +} + +uint64_t CreatePSPSObjectID(double ra, double dec) { Index: trunk/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 34259) +++ trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 34260) @@ -23,4 +23,7 @@ sky_tessellation_rings (db, level, Nmax); return TRUE; + case TAMAS: + sky_tessellation_tamas (db, level, Nmax); + return TRUE; default: break; @@ -284,4 +287,72 @@ free (ring); free (image); + } + return (TRUE); +} + +// the RINGS tessellation uses the declination zones proposed by Tamas Budavari, +// based on code supplied by Tamas 2012.07.23 +int sky_tessellation_tamas (FITS_DB *db, int level, int Nmax) { + + int j, nDEC, Nimage, Nring, Ntotal, Ndigit; + double dec, dDEC; + SkyRectangle *ring; + Image *image; + char format[16]; + + // The tessellation has one input parameter: the approximate cell size. Starting with + // the cell size, determine the optimal projection cell height (dDEC) that results in an + // integer number of dec zones between -90 and +90 + + // in fact, we place a single image on each pole, so the real range of dec is 180.0 - CELLSIZE: + + nDEC = (180.0 - CELLSIZE) / CELLSIZE; + dDEC = (180.0 - CELLSIZE) / nDEC; + nDEC += 2; + + // how many total projection cells for this realization? divide sky area by cell area: + // this is used to set the number of digits, so it does not need to be very accurate... + Ntotal = 41254.2 / (dDEC*dDEC); + Ndigit = (int)(log10(Ntotal)) + 1 ; + snprintf (format, 16, "skycell.%%0%dd", Ndigit); + + double d2r = M_PI / 180; // is RAD_DEG + + // parameter 'a' is the cell size in degrees + double adeg = 3.955; + + // half of 'a' in radians and its atan + double halfa = adeg / 2 * d2r; + double halftheta = atan(halfa); + + // loop init + dec = 0; // starting Decl. - could change this... + + while (dec < M_PI / 2 - halftheta) { + double dm = dec - halftheta; // eq.5 + if (dec == 0) dm = 0; // initial + + // dec is modified by the call below + ring = sky_rectangle_tamas (&dec, dm, halfa, halftheta, &Nring, format); + if (!ring) continue; + + // subdivide each image (Nx x Ny subcells) + Nimage = NX_SUB*NY_SUB*Nring; + ALLOCATE (image, Image, Nimage); + for (j = 0; j < Nring; j++) { + // convert the SkyRectangles to Images for output + sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB); + // printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2); + } + + /* add the new images and save */ + dvo_image_addrows (db, image, Nimage); + SetProtect (TRUE); + dvo_image_update (db, VERBOSE); + SetProtect (FALSE); + dvo_image_clear_vtable (db); + + free (ring); + free (image); } return (TRUE); @@ -666,4 +737,79 @@ } +// define the parameters of a projection centers for this ring +// dec : ~ center of ring in Dec +// dDEC : approximate height +// nring : number of cells generated for this ring +// format : guide to generate the filenames (c-type string format) +SkyRectangle *sky_rectangle_tamas (double *Dec, double dm, double halfa, double halftheta, int *nring, char *format) { + + static int Nname = 0; + int i, j, NX, NY; + SkyRectangle *ring; + + double d2r = M_PI / 180; // is RAD_DEG + double dec = *Dec; + + int nRA = (int)ceil(M_PI * cos(dm) / halftheta); // eq.6 + double dRA = 2 * M_PI / nRA; // eq.7 + double dp = atan(tan(dec + halftheta) * cos(dRA / 2)); // eq.9 + + if (dec == 0.0) { + ALLOCATE (ring, SkyRectangle, nRA); + } else { + ALLOCATE (ring, SkyRectangle, 2*nRA); + } + + for (i = 0; i < nRA; i++) { + // R.A. can use different phase per ring + double ra = i * dRA; // + phase (watch wraparound) + + int npass = (dec == 0.0) ? 1 : 2; + for (j = 0; j < npass; j++) { + + int N = j*nRA + i; + + memset (&ring[N], 0, sizeof(SkyRectangle)); + memset (&ring[N].coords, 0, sizeof(Coords)); + + ring[N].coords.crval1 = ra / d2r; + ring[N].coords.crval2 = (j == 0) ? dec / d2r : -dec / d2r; + + printf(" \t %d %25.20f %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1); + + ring[N].coords.pc1_1 = +1.0 * X_PARITY; + ring[N].coords.pc1_2 = +0.0; + ring[N].coords.pc2_1 = -0.0; + ring[N].coords.pc2_2 = +1.0; + + // range values are in projected degrees + NX = cos(dec - halftheta) * dRA * 3600.0 / SCALE / d2r; + NY = 2 * halftheta * 3600.0 / SCALE / d2r; + + // crpix1,crpix2 is the projection center + ring[N].coords.crpix1 = 0.5*NX; + ring[N].coords.crpix2 = 0.5*NY; + + ring[N].coords.cdelt1 = SCALE / 3600.0; + ring[N].coords.cdelt2 = SCALE / 3600.0; + + strcpy (ring[N].coords.ctype, "DEC--TAN"); + + ring[N].NX = NX*(1.0 + PADDING); + ring[N].NY = NY*(1.0 + PADDING); + ring[N].photcode = 1; // this needs to be set more sensibly + + snprintf (ring[N].name, DVO_IMAGE_NAME_LEN, format, Nname); + Nname++; + } + } + + // advance to next ring + *Dec = halftheta + dp; + + *nring = (dec == 0.0) ? nRA : 2*nRA; + return ring; +} + // an allocated image set is supplied, we fill in the values int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) { @@ -682,6 +828,10 @@ } - Ndigit = (int)(log10(Nx*Ny)) + 1 ; - snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); + if (Nx * Ny > 1) { + Ndigit = (int)(log10(Nx*Ny)) + 1 ; + snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); + } else { + snprintf (format, 24, "%s", input[0].name); + } // if requested extend, the skycell boundaries so that skycells overlap @@ -696,5 +846,10 @@ memcpy (&output[N].coords, &input[0].coords, sizeof(Coords)); - snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N); + if (Nx + Ny > 1) { + snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N); + } else { + snprintf (output[N].name, DVO_IMAGE_NAME_LEN, "%s", format); + } + output[N].NX = NX + 2 * pad_x; output[N].NY = NY + 2 * pad_y; Index: trunk/Ohana/src/addstar/test/simple.dvo =================================================================== --- trunk/Ohana/src/addstar/test/simple.dvo (revision 34259) +++ trunk/Ohana/src/addstar/test/simple.dvo (revision 34260) @@ -21,4 +21,10 @@ test.fields PS1_V2 PS1_V3 test.fields PS1_V3 PS1_V3 + + test.fields PS1_DEV_0 PS1_V4 + test.fields PS1_DEV_1 PS1_V4 + test.fields PS1_V1 PS1_V4 + test.fields PS1_V2 PS1_V4 + test.fields PS1_V3 PS1_V4 end @@ -83,4 +89,13 @@ sort id1 v1 sort id2 v2 + + # some fields require arithmetic manipulations + if ("$name:0" == "KRON_FLUX") + set v1 = -2.5*log(v1) + end + if ("$name:0" == "KRON_FLUX_ERR") + set v1 = KRON_FLUX_ERR / KRON_FLUX + end + set d = v1 - v2 vstat -q d @@ -88,4 +103,12 @@ #echo tapOK fabs($MEAN) < 0.001 "$name:0 vs $name:2 (MEAN)" #echo tapOK fabs($SIGMA) < 0.001 "$name:0 vs $name:2 (SIGMA)" + + # THETA is stored to only (360/65536) deg accuracy + if ("$name:0" == "PSF_THETA") + echo $MEAN + tapOK {abs($MEAN) < 0.006} "$name:0 vs $name:2 (MEAN)" + tapOK {abs($SIGMA) < 0.001} "$name:0 vs $name:2 (SIGMA)" + continue + end tapOK {abs($MEAN) < 0.001} "$name:0 vs $name:2 (MEAN)" @@ -111,4 +134,7 @@ output stdout end + +# the following lists define fields in the cmf files which can be compared to their equivalents in DVO +# the left column is the cmf field name, the right column is the dvo field name # list of cmf fields to test matched to mextract fields @@ -122,5 +148,5 @@ PSF_INST_MAG : mag:inst PSF_INST_MAG_SIG : mag:err - PEAK_FLUX_AS_MAG : SKIP + PEAK_FLUX_AS_MAG : SKIP # not ingested into DVO SKY : sky SKY_SIG : sky_err @@ -130,5 +156,5 @@ PSF_THETA : THETA PSF_QF : PSF_QF - N_FRAMES : SKIP + N_FRAMES : SKIP # not ingested into DVO end @@ -143,5 +169,5 @@ PSF_INST_MAG : mag:inst PSF_INST_MAG_SIG : mag:err - PEAK_FLUX_AS_MAG : SKIP + PEAK_FLUX_AS_MAG : SKIP # not ingested into DVO SKY : sky SKY_SIG : sky_err @@ -153,5 +179,5 @@ PSF_THETA : THETA PSF_QF : PSF_QF - N_FRAMES : SKIP + N_FRAMES : SKIP # not ingested into DVO FLAGS : phot_flags end @@ -236,15 +262,18 @@ X_PSF_SIG : xccd:err # FAIL Y_PSF_SIG : yccd:err # FAIL - RA_PSF : SKIP # astrometry is not calibrated in the cmf - DEC_PSF : SKIP # astrometry is not calibrated in the cmf POSANGLE : SKIP # astrometry is not calibrated in the cmf PLTSCALE : SKIP # astrometry is not calibrated in the cmf PSF_INST_MAG : mag:inst PSF_INST_MAG_SIG : mag:err - AP_MAG_STANDARD : mag:ap # FAIL - AP_MAG_RADIUS : SKIP # no accessor - PEAK_FLUX_AS_MAG : SKIP # no accessor + PSF_INST_FLUX : SKIP # not ingested into DVO + PSF_INST_FLUX_SIG : SKIP # not ingested into DVO + AP_MAG_STANDARD : mag:aperinst # FAIL + AP_MAG_RAW : SKIP # not ingested into DVO + AP_MAG_RADIUS : SKIP # not ingested into DVO CAL_PSF_MAG : SKIP # photometry is not calibrated in the cmf CAL_PSF_MAG_SIG : SKIP # photometry is not calibrated in the cmf + RA_PSF : SKIP # astrometry is not calibrated in the cmf + DEC_PSF : SKIP # astrometry is not calibrated in the cmf + PEAK_FLUX_AS_MAG : SKIP # not ingested into DVO SKY : sky SKY_SIG : sky_err @@ -256,10 +285,21 @@ PSF_THETA : THETA # FAIL PSF_QF : PSF_QF - PSF_NDOF : SKIP # no accessor - PSF_NPIX : SKIP # no accessor - MOMENTS_XX : SKIP # no accessor - MOMENTS_XY : SKIP # no accessor - MOMENTS_YY : SKIP # no accessor + PSF_QF_PERFECT : SKIP # not ingested into DVO + PSF_NDOF : PSF_NDOF + PSF_NPIX : PSF_NPIX + MOMENTS_XX : MXX + MOMENTS_XY : MXY + MOMENTS_YY : MYY + MOMENTS_M3C : SKIP # not ingested into DVO + MOMENTS_M3S : SKIP # not ingested into DVO + MOMENTS_M4C : SKIP # not ingested into DVO + MOMENTS_M4S : SKIP # not ingested into DVO + MOMENTS_R1 : SKIP # not ingested into DVO + MOMENTS_RH : SKIP # not ingested into DVO + KRON_FLUX : mag:kroninst + KRON_FLUX_ERR : mag:kronerr + KRON_FLUX_INNER : SKIP # not ingested into DVO + KRON_FLUX_OUTER : SKIP # not ingested into DVO FLAGS : phot_flags - N_FRAMES : SKIP # no accessor -end + N_FRAMES : SKIP # not ingested into DVO +end Index: trunk/Ohana/src/dvomerge/Makefile =================================================================== --- trunk/Ohana/src/dvomerge/Makefile (revision 34259) +++ trunk/Ohana/src/dvomerge/Makefile (revision 34260) @@ -20,6 +20,8 @@ dvomerge_client : $(BIN)/dvomerge_client.$(ARCH) dvoconvert : $(BIN)/dvoconvert.$(ARCH) +dvoverify : $(BIN)/dvoverify.$(ARCH) +dvoverify_client : $(BIN)/dvoverify_client.$(ARCH) -all: dvomerge dvomerge_client dvoconvert dvosecfilt +all: dvomerge dvomerge_client dvoconvert dvosecfilt dvoverify dvoverify_client # $(SRC)/dvomergeContinue.$(ARCH).o @@ -49,5 +51,4 @@ $(DVOMERGE) : $(INC)/dvomerge.h - $(BIN)/dvomerge.$(ARCH) : $(DVOMERGE) @@ -72,5 +73,4 @@ $(DVOMERGE_CLIENT) : $(INC)/dvomerge.h - $(BIN)/dvomerge_client.$(ARCH) : $(DVOMERGE_CLIENT) @@ -87,5 +87,4 @@ $(DVOCONVERT) : $(INC)/dvomerge.h - $(BIN)/dvoconvert.$(ARCH) : $(DVOCONVERT) @@ -99,5 +98,4 @@ $(DVOSECFILT) : $(INC)/dvomerge.h - $(BIN)/dvosecfilt.$(ARCH) : $(DVOSECFILT) @@ -116,15 +114,25 @@ $(DVOREPAIR) : $(INC)/dvomerge.h - $(BIN)/dvorepair.$(ARCH) : $(DVOREPAIR) DVOVERIFY = \ -$(SRC)/dvoverify.$(ARCH).o +$(SRC)/dvoverify.$(ARCH).o \ +$(SRC)/dvoverify_args.$(ARCH).o \ +$(SRC)/dvoverify_catalogs.$(ARCH).o \ +$(SRC)/dvoverify_utils.$(ARCH).o -$(DVOVERIFY) : $(INC)/dvomerge.h - +$(DVOVERIFY) : $(INC)/dvoverify.h $(BIN)/dvoverify.$(ARCH) : $(DVOVERIFY) -INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvorepair dvoverify +DVOVERIFY_CLIENT = \ +$(SRC)/dvoverify_client.$(ARCH).o \ +$(SRC)/dvoverify_args.$(ARCH).o \ +$(SRC)/dvoverify_catalogs.$(ARCH).o \ +$(SRC)/dvoverify_utils.$(ARCH).o + +$(DVOVERIFY_CLIENT) : $(INC)/dvoverify.h +$(BIN)/dvoverify_client.$(ARCH) : $(DVOVERIFY_CLIENT) + +INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvorepair dvoverify dvoverify_client # dependancy rules for binary code ######################### Index: trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 34259) +++ trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 34260) @@ -17,6 +17,4 @@ # include -# define DVO_MAX_PATH 1024 - int PARALLEL; int PARALLEL_MANUAL; @@ -25,5 +23,4 @@ int HOST_ID; char *HOSTDIR; - int VERBOSE; Index: trunk/Ohana/src/dvomerge/include/dvoverify.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvoverify.h (revision 34260) +++ trunk/Ohana/src/dvomerge/include/dvoverify.h (revision 34260) @@ -0,0 +1,53 @@ +# include +# include +# include +# include +# include +# include + +/* solaris requires both of these instead of ip.h: + # include + # include +*/ + +/* linux is happy with this, not solaris */ +# include +# include +# include +# include + +# define DEBUG 0 + +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +char *CATDIR; + +int HOST_ID; +char *HOSTDIR; + +char *RESULTS; + +int CHECKSORTED; +int VERBOSE; +int NNotSorted; +int CHECK_TOPLEVEL; +int LIST_MISSING; + +SkyRegion UserPatch; + +int dvoverify_args (int *argc, char **argv); +int dvoverify_client_args (int *argc, char **argv); + +int dvoverify_catalogs (SkyList *skylist, int *Nbad); +int dvoverify_parallel (SkyList *skylist, int *Nbad); + +int dvoverify_single (char *filename); + +int VerifyTableFile (char *filename); +int CheckCatalogIndexes (char *filename, SkyRegion *region); + +void InitFailures (); +void AddFailures (char *filename); +char **GetFailures (int *N); Index: trunk/Ohana/src/dvomerge/src/args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/args.c (revision 34259) +++ trunk/Ohana/src/dvomerge/src/args.c (revision 34260) @@ -50,6 +50,4 @@ } - // XXX for the moment, make this selection manual. it needs to be automatic - // based on the state of the SkyTable PARALLEL = FALSE; if ((N = get_argument (*argc, argv, "-parallel"))) { Index: trunk/Ohana/src/dvomerge/src/dvomergeContinue.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 34259) +++ trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 34260) @@ -1,2 +1,5 @@ +NOTE: +/// this is now not used; it has been merged with dvomergeUpdate + # include "dvomerge.h" Index: trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 34259) +++ trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 34260) @@ -1,2 +1,5 @@ +NOTE: +/// this is now not used; it has been merged with dvomergeUpdate + # include "dvomerge.h" # include Index: trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 34259) +++ trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 34260) @@ -94,9 +94,10 @@ fclose (fout); - // update header of output catalog + // XXX note that we are hardwired to v 2 + // check the header of output catalog for an existing merge long long last_size; char last_moddate[80]; - gfits_scan (&outheader, "LMRG_SZ", "%lld", 1, &last_size); - gfits_scan (&outheader, "LMRG_DT", "%s", 1, last_moddate); + gfits_scan (&outheader, "LMRG_SZ2", "%lld", 1, &last_size); + gfits_scan (&outheader, "LMRG_DT2", "%s", 1, last_moddate); time_t last_mod = ohana_date_to_sec (last_moddate); @@ -150,5 +151,6 @@ outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; - // update header of output catalog +# if (0) + // get the LMRG data from the previous pass long long last_size; char last_moddate[80]; @@ -156,5 +158,5 @@ int status_date = gfits_scan (&outcatalog.header, "LMRG_DT", "%s", 1, last_moddate); - // save the previous size of merged entry + // save the LMRG data from the previous pass if (status_size && status_date) { gfits_modify (&outcatalog.header, "LMRG_SZ1", "%lld", 1, (long long) last_size); @@ -165,4 +167,10 @@ gfits_modify (&outcatalog.header, "LMRG_SZ", "%lld", 1, (long long) instats.st_size); gfits_modify (&outcatalog.header, "LMRG_DT", "%s", 1, moddate); +# endif + + // update header of output catalog + // XXX note that we are hardwired to v 2 + gfits_modify (&outcatalog.header, "LMRG_SZ2", "%lld", 1, (long long) instats.st_size); + gfits_modify (&outcatalog.header, "LMRG_DT2", "%s", 1, moddate); if (!dvo_catalog_backup (&outcatalog, TRUE)) { Index: trunk/Ohana/src/dvomerge/src/dvoverify.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 34259) +++ trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 34260) @@ -1,3 +1,3 @@ -# include "dvomerge.h" +# include "dvoverify.h" /* things we can verify easily: @@ -8,108 +8,21 @@ */ -int VerifyTableFile (char *filename); -int CheckCatalogIndexes (char *catdir, char *filename, SkyRegion *region); - -# define DEBUG 0 - -int VERBOSE = FALSE; -int NNotSorted = 0; - int main (int argc, char **argv) { - char filename[1024]; + SkyTable *sky; + SkyList *skylist; - int N, Nbad; - off_t i; - SkyTable *insky; - SkyList *inlist; - SkyRegion UserPatch; - // Catalog catalog; + // check various options + dvoverify_args (&argc, argv); + CATDIR = argv[1]; - int CHECKSORTED; - - if ((N = get_argument (argc, argv, "-v"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-verbose"))) { - VERBOSE = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-s"))) { - CHECKSORTED = TRUE; - remove_argument (N, &argc, argv); - } - if ((N = get_argument (argc, argv, "-sorted"))) { - CHECKSORTED = TRUE; - remove_argument (N, &argc, argv); - } - - if ((N = get_argument (argc, argv, "-cpt"))) { - remove_argument (N, &argc, argv); - char *filename = strcreate (argv[N]); - remove_argument (N, &argc, argv); - - int isBad = FALSE; - - if (!VerifyTableFile (filename)) { - fprintf (stderr, "bad average table %s\n", filename); - isBad = TRUE; - } - - // change last 't' to 's': - int Nlast; - Nlast = strlen(filename) - 1; - filename[Nlast] = 's'; - if (!VerifyTableFile (filename)) { - fprintf (stderr, "bad secfilt table %s\n", filename); - isBad = TRUE; - } - - // change last 't' to 's': - Nlast = strlen(filename) - 1; - filename[Nlast] = 'm'; - if (!VerifyTableFile (filename)) { - fprintf (stderr, "bad measure table %s\n", filename); - isBad = TRUE; - } - - if (isBad) exit (1); - exit (0); - } - - // restrict to a portion of the sky - UserPatch.Rmin = 0; - UserPatch.Rmax = 360; - 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); - 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); - } - - if (argc != 2) { - fprintf (stderr, "USAGE: dvoverify (catdir) [-region Rmin Rmax Dmin Dmax] [-v] [-s]\n\n"); - fprintf (stderr, " -v : VERBOSE\n"); - fprintf (stderr, " -s : checks if sorted, return error if not\n"); - fprintf (stderr, " (catdir) : database of interest\n"); - exit (2); - } - - char *catdir = argv[1]; - - Nbad = 0; + int Nbad = 0; // XXX make this step optional - if (1) { + if (CHECK_TOPLEVEL) { + char filename[DVO_MAX_PATH]; + // check the photcode table - sprintf (filename, "%s/Photcodes.dat", catdir); + sprintf (filename, "%s/Photcodes.dat", CATDIR); if (!VerifyTableFile (filename)) { Nbad ++; @@ -117,5 +30,5 @@ // check the skytable - char *skyfile = SkyTableFilename (catdir); + char *skyfile = SkyTableFilename (CATDIR); if (!VerifyTableFile (skyfile)) { Nbad ++; @@ -123,5 +36,5 @@ // check the image table - sprintf (filename, "%s/Images.dat", catdir); + sprintf (filename, "%s/Images.dat", CATDIR); if (!VerifyTableFile (filename)) { Nbad ++; @@ -130,38 +43,17 @@ // load the sky table for the existing database - insky = SkyTableLoadOptimal (catdir, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); - myAssert(insky, "can't read SkyTable"); - SkyTableSetFilenames (insky, catdir, "cpt"); - inlist = SkyListByPatch (insky, -1, &UserPatch); + sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + myAssert(sky, "can't read SkyTable"); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &UserPatch); - // loop over all catalogs, save to output catalogs - for (i = 0; i < inlist[0].Nregions; i++) { - if (!inlist[0].regions[i][0].table) continue; - if (i % 1000 == 0) fprintf (stderr, "."); + dvoverify_catalogs (skylist, &Nbad); - // sprintf (filename, "%s/%s.cpt", catdir, inlist[0].regions[i][0].name); - if (!VerifyTableFile (inlist[0].filename[i])) { - Nbad ++; - } + int i, Nfailures; + char **failures = GetFailures (&Nfailures); - sprintf (filename, "%s/%s.cps", catdir, inlist[0].regions[i][0].name); - if (!VerifyTableFile (filename)) { - Nbad ++; - } - - sprintf (filename, "%s/%s.cpm", catdir, inlist[0].regions[i][0].name); - if (!VerifyTableFile (filename)) { - Nbad ++; - } - - if (!CheckCatalogIndexes(catdir, inlist[0].filename[i], inlist[0].regions[i])){ - Nbad ++; - } - - // exit immediately if any file are unsorted and we require sorted tables - if (CHECKSORTED && NNotSorted) { - fprintf (stderr, "ERROR: files are not sorted\n"); - exit (1); - } + fprintf (stderr, "---- files with errors ---- \n"); + for (i = 0; i < Nfailures; i++) { + fprintf (stderr, "%s\n", failures[i]); } @@ -178,237 +70,2 @@ } -// is this file a consistent FITS file? -int VerifyTableFile (char *filename) { - - int status, Next; - off_t Nbytes; - Header header; - - struct stat fileStats; - FILE *file; - - // does the file exist? - status = stat (filename, &fileStats); - if (status) { - // some error accessing the file. there is only one acceptable error: file not found - switch (errno) { - case ENOENT: - if (DEBUG) fprintf (stderr, "file does not exist, skipping %s\n", filename); - return TRUE; - case ENOMEM: - fprintf (stderr, "Out of memory: %s\n", filename); - return TRUE; - case EACCES: - fprintf (stderr, "Permission error on %s\n", filename); - return FALSE; - case EFAULT: - fprintf (stderr, "Bad address: %s\n", filename); - return FALSE; - case ELOOP: - fprintf (stderr, "Too many symbolic links encountered while traversing the path: %s\n", filename); - return FALSE; - case ENAMETOOLONG: - fprintf (stderr, "File name too long: %s\n", filename); - return FALSE; - case ENOTDIR: - fprintf (stderr, "A component of the path is not a directory: %s\n", filename); - return FALSE; - case EOVERFLOW: - fprintf (stderr, "file too large for program version: %s\n", filename); - return FALSE; - default: - fprintf (stderr, "unknown error: %s\n", filename); - return FALSE; - } - } - - // does it have any data? - if (fileStats.st_size == 0) { - fprintf (stderr, "file is empty: %s\n", filename); - return FALSE; - } - - // can we open it? - file = fopen(filename, "r"); - if (!file) { - fprintf (stderr, "unable to open valid file: %s\n", filename); - return FALSE; - } - - // scan all extentions - Nbytes = 0; - Next = -1; - if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); - while (Nbytes < fileStats.st_size) { - - // Check on the PHU - if (!gfits_fread_header (file, &header)) { - if (Next == -1) { - fprintf (stderr, "unable to read PHU header for %s\n", filename); - } else { - fprintf (stderr, "unable to read header for %s, extension %d (or file has excess bytes)\n", filename, Next); - } - fclose(file); - return (FALSE); - } - - // move to TBL header - Nbytes += header.datasize + gfits_data_size (&header); - if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); - if (Nbytes > fileStats.st_size) { - fprintf (stderr, "file is short ("OFF_T_FMT" vs "OFF_T_FMT"): %s\n", Nbytes, fileStats.st_size, filename); - gfits_free_header(&header); - fclose (file); - return FALSE; - } - gfits_free_header(&header); - - status = fseeko (file, Nbytes, SEEK_SET); - if (status) { - switch (errno) { - case EBADF: - fprintf (stderr, "something wrong with file handle: %s\n", filename); - fclose (file); - return FALSE; - case EINVAL: - fprintf (stderr, "invalid offset: %s\n", filename); - fclose (file); - return FALSE; - default: - fprintf (stderr, "other error in fseeko: %s\n", filename); - fclose (file); - return FALSE; - } - } - Next ++; - } - if (DEBUG) fprintf (stderr, "file is good: %s\n", filename); - fclose (file); - return TRUE; -} - -// CheckCatalogIndexes(catdir, inlist[0].regions[i][0].name); - -int CheckCatalogIndexes (char *catdir, char *filename, SkyRegion *region) { - - Catalog catalog; - int i, j, m, status; - - status = TRUE; - - // set the parameters which guide catalog open/load/create - catalog.filename = filename; - 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; - catalog.Nsecfilt = 0; - - // an error exit status here is a significant error (disk I/O or file access) - if (!dvo_catalog_open (&catalog, region, VERBOSE, "r")) { - fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); - return FALSE; - } - - // Naves_disk == 0 implies an empty catalog file, skip empty catalogs - if (catalog.Naves_disk == 0) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - return TRUE; - } - - // if the table is SORTED, then the following can be checked - // check the following: - // measure[j].averef -> average[averef] - // measure[j].objID = average[averef].objID - // measure[j].catID = average[averef].catID - // measure[j].measureOffset < Nmeasure - // \sum average[].Nmeasure = Nmeasure - - // if the table is NOT SORTED, we have a subset of checks we can make - if (!catalog.sorted) { - fprintf (stderr, "!"); - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - NNotSorted++; - if (VERBOSE) fprintf (stderr, "file is not sorted: %s\n", filename); - return TRUE; - } - - int NmeasureTotal = 0; - int measureOffsetOK = TRUE; - for (i = 0; i < catalog.Naverage; i++) { - NmeasureTotal += catalog.average[i].Nmeasure; - if (VERBOSE && !(NmeasureTotal <= catalog.Nmeasure)) { - fprintf (stderr, "NmeasureTotal > catalog.Nmeasure: %d %d %d\n", i, catalog.average[i].Nmeasure, (int) catalog.Nmeasure); - } - measureOffsetOK &= (catalog.average[i].measureOffset < catalog.Nmeasure); - if (VERBOSE && !(catalog.average[i].measureOffset < catalog.Nmeasure)) { - fprintf (stderr, "measureOffset >= catalog.Nmeasure: %d %d %d\n", i, catalog.average[i].measureOffset, (int) catalog.Nmeasure); - } - measureOffsetOK &= (catalog.average[i].measureOffset + catalog.average[i].Nmeasure <= catalog.Nmeasure); - if (VERBOSE && !(catalog.average[i].measureOffset + catalog.average[i].Nmeasure <= catalog.Nmeasure)) { - fprintf (stderr, "measureOffset + Nmeasure > catalog.Nmeasure : %d %d %d\n", i, catalog.average[i].Nmeasure, (int) catalog.Nmeasure); - } - } - - if (!measureOffsetOK) { - fprintf (stderr, "ERROR: catalog %s has an invalid measureOffset\n", catalog.filename); - status = FALSE; - } - - if (NmeasureTotal != catalog.Nmeasure) { - fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog.filename); - status = FALSE; - } - - if (!status) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - return (status); - } - - int objIDsOK = TRUE; - int catIDsOK = TRUE; - int averefOK = TRUE; - - for (i = 0; i < catalog.Naverage; i++) { - m = catalog.average[i].measureOffset; - for (j = 0; j < catalog.average[i].Nmeasure; j++) { - objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID); - catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID); - averefOK &= (catalog.measure[m+j].averef == i); - } - } - - if (!objIDsOK) { - fprintf (stderr, "ERROR: catalog %s has invalid obj IDs\n", catalog.filename); - status = FALSE; - } - if (!catIDsOK) { - fprintf (stderr, "ERROR: catalog %s has invalid cat IDs\n", catalog.filename); - status = FALSE; - } - if (!averefOK) { - fprintf (stderr, "ERROR: catalog %s has invalid averef values\n", catalog.filename); - status = FALSE; - } - -// for (i = 0; i < catalog.Naverage; i++) { -// m = catalog.average[i].measureOffset; -// for (j = 0; i < catalog.Nmeasure; i++) { -// objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID); -// catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID); -// averefOK &= (catalog.measure[m+j].averef = i); -// } -// } - - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - - return status; -} - -// gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow); -// gfits_scan(&cpmHeaderTBL, "NAXIS2", "%d", 1, &Nrows); - - Index: trunk/Ohana/src/dvomerge/src/dvoverify_args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify_args.c (revision 34260) +++ trunk/Ohana/src/dvomerge/src/dvoverify_args.c (revision 34260) @@ -0,0 +1,194 @@ +# include "dvoverify.h" + +int dvoverify_args (int *argc, char **argv) { + + int N; + + NNotSorted = 0; + + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-verbose"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + CHECKSORTED = FALSE; + if ((N = get_argument (*argc, argv, "-s"))) { + CHECKSORTED = TRUE; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-sorted"))) { + CHECKSORTED = TRUE; + remove_argument (N, argc, argv); + } + + LIST_MISSING = FALSE; + if ((N = get_argument (*argc, argv, "-list-missing"))) { + LIST_MISSING = TRUE; + remove_argument (N, argc, argv); + } + + if ((N = get_argument (*argc, argv, "-cpt"))) { + remove_argument (N, argc, argv); + char *filename = strcreate (argv[N]); + remove_argument (N, argc, argv); + + if (*argc != 1) { + fprintf (stderr, "USAGE: dvoverify -cpt filename.cpt\n"); + fprintf (stderr, "OPTIONS: -v -verbose -s -sorted\n"); + fprintf (stderr, "NOTE: -parallel and other options not allowed for -cpt mode\n"); + exit (2); + } + + int isGood = dvoverify_single (filename); + if (!isGood) exit (1); + exit (0); + } + + CHECK_TOPLEVEL = TRUE; + if ((N = get_argument (*argc, argv, "-skip-toplevel"))) { + CHECK_TOPLEVEL = FALSE; + remove_argument (N, argc, argv); + } + + // restrict to a portion of the sky + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + 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); + 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); + } + + HOST_ID = 0; + HOSTDIR = NULL; + + PARALLEL = FALSE; + if ((N = get_argument (*argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // relphot 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"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, argc, argv); + } + // this is a test mode : rather than launching the relphot_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (*argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, argc, argv); + } + + if (*argc != 2) { + fprintf (stderr, "USAGE: dvoverify (catdir) [-region Rmin Rmax Dmin Dmax] [-v] [-s]\n\n"); + fprintf (stderr, " (catdir) : database of interest\n"); + fprintf (stderr, " -v : VERBOSE\n"); + fprintf (stderr, " -s : checks if sorted, return error if not\n"); + fprintf (stderr, " -region : limit checks to specified region\n"); + fprintf (stderr, " -parallel : run in parallel across cluster\n"); + fprintf (stderr, " -skip-toplevel : do not check top-level files (Images, Photcodes, etc)\n\n"); + fprintf (stderr, " OR : -cpt (filename)\n"); + exit (2); + } + + return TRUE; +} + +int dvoverify_client_args (int *argc, char **argv) { + + int N; + + VERBOSE = FALSE; + CHECKSORTED = FALSE; + NNotSorted = 0; + + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-s"))) { + CHECKSORTED = TRUE; + remove_argument (N, argc, argv); + } + + LIST_MISSING = FALSE; + if ((N = get_argument (*argc, argv, "-list-missing"))) { + LIST_MISSING = TRUE; + remove_argument (N, argc, argv); + } + + // restrict to a portion of the sky + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + 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); + 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); + } + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + HOST_ID = 0; + if ((N = get_argument (*argc, argv, "-hostID"))) { + remove_argument (N, argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, argc, argv); + } + + HOSTDIR = NULL; + if ((N = get_argument (*argc, argv, "-hostdir"))) { + remove_argument (N, argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + RESULTS = NULL; + if ((N = get_argument (*argc, argv, "-results"))) { + remove_argument (N, argc, argv); + RESULTS = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + if (!HOST_ID || !HOSTDIR || (*argc != 2)) { + fprintf (stderr, "USAGE: dvoverify_client (catdir) -results (file) -hostID ID -hostdir DIR [-region Rmin Rmax Dmin Dmax] [-v] [-s]\n\n"); + fprintf (stderr, " (catdir) : database of interest\n"); + fprintf (stderr, " -v : VERBOSE\n"); + fprintf (stderr, " -s : checks if sorted, return error if not\n"); + exit (2); + } + + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvoverify_catalogs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify_catalogs.c (revision 34260) +++ trunk/Ohana/src/dvomerge/src/dvoverify_catalogs.c (revision 34260) @@ -0,0 +1,151 @@ +# include "dvoverify.h" + +int dvoverify_catalogs (SkyList *skylist, int *nbad) { + + int i; + char filename[DVO_MAX_PATH]; + + InitFailures (); + + if (PARALLEL && !HOST_ID) { + int status = dvoverify_parallel (skylist, nbad); + return status; + } + + int Nbad = *nbad; + + char *mycatdir = HOST_ID ? HOSTDIR : CATDIR; + + // loop over all catalogs, save to output catalogs + for (i = 0; i < skylist[0].Nregions; i++) { + if (!skylist[0].regions[i][0].table) continue; + if (i % 1000 == 0) fprintf (stderr, "."); + + int skipIndexCheck = FALSE; + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + sprintf (filename, "%s/%s.cpt", mycatdir, skylist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + skipIndexCheck = TRUE; + AddFailures (filename); + } + + sprintf (filename, "%s/%s.cps", mycatdir, skylist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + AddFailures (filename); + } + + sprintf (filename, "%s/%s.cpm", mycatdir, skylist[0].regions[i][0].name); + if (!VerifyTableFile (filename)) { + Nbad ++; + skipIndexCheck = TRUE; + AddFailures (filename); + } + + sprintf (filename, "%s/%s.cpt", mycatdir, skylist[0].regions[i][0].name); + if (!skipIndexCheck) { + if (!CheckCatalogIndexes(filename, skylist[0].regions[i])){ + Nbad ++; + AddFailures (filename); + } + } + + // exit immediately if any file are unsorted and we require sorted tables + // this probably goes elsewhere... + if (CHECKSORTED && NNotSorted) { + fprintf (stderr, "ERROR: files are not sorted\n"); + return FALSE; + } + } + + *nbad = Nbad; + + return TRUE; +} + +// launch the dvoverify_client jobs to the parallel hosts +int dvoverify_parallel (SkyList *skylist, int *Nbad) { + + // ensure that the paths are absolute path names + char *abscatdir = abspath (CATDIR, DVO_MAX_PATH); + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, skylist->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", skylist->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/dvoverify.dat", table->hosts[i].pathname); + + // options / arguments that can affect relastro_client -update-objects: + char command[DVO_MAX_PATH]; + snprintf (command, DVO_MAX_PATH, "dvoverify_client %s -results %s -hostID %d -hostdir %s -region %f %f %f %f", + abscatdir, table->hosts[i].results, table->hosts[i].hostID, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax + ); + + char tmpline[DVO_MAX_PATH]; + if (VERBOSE) { snprintf (tmpline, DVO_MAX_PATH, "%s -v", command); strcpy (command, tmpline); } + if (CHECKSORTED) { snprintf (tmpline, DVO_MAX_PATH, "%s -s", command); strcpy (command, tmpline); } + if (LIST_MISSING) { snprintf (tmpline, DVO_MAX_PATH, "%s -list-missing", command); strcpy (command, tmpline); } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running photdbc_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the photdbc_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + int NbadHost = 0; + int NNotSortedHost = 0; + for (i = 0; i < table->Nhosts; i++) { + FILE *results = fopen (table->hosts[i].results, "r"); + if (!results) { + fprintf (stderr, "cannot read results for %d: %s\n", i, table->hosts[i].results); + continue; + } + fscanf (results, "%*s %d %*s %d", &NbadHost, &NNotSortedHost); + *Nbad += NbadHost; + NNotSorted += NNotSortedHost; + fclose (results); + } + + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvoverify_client.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify_client.c (revision 34260) +++ trunk/Ohana/src/dvomerge/src/dvoverify_client.c (revision 34260) @@ -0,0 +1,46 @@ +# include "dvoverify.h" + +int main (int argc, char **argv) { + + int Nbad; + SkyTable *sky; + SkyList *skylist; + + // check various options + dvoverify_client_args (&argc, argv); + CATDIR = argv[1]; + + Nbad = 0; + + // load the sky table for the existing database + sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + myAssert(sky, "can't read SkyTable"); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &UserPatch); + + dvoverify_catalogs (skylist, &Nbad); + + // write result to file + FILE *results = fopen (RESULTS, "w"); + fprintf (results, "Nbad: %d NNotSorted: %d\n", Nbad, NNotSorted); + + int i, Nfailures; + char **failures = GetFailures (&Nfailures); + fprintf (results, "---- files with errors ---- \n"); + for (i = 0; i < Nfailures; i++) { + fprintf (results, "%s\n", failures[i]); + } + fclose (results); + + if (Nbad > 0) { + fprintf (stderr, "ERROR: %d files are bad\n", Nbad); + exit (1); + } + + fprintf (stderr, "SUCCESS on %d: no files are bad\n", HOST_ID); + if (NNotSorted) { + fprintf (stderr, "NOTE: %d files are not sorted\n", NNotSorted); + } + exit (0); +} + Index: trunk/Ohana/src/dvomerge/src/dvoverify_utils.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify_utils.c (revision 34260) +++ trunk/Ohana/src/dvomerge/src/dvoverify_utils.c (revision 34260) @@ -0,0 +1,283 @@ +# include "dvoverify.h" + +static int Nfailures = 0; +static int NFAILURES = 100; +static char **failures = NULL; + +void InitFailures () { + ALLOCATE (failures, char *, NFAILURES); +} + +void AddFailures (char *filename) { + failures[Nfailures] = strcreate (filename); + Nfailures ++; + CHECK_REALLOCATE (failures, char *, NFAILURES, Nfailures, 100); +} + +char **GetFailures (int *N) { + *N = Nfailures; + return failures; +} + +int dvoverify_single (char *filename) { + + int isGood = TRUE; + + if (!VerifyTableFile (filename)) { + fprintf (stderr, "bad average table %s\n", filename); + isGood = FALSE; + } + + // change last 't' to 's': + int Nlast; + Nlast = strlen(filename) - 1; + filename[Nlast] = 's'; + if (!VerifyTableFile (filename)) { + fprintf (stderr, "bad secfilt table %s\n", filename); + isGood = FALSE; + } + + // change last 't' to 's': + Nlast = strlen(filename) - 1; + filename[Nlast] = 'm'; + if (!VerifyTableFile (filename)) { + fprintf (stderr, "bad measure table %s\n", filename); + isGood = FALSE; + } + return (isGood); +} + +// is this file a consistent FITS file? +int VerifyTableFile (char *filename) { + + int status, Next; + off_t Nbytes; + Header header; + + struct stat fileStats; + FILE *file; + + // does the file exist? + status = stat (filename, &fileStats); + if (status) { + // some error accessing the file. there is only one acceptable error: file not found + switch (errno) { + case ENOENT: + if (DEBUG) fprintf (stderr, "file does not exist, skipping %s\n", filename); + if (LIST_MISSING) return FALSE; + return TRUE; + case ENOMEM: + fprintf (stderr, "Out of memory: %s\n", filename); + return TRUE; + case EACCES: + fprintf (stderr, "Permission error on %s\n", filename); + return FALSE; + case EFAULT: + fprintf (stderr, "Bad address: %s\n", filename); + return FALSE; + case ELOOP: + fprintf (stderr, "Too many symbolic links encountered while traversing the path: %s\n", filename); + return FALSE; + case ENAMETOOLONG: + fprintf (stderr, "File name too long: %s\n", filename); + return FALSE; + case ENOTDIR: + fprintf (stderr, "A component of the path is not a directory: %s\n", filename); + return FALSE; + case EOVERFLOW: + fprintf (stderr, "file too large for program version: %s\n", filename); + return FALSE; + default: + fprintf (stderr, "unknown error: %s\n", filename); + return FALSE; + } + } + + // does it have any data? + if (fileStats.st_size == 0) { + fprintf (stderr, "file is empty: %s\n", filename); + return FALSE; + } + + // can we open it? + file = fopen(filename, "r"); + if (!file) { + fprintf (stderr, "unable to open valid file: %s\n", filename); + return FALSE; + } + + // scan all extentions + Nbytes = 0; + Next = -1; + if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); + while (Nbytes < fileStats.st_size) { + + // Check on the PHU + if (!gfits_fread_header (file, &header)) { + if (Next == -1) { + fprintf (stderr, "unable to read PHU header for %s\n", filename); + } else { + fprintf (stderr, "unable to read header for %s, extension %d (or file has excess bytes)\n", filename, Next); + } + fclose(file); + return (FALSE); + } + + // move to TBL header + Nbytes += header.datasize + gfits_data_size (&header); + if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); + if (Nbytes > fileStats.st_size) { + fprintf (stderr, "file is short ("OFF_T_FMT" vs "OFF_T_FMT"): %s\n", Nbytes, fileStats.st_size, filename); + gfits_free_header(&header); + fclose (file); + return FALSE; + } + gfits_free_header(&header); + + status = fseeko (file, Nbytes, SEEK_SET); + if (status) { + switch (errno) { + case EBADF: + fprintf (stderr, "something wrong with file handle: %s\n", filename); + fclose (file); + return FALSE; + case EINVAL: + fprintf (stderr, "invalid offset: %s\n", filename); + fclose (file); + return FALSE; + default: + fprintf (stderr, "other error in fseeko: %s\n", filename); + fclose (file); + return FALSE; + } + } + Next ++; + } + if (DEBUG) fprintf (stderr, "file is good: %s\n", filename); + fclose (file); + return TRUE; +} + +int CheckCatalogIndexes (char *filename, SkyRegion *region) { + + Catalog catalog; + int i, j, m, status; + + status = TRUE; + + // set the parameters which guide catalog open/load/create + catalog.filename = filename; + catalog.catformat = DVO_FORMAT_UNDEF; // read-only,do not set the catformat + catalog.catmode = DVO_MODE_UNDEF; // read-only, do not set the catmode + catalog.catflags = LOAD_AVES | LOAD_MEAS; + catalog.Nsecfilt = 0; + + // an error exit status here is a significant error (disk I/O or file access) + if (!dvo_catalog_open (&catalog, region, VERBOSE, "r")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + return FALSE; + } + + // Naves_disk == 0 implies an empty catalog file, skip empty catalogs + if (catalog.Naves_disk == 0) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + return TRUE; + } + + // if the table is SORTED, then the following can be checked + // check the following: + // measure[j].averef -> average[averef] + // measure[j].objID = average[averef].objID + // measure[j].catID = average[averef].catID + // measure[j].measureOffset < Nmeasure + // \sum average[].Nmeasure = Nmeasure + + // if the table is NOT SORTED, we have a subset of checks we can make + if (!catalog.sorted) { + fprintf (stderr, "!"); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + NNotSorted++; + if (VERBOSE) fprintf (stderr, "file is not sorted: %s\n", filename); + return TRUE; + } + + int NmeasureTotal = 0; + int measureOffsetOK = TRUE; + for (i = 0; i < catalog.Naverage; i++) { + NmeasureTotal += catalog.average[i].Nmeasure; + if (VERBOSE && !(NmeasureTotal <= catalog.Nmeasure)) { + fprintf (stderr, "NmeasureTotal > catalog.Nmeasure: %d %d %d\n", i, catalog.average[i].Nmeasure, (int) catalog.Nmeasure); + } + measureOffsetOK &= (catalog.average[i].measureOffset < catalog.Nmeasure); + if (VERBOSE && !(catalog.average[i].measureOffset < catalog.Nmeasure)) { + fprintf (stderr, "measureOffset >= catalog.Nmeasure: %d %d %d\n", i, catalog.average[i].measureOffset, (int) catalog.Nmeasure); + } + measureOffsetOK &= (catalog.average[i].measureOffset + catalog.average[i].Nmeasure <= catalog.Nmeasure); + if (VERBOSE && !(catalog.average[i].measureOffset + catalog.average[i].Nmeasure <= catalog.Nmeasure)) { + fprintf (stderr, "measureOffset + Nmeasure > catalog.Nmeasure : %d %d %d\n", i, catalog.average[i].Nmeasure, (int) catalog.Nmeasure); + } + } + + if (!measureOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid measureOffset\n", catalog.filename); + status = FALSE; + } + + if (NmeasureTotal != catalog.Nmeasure) { + fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog.filename); + status = FALSE; + } + + if (!status) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + return (status); + } + + int objIDsOK = TRUE; + int catIDsOK = TRUE; + int averefOK = TRUE; + + for (i = 0; i < catalog.Naverage; i++) { + m = catalog.average[i].measureOffset; + for (j = 0; j < catalog.average[i].Nmeasure; j++) { + objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID); + catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID); + averefOK &= (catalog.measure[m+j].averef == i); + } + } + + if (!objIDsOK) { + fprintf (stderr, "ERROR: catalog %s has invalid obj IDs\n", catalog.filename); + status = FALSE; + } + if (!catIDsOK) { + fprintf (stderr, "ERROR: catalog %s has invalid cat IDs\n", catalog.filename); + status = FALSE; + } + if (!averefOK) { + fprintf (stderr, "ERROR: catalog %s has invalid averef values\n", catalog.filename); + status = FALSE; + } + +// for (i = 0; i < catalog.Naverage; i++) { +// m = catalog.average[i].measureOffset; +// for (j = 0; i < catalog.Nmeasure; i++) { +// objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID); +// catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID); +// averefOK &= (catalog.measure[m+j].averef = i); +// } +// } + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + return status; +} + +// gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow); +// gfits_scan(&cpmHeaderTBL, "NAXIS2", "%d", 1, &Nrows); + + Index: trunk/Ohana/src/libautocode/Makefile.Targets =================================================================== --- trunk/Ohana/src/libautocode/Makefile.Targets (revision 34259) +++ trunk/Ohana/src/libautocode/Makefile.Targets (revision 34260) @@ -11,4 +11,5 @@ $(ASRC)/average-ps1-v2.$(ARCH).o \ $(ASRC)/average-ps1-v3.$(ARCH).o \ +$(ASRC)/average-ps1-v4.$(ARCH).o \ $(ASRC)/average-ps1-ref.$(ARCH).o \ $(ASRC)/secfilt.$(ARCH).o \ @@ -22,4 +23,5 @@ $(ASRC)/secfilt-ps1-v2.$(ARCH).o \ $(ASRC)/secfilt-ps1-v3.$(ARCH).o \ +$(ASRC)/secfilt-ps1-v4.$(ARCH).o \ $(ASRC)/secfilt-ps1-ref.$(ARCH).o \ $(ASRC)/measure.$(ARCH).o \ @@ -33,4 +35,5 @@ $(ASRC)/measure-ps1-v2.$(ARCH).o \ $(ASRC)/measure-ps1-v3.$(ARCH).o \ +$(ASRC)/measure-ps1-v4.$(ARCH).o \ $(ASRC)/measure-ps1-ref.$(ARCH).o \ $(ASRC)/missing.$(ARCH).o \ @@ -43,4 +46,5 @@ $(ASRC)/photcode-ps1-v2.$(ARCH).o \ $(ASRC)/photcode-ps1-v3.$(ARCH).o \ +$(ASRC)/photcode-ps1-v4.$(ARCH).o \ $(ASRC)/photcode-ps1-ref.$(ARCH).o \ $(ASRC)/image.$(ARCH).o \ @@ -55,4 +59,5 @@ $(ASRC)/image-ps1-v2.$(ARCH).o \ $(ASRC)/image-ps1-v3.$(ARCH).o \ +$(ASRC)/image-ps1-v4.$(ARCH).o \ $(ASRC)/image-ps1-ref.$(ARCH).o \ $(ASRC)/regimage.$(ARCH).o \ @@ -90,4 +95,5 @@ $(AINC)/average-ps1-v2.h \ $(AINC)/average-ps1-v3.h \ +$(AINC)/average-ps1-v4.h \ $(AINC)/average-ps1-ref.h \ $(AINC)/secfilt.h \ @@ -101,4 +107,5 @@ $(AINC)/secfilt-ps1-v2.h \ $(AINC)/secfilt-ps1-v3.h \ +$(AINC)/secfilt-ps1-v4.h \ $(AINC)/secfilt-ps1-ref.h \ $(AINC)/measure.h \ @@ -112,4 +119,5 @@ $(AINC)/measure-ps1-v2.h \ $(AINC)/measure-ps1-v3.h \ +$(AINC)/measure-ps1-v4.h \ $(AINC)/measure-ps1-ref.h \ $(AINC)/missing.h \ @@ -122,4 +130,5 @@ $(AINC)/photcode-ps1-v2.h \ $(AINC)/photcode-ps1-v3.h \ +$(AINC)/photcode-ps1-v4.h \ $(AINC)/photcode-ps1-ref.h \ $(AINC)/image.h \ @@ -134,4 +143,5 @@ $(AINC)/image-ps1-v2.h \ $(AINC)/image-ps1-v3.h \ +$(AINC)/image-ps1-v4.h \ $(AINC)/image-ps1-ref.h \ $(AINC)/regimage.h \ Index: trunk/Ohana/src/libautocode/def/average-ps1-v4.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-v4.d (revision 34260) +++ trunk/Ohana/src/libautocode/def/average-ps1-v4.d (revision 34260) @@ -0,0 +1,62 @@ +STRUCT Average_PS1_V4 +EXTNAME DVO_AVERAGE_PS1_V4 +TYPE BINTABLE +SIZE 120 +DESCRIPTION DVO Average Object Table + +# elements of data structure / FITS table + +FIELD R, RA, double, RA, decimal degrees +FIELD D, DEC, double, DEC, decimal degrees +FIELD dR, RA_ERR, float, RA error, arcsec +FIELD dD, DEC_ERR, float, DEC error, arcsec + +FIELD uR, U_RA, float, RA*cos(D) proper-motion, arcsec/year +FIELD uD, U_DEC, float, DEC proper-motion, arcsec/year +FIELD duR, V_RA_ERR, float, RA*cos(D) p-m error, arcsec/year +FIELD duD, V_DEC_ERR, float, DEC p-m error, arcsec/year +FIELD P, PAR, float, parallax, arcsec +FIELD dP, PAR_ERR, float, parallax error, arcsec + +FIELD ChiSqAve, CHISQ_POS, float, astrometry analysis chisq +FIELD ChiSqPM, CHISQ_PM, float, astrometry analysis chisq +FIELD ChiSqPar, CHISQ_PAP, float, astrometry analysis chisq +FIELD Tmean, MEAN_EPOCH, int, mean epoch (PM,PAR ref), unix time seconds +FIELD Trange, TIME_RANGE, int, mean epoch (PM,PAR ref), unix time seconds + +FIELD Xp, SIGMA_POS, float, position scatter, 1/100 arcsec +FIELD Npos, NUMBER_POS, unsigned short, number of detections used for astrometry + +# this limits us to a max of 64k measurements per object +FIELD Nmeasure, NMEASURE, unsigned short, number of psf measurements +FIELD Nmissing, NMISSING, unsigned short, number of missings +FIELD Nextend, NEXTEND, unsigned short, number of extended measurements + +FIELD measureOffset, OFF_MEASURE, uint32_t, offset to first psf measurement +FIELD missingOffset, OFF_MISSING, uint32_t, offset to first missing obs +FIELD extendOffset, OFF_EXTEND, uint32_t, offset to first extended measurement + +# 'flags' was called 'code' prior to 2009.02.07 +FIELD flags, FLAGS, uint32_t, average object flags (star; ghost; etc) +FIELD photFlagsUpper, PHOTFLAGS_U, uint32_t, upper bit of 2 bit summary of per-measure photflags +FIELD photFlagsLower, PHOTFLAGS_L, uint32_t, lower bit of 2 bit summary of per-measure photflags + +# objID + catID gives a unique ID for all objects in the database +FIELD objID, OBJ_ID, unsigned int, unique ID for object in table +FIELD catID, CAT_ID, unsigned int, unique ID for table in which object was first realized +FIELD extID, EXT_ID, uint64_t, external ID for object (eg PSPS objID) + +# 2*8 + 1*8 + 12*4 + 7*4 + 4*2 + +# photflagsUpper & photflagsLower: we have Nmeasures of a given source +# using psphot. each of these as a photFlag field, with bits +# describing the detection quality. photflagsUpper and photFlagsLower +# represent a 2 bit value (0-3) which defines the frequency of a given +# bit in the measurements: +# for a given bit, if that bit is raised for these percentile ranges, +# the following bits are set: +# min | max | L | U +# 0% | 25% | 0 | 0 +# 25% | 50% | 1 | 0 +# 50% | 75% | 0 | 1 +# 75% | 100% | 1 | 1 Index: trunk/Ohana/src/libautocode/def/cmf-ps1-sv1.d =================================================================== --- trunk/Ohana/src/libautocode/def/cmf-ps1-sv1.d (revision 34259) +++ trunk/Ohana/src/libautocode/def/cmf-ps1-sv1.d (revision 34260) @@ -4,4 +4,14 @@ TYPE BINTABLE SIZE 200 + +## note that this definition is inconsistent with the definition in +## pmModules/src/objects/pmSourceIO_CMF.c.in. The version in +## pmSourceIO*.c creates a table with data not properly aligned with +## the 8-byte boundaries. The structure defined by libautocode does +## this, but has a different order of elements (and adds padding2 to +## fix things). We have generated may files with PS1_SV1 as is, so +## I'll leave it. But in future a PS1_SV2 should be forced to match +## libautocode. Note that addstar knows to detect the alternate +## version of PS1_SV1 and correctly interpret its fields. # elements of data structure / FITS table Index: trunk/Ohana/src/libautocode/def/cmf-ps1-v1.d =================================================================== --- trunk/Ohana/src/libautocode/def/cmf-ps1-v1.d (revision 34259) +++ trunk/Ohana/src/libautocode/def/cmf-ps1-v1.d (revision 34260) @@ -44,35 +44,2 @@ # model shape parameters: F_major, F_minor, F_theta # centroid errors: sigma_X, sigma_Y, sigma_XY - -# IPP_IDET 1J IPP_IDET, unsigned i -# X_PSF 1E X_PSF, float, -# Y_PSF 1E Y_PSF, float, -# X_PSF_SIG 1E X_PSF_SIG, float, -# Y_PSF_SIG 1E Y_PSF_SIG, float, -# RA_PSF 1E RA_PSF, float, -# DEC_PSF 1E DEC_PSF, float, -# POSANGLE 1E POSANGLE, float, -# PLTSCALE 1E PLTSCALE, float, -# PSF_INST_MAG 1E PSF_INST_MAG, float, -# PSF_INST_MAG_SIG 1E PSF_INST_MAG_SIG, float, -# AP_MAG_STANDARD 1E AP_MAG_STANDARD, float, -# AP_MAG_RADIUS 1E -# PEAK_FLUX_AS_MAG 1E PEAK_FLUX_AS_MAG, float, -# CAL_PSF_MAG 1E CAL_PSF_MAG, float, -# CAL_PSF_MAG_SIG 1E CAL_PSF_MAG_SIG, float, -# SKY 1E SKY, float, -# SKY_SIGMA 1E SKY_SIG, float, -# PSF_CHISQ 1E PSF_CHISQ, float, -# CR_NSIGMA 1E CR_NSIGMA, float, -# EXT_NSIGMA 1E EXT_NSIGMA, float, -# PSF_MAJOR 1E PSF_MAJOR, float, -# PSF_MINOR 1E PSF_MINOR, float, -# PSF_THETA 1E PSF_THETA, float, -# PSF_QF 1E PSF_QF, float, -# PSF_NDOF 1J PSF_NDOF, int, -# PSF_NPIX 1J PSF_NPIX, int, -# MOMENTS_XX 1E MOMENTS_XX, float, -# MOMENTS_XY 1E MOMENTS_XY, float, -# MOMENTS_YY 1E MOMENTS_YY, float, -# N_FRAMES 1I N_FRAMES, short, -# FLAGS 1J FLAGS, int, Index: trunk/Ohana/src/libautocode/def/cmf-ps1-v2.d =================================================================== --- trunk/Ohana/src/libautocode/def/cmf-ps1-v2.d (revision 34259) +++ trunk/Ohana/src/libautocode/def/cmf-ps1-v2.d (revision 34260) @@ -40,39 +40,4 @@ FIELD padding, PADDING, short, padding for 8byte records -# for an object in an image, we have three triplets that tell us about the shape: -# second moments: Mxx, Mxy, Myy -# model shape parameters: F_major, F_minor, F_theta -# centroid errors: sigma_X, sigma_Y, sigma_XY - -# IPP_IDET 1J IPP_IDET, unsigned i -# X_PSF 1E X_PSF, float, -# Y_PSF 1E Y_PSF, float, -# X_PSF_SIG 1E X_PSF_SIG, float, -# Y_PSF_SIG 1E Y_PSF_SIG, float, -# RA_PSF 1E RA_PSF, float, -# DEC_PSF 1E DEC_PSF, float, -# POSANGLE 1E POSANGLE, float, -# PLTSCALE 1E PLTSCALE, float, -# PSF_INST_MAG 1E PSF_INST_MAG, float, -# PSF_INST_MAG_SIG 1E PSF_INST_MAG_SIG, float, -# AP_MAG_STANDARD 1E AP_MAG_STANDARD, float, -# AP_MAG_RADIUS 1E -# PEAK_FLUX_AS_MAG 1E PEAK_FLUX_AS_MAG, float, -# CAL_PSF_MAG 1E CAL_PSF_MAG, float, -# CAL_PSF_MAG_SIG 1E CAL_PSF_MAG_SIG, float, -# SKY 1E SKY, float, -# SKY_SIGMA 1E SKY_SIG, float, -# PSF_CHISQ 1E PSF_CHISQ, float, -# CR_NSIGMA 1E CR_NSIGMA, float, -# EXT_NSIGMA 1E EXT_NSIGMA, float, -# PSF_MAJOR 1E PSF_MAJOR, float, -# PSF_MINOR 1E PSF_MINOR, float, -# PSF_THETA 1E PSF_THETA, float, -# PSF_QF 1E PSF_QF, float, -# PSF_NDOF 1J PSF_NDOF, int, -# PSF_NPIX 1J PSF_NPIX, int, -# MOMENTS_XX 1E MOMENTS_XX, float, -# MOMENTS_XY 1E MOMENTS_XY, float, -# MOMENTS_YY 1E MOMENTS_YY, float, -# N_FRAMES 1I N_FRAMES, short, -# FLAGS 1J FLAGS, int, +# fields added since PS1_V1: +# RA & DEC promoted to double Index: trunk/Ohana/src/libautocode/def/cmf-ps1-v3.d =================================================================== --- trunk/Ohana/src/libautocode/def/cmf-ps1-v3.d (revision 34259) +++ trunk/Ohana/src/libautocode/def/cmf-ps1-v3.d (revision 34260) @@ -64,54 +64,18 @@ # centroid errors: sigma_X, sigma_Y, sigma_XY -# # elements of data structure / FITS table -# FIELD detID, IPP_IDET, IPP_IDET 1J unsigned int, detection ID -# FIELD X, X_PSF, X_PSF 1E float, x coord, pixels -# FIELD Y, Y_PSF, Y_PSF 1E float, y coord, pixels -# FIELD dX, X_PSF_SIG, X_PSF_SIG 1E float, x coord error, pixels -# FIELD dY, Y_PSF_SIG, Y_PSF_SIG 1E float, y coord error, pixels -# FIELD posangle, POSANGLE, POSANGLE 1E float, Posangle at source, degrees -# FIELD pltscale, PLTSCALE, PLTSCALE 1E float, Plate Scale at source, arcsec/pixel -# FIELD M, PSF_INST_MAG, PSF_INST_MAG 1E float, inst mags, mags -# FIELD dM, PSF_INST_MAG_SIG, PSF_INST_MAG_SIG 1E float, inst mag error, mags -# FIELD flux, PSF_INST_FLUX, PSF_INST_FLUX 1E float, psf flux, counts -# FIELD flux, PSF_INST_FLUX_SIG, PSF_INST_FLUX_SIG 1E float, psf flux error, counts -# FIELD Map, AP_MAG_STANDARD, AP_MAG 1E float, standard aperture mag, mags -# FIELD Map, AP_MAG_RAW, AP_MAG_RAW 1E float, raw aperture mag, mags -# FIELD apRadius, AP_MAG_RADIUS, AP_MAG_RADIUS 1E float, radius used for fit, pixels -# FIELD Mpeak, PEAK_FLUX_AS_MAG, PEAK_FLUX_AS_MAG 1E float, peak flux as a mag, mags -# FIELD Mcalib, CAL_PSF_MAG, CAL_PSF_MAG 1E float, calibrated psf mag, mags -# FIELD dMcal, CAL_PSF_MAG_SIG, CAL_PSF_MAG_SIG 1E float, zero point scatter, mags -# FIELD RA, RA_PSF, RA_PSF 1D double, PSF RA coord, degrees -# FIELD DEC, DEC_PSF, DEC_PSF 1D double, PSF DEC coord, degrees -# FIELD sky, SKY, SKY 1E float, sky flux, cnts/sec -# FIELD dSky, SKY_SIG, SKY_SIGMA 1E float, sky flux error, cnts/sec -# FIELD psfChisq, PSF_CHISQ, PSF_CHISQ 1E float, psf fit chisq -# FIELD crNsigma, CR_NSIGMA, CR_NSIGMA 1E float, Nsigma deviations from PSF to CF -# FIELD extNsigma, EXT_NSIGMA, EXT_NSIGMA 1E float, Nsigma deviations from PSF to EXT -# FIELD fx, PSF_MAJOR, PSF_MAJOR 1E float, psf fit major axis, pixels -# FIELD fy, PSF_MINOR, PSF_MINOR 1E float, psf fit minor axis, pixels -# FIELD df, PSF_THETA, PSF_THETA 1E float, ellipse angle, degrees -# FIELD psfQual, PSF_QF, PSF_QF 1E float, quality factor -# FIELD psfQualPerfect, PSF_QF_PERFECT, PSF_QF_PERFECT 1E float, quality factor perfect -# FIELD psfNdof, PSF_NDOF, PSF_NDOF 1J int, psf degrees of freedom -# FIELD psfNpix, PSF_NPIX, PSF_NPIX 1J int, psf number of pixels -# FIELD Mxx, MOMENTS_XX, MOMENTS_XX 1E float, second moment X, pixels^2 -# FIELD Mxy, MOMENTS_XY, MOMENTS_XY 1E float, second moment Y, pixels^2 -# FIELD Myy, MOMENTS_YY, MOMENTS_YY 1E float, second moment XY, pixels^2 -# FIELD M3c, MOMENTS_M3C, MOMENTS_M3C 1E float, third moment cos(t), pixels^3 -# FIELD M3s, MOMENTS_M3S, MOMENTS_M3S 1E float, third moment sin(t), pixels^3 -# FIELD M4c, MOMENTS_M4C, MOMENTS_M4C 1E float, fourth moment cos(t), pixels^4 -# FIELD M4s, MOMENTS_M4S, MOMENTS_M4S 1E float, fourth moment sin(t), pixels^4 -# FIELD Mr1, MOMENTS_R1, MOMENTS_R1 1E float, first radial moment, pixels -# FIELD Mrh, MOMENTS_RH, MOMENTS_RH 1E float, half radial moment, pixels^1/2 -# FIELD kronFlux, KRON_FLUX, KRON_FLUX 1E float, kron flux, counts -# FIELD kronFluxErr, KRON_FLUX_ERR, KRON_FLUX_ERR 1E float, kron flux error, counts -# FIELD kronInner, KRON_FLUX_INNER, KRON_FLUX_INNER 1E float, kron flux 1DEC_origin); + gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset); + + ftable.header = &theader; + + /*** zone information table ***/ + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN_RAW(tree->Nband, "NBAND", int); + GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); + tree->Nzone = Nrow; + + // allocate the storage arrays + ALLOCATE (tree->ra, double *, tree->Nzone); + ALLOCATE (tree->dec, double *, tree->Nzone); + ALLOCATE (tree->cell, int *, tree->Nzone); + ALLOCATE (tree->name, char **, tree->Nzone); + for (i = 0; i < tree->Nzone; i++) { + ALLOCATE (tree->ra[i], double, tree->Nband[i]); + ALLOCATE (tree->dec[i], double, tree->Nband[i]); + ALLOCATE (tree->cell[i], int, tree->Nband[i]); + ALLOCATE (tree->name[i], char *, tree->Nband[i]); + for (j = 0; j < tree->Nband[i]; j++) { + ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH); + } + } + } + + /*** cell information table ***/ + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN_NEW(R, "RA", double); + GET_COLUMN_NEW(D, "DEC", double); + GET_COLUMN_NEW(zone, "ZONE", int); + GET_COLUMN_NEW(band, "BAND", int); + GET_COLUMN_NEW(index, "INDEX", int); + GET_COLUMN_NEW(name, "NAME", char); // XXX how is this done? + gfits_free_header (&theader); + gfits_free_table (&ftable); + fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow); + + // assign the storage arrays + for (i = 0; i < Nrow; i++) { + nz = zone[i]; + nb = band[i]; + tree->ra[nz][nb] = R[i]; + tree->dec[nz][nb] = D[i]; + tree->cell[nz][nb] = i; // XXX ? + memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); + } + + free (R ); + free (D ); + free (zone ); + free (band ); + free (index ); + free (name ); + } + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return tree; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (tree) free (tree); + + fclose (f); + return NULL; +} + +// we are passed a BoundaryTree structure, write it to a FITS table (3 ext) +int BoundaryTreeSave(char *filename, BoundaryTree *tree) { + + int i, nz, nb; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open boundary tree file for output %s\n", filename); + return FALSE; + } + + // we need some information in the header to define the layout + gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin); + gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset); + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + /*** zone information table ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "ZONE_DATA"); + + gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NBAND", "number of cells in each zone", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA_ORIGIN", "origin of ra cell sequence", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA_OFFSET", "offset per cell of ra cell sequence", "degree/cell", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + int *zone = NULL; ALLOCATE (zone, int, tree->Nzone); + + // assign the storage arrays + for (i = 0; i < tree->Nzone; i++) { + zone[i] = i; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "NBAND", tree->Nband, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN", tree->RA_origin, tree->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET", tree->RA_offset, tree->Nzone); + free (zone); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** cell information table ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "CELL_DATA"); + + char fmt[16]; + snprintf (fmt, 16, "%dA", BOUNDARY_TREE_NAME_LENGTH); + gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000) of cell center", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + int Ncell = 0; + for (i = 0; i < tree->Nzone; i++) { + Ncell += tree->Nband[i]; + } + + // create intermediate storage arrays + // NOTE: we have to unroll the 2D arrays in tree into 1D arrays + double *R ; ALLOCATE (R, double, Ncell); + double *D ; ALLOCATE (D, double, Ncell); + int *zone ; ALLOCATE (zone, int, Ncell); + int *band ; ALLOCATE (band, int, Ncell); + int *index ; ALLOCATE (index, int, Ncell); + char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TREE_NAME_LENGTH); + + // NOTE: a table column of characters must be fixed width, and is passed as a + // contiguous array of Nchar * Nrow values + + // assign the storage arrays + i = 0; + for (nz = 0; nz < tree->Nzone; nz++) { + for (nb = 0; nb < tree->Nband[nz]; nb++) { + R[i] = tree->ra[nz][nb]; + D[i] = tree->dec[nz][nb]; + zone[i] = nz; + band[i] = nb; + index[i] = i; // or tree->cells[nz][nb] ? + memcpy(&name[i*BOUNDARY_TREE_NAME_LENGTH], tree->name[nz][nb], BOUNDARY_TREE_NAME_LENGTH); + i++; + } + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Ncell); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Ncell); + gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, Ncell); + gfits_set_bintable_column (&theader, &ftable, "BAND", band, Ncell); + gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell); + gfits_set_bintable_column (&theader, &ftable, "NAME", name, Ncell); + + free (R ); + free (D ); + free (zone ); + free (band ); + free (index ); + free (name ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +// the boundary tree... +// given an (ra,dec) pair, find the containing projection cell. + +int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { + + // first, find the containing zone + + // if we know dDEC, we can get the bin instantly: + *zone = (dec - tree->DEC_origin) / tree->DEC_offset; + + if (*zone < 0) return FALSE; + if (*zone >= tree->Nzone) return FALSE; + + // now select the RA bin for that zone + *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone]; + + if (*band < 0) return FALSE; + if (*band >= tree->Nband[*zone]) *band = 0; + + return TRUE; +} + Index: trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c =================================================================== --- trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c (revision 34260) @@ -63,4 +63,5 @@ CONVERT_FORMAT("DVO_PHOTCODE_PS1_V2", PS1_V2); CONVERT_FORMAT("DVO_PHOTCODE_PS1_V3", PS1_V3); + CONVERT_FORMAT("DVO_PHOTCODE_PS1_V4", PS1_V4); table = GetPhotcodeTable (); Index: trunk/Ohana/src/libdvo/src/SavePhotcodesFITS.c =================================================================== --- trunk/Ohana/src/libdvo/src/SavePhotcodesFITS.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/SavePhotcodesFITS.c (revision 34260) @@ -4,4 +4,5 @@ /* locking is used to avoid collisions with programs trying to update the photcodes values */ /* XXX better distinction between NOT FOUND and FAILURE */ + int SavePhotcodesFITS (char *filename) { @@ -29,9 +30,9 @@ // for the moment, we simply support the latest photcode format for output // XXX update this as needed as new formats are defined - PhotCode_PS1_V3 *photcode_output = PhotCode_Internal_To_PS1_V3 (table[0].code, table[0].Ncode); + PhotCode_PS1_V4 *photcode_output = PhotCode_Internal_To_PS1_V4 (table[0].code, table[0].Ncode); /* convert FITS format data to internal format (byteswaps & EXTNAME) */ if (!gfits_db_create (&db)) return (FALSE); - if (!gfits_table_set_PhotCode_PS1_V3 (&db.ftable, photcode_output, table[0].Ncode)) return (FALSE); + if (!gfits_table_set_PhotCode_PS1_V4 (&db.ftable, photcode_output, table[0].Ncode)) return (FALSE); if (!gfits_db_save (&db)) return (FALSE); if (!gfits_db_close (&db)) return (FALSE); Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 34260) @@ -233,4 +233,8 @@ break; + case MAG_KRON: + value.Flt = PhotKronAve (field->photcode, average, secfilt); + break; + case MAG_20: value.Flt = PhotM20 (field->photcode, average, secfilt); @@ -241,4 +245,20 @@ case MAG_UC_DIST: value.Flt = PhotUCdist (field->photcode, average, secfilt); + break; + case MAG_STACK_ID: + value.Int = PhotStackID (field->photcode, average, secfilt); + break; + + case MAG_FLUX_PSF: + value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt); + break; + case MAG_FLUX_PSF_ERR: + value.Flt = PhotAvedFluxPSF (field->photcode, average, secfilt); + break; + case MAG_FLUX_KRON: + value.Flt = PhotAveFluxKron (field->photcode, average, secfilt); + break; + case MAG_FLUX_KRON_ERR: + value.Flt = PhotAvedFluxKron (field->photcode, average, secfilt); break; } Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 34260) @@ -117,5 +117,17 @@ break; case MAG_APER: - value.Flt = PhotAper (measure); + value.Flt = PhotAper (measure); + break; + case MAG_APER_INST: + value.Flt = PhotAperInst (measure); + break; + case MAG_KRON: + value.Flt = PhotKron (measure); + break; + case MAG_KRON_INST: + value.Flt = PhotKronInst (measure); + break; + case MAG_KRON_ERR: + value.Flt = measure[0].dMkron; break; case MAG_ERR: @@ -152,4 +164,16 @@ if (Nsec == -1) break; value.Int = secfilt[Nsec].Nused; + break; + case MAG_FLUX_PSF: + value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt); + break; + case MAG_FLUX_PSF_ERR: + value.Flt = PhotAvedFluxPSF (field->photcode, average, secfilt); + break; + case MAG_FLUX_KRON: + value.Flt = PhotAveFluxKron (field->photcode, average, secfilt); + break; + case MAG_FLUX_KRON_ERR: + value.Flt = PhotAvedFluxKron (field->photcode, average, secfilt); break; } @@ -498,5 +522,33 @@ break; - // add the star/galaxy sep -- + case MEAS_CENTER_OFFSET: /* OK */ + { + Image *image; + image = MatchImageDVO (measure[0].t, measure[0].photcode, measure[0].imageID); + if (image == NULL) break; + + // we have measure[0].Xccd,Yccd and image[0].NX,NY. Find the distance to the center + + // XXX we may hypotetically have images with -NX to +NX here (eg, projection center), but + // we do not get a detection from that type of image + float Xcenter = 0.5*image[0].NX; + float Ycenter = 0.5*image[0].NY; + float distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter); + value.Flt = distance; + } + break; + + case MEAS_FLUX_PSF: /* OK */ + value.Flt = measure[0].FluxPSF; + break; + case MEAS_FLUX_PSF_ERR: /* OK */ + value.Flt = measure[0].dFluxPSF; + break; + case MEAS_FLUX_KRON: /* OK */ + value.Flt = measure[0].FluxKron; + break; + case MEAS_FLUX_KRON_ERR: /* OK */ + value.Flt = measure[0].dFluxKron; + break; } return (value); Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 34260) @@ -30,24 +30,35 @@ int GetMagMode (char *string) { - if (!strcasecmp (string, "inst")) return (MAG_INST); - if (!strcasecmp (string, "cat")) return (MAG_CAT); - if (!strcasecmp (string, "sys")) return (MAG_SYS); - if (!strcasecmp (string, "rel")) return (MAG_REL); - if (!strcasecmp (string, "cal")) return (MAG_CAL); - if (!strcasecmp (string, "ave")) return (MAG_AVE); - if (!strcasecmp (string, "ref")) return (MAG_REF); - if (!strcasecmp (string, "ap")) return (MAG_APER); - if (!strcasecmp (string, "aper")) return (MAG_APER); - if (!strcasecmp (string, "err")) return (MAG_ERR); - if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); - if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); - if (!strcasecmp (string, "flags")) return (MAG_PHOT_FLAGS); - if (!strcasecmp (string, "chisq")) return (MAG_CHISQ); - if (!strcasecmp (string, "ncode")) return (MAG_NCODE); - if (!strcasecmp (string, "nphot")) return (MAG_NPHOT); - if (!strcasecmp (string, "stdev")) return (MAG_STDEV); - if (!strcasecmp (string, "20")) return (MAG_20); - if (!strcasecmp (string, "80")) return (MAG_80); - if (!strcasecmp (string, "ucdist")) return (MAG_UC_DIST); + if (!strcasecmp (string, "inst")) return (MAG_INST); + if (!strcasecmp (string, "cat")) return (MAG_CAT); + if (!strcasecmp (string, "sys")) return (MAG_SYS); + if (!strcasecmp (string, "rel")) return (MAG_REL); + if (!strcasecmp (string, "cal")) return (MAG_CAL); + if (!strcasecmp (string, "ave")) return (MAG_AVE); + if (!strcasecmp (string, "ref")) return (MAG_REF); + if (!strcasecmp (string, "ap")) return (MAG_APER); + if (!strcasecmp (string, "aper")) return (MAG_APER); + if (!strcasecmp (string, "aperinst")) return (MAG_APER_INST); + if (!strcasecmp (string, "aper_inst")) return (MAG_APER_INST); + if (!strcasecmp (string, "kron")) return (MAG_KRON); + if (!strcasecmp (string, "kroninst")) return (MAG_KRON_INST); + if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); + if (!strcasecmp (string, "kronerr")) return (MAG_KRON_ERR); + if (!strcasecmp (string, "err")) return (MAG_ERR); + if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); + if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); + if (!strcasecmp (string, "flags")) return (MAG_PHOT_FLAGS); + if (!strcasecmp (string, "chisq")) return (MAG_CHISQ); + if (!strcasecmp (string, "ncode")) return (MAG_NCODE); + if (!strcasecmp (string, "nphot")) return (MAG_NPHOT); + if (!strcasecmp (string, "stdev")) return (MAG_STDEV); + if (!strcasecmp (string, "20")) return (MAG_20); + if (!strcasecmp (string, "80")) return (MAG_80); + if (!strcasecmp (string, "ucdist")) return (MAG_UC_DIST); + if (!strcasecmp (string, "stackID")) return (MAG_STACK_ID); + if (!strcasecmp (string, "fluxpsf")) return (MAG_FLUX_PSF); + if (!strcasecmp (string, "fluxpsferr")) return (MAG_FLUX_PSF_ERR); + if (!strcasecmp (string, "fluxkron")) return (MAG_FLUX_KRON); + if (!strcasecmp (string, "fluxkronerr")) return (MAG_FLUX_KRON_ERR); return (MAG_NONE); } @@ -229,4 +240,11 @@ if (!strcasecmp (fieldName, "MCAL_OFFSET")) ESCAPE (MEAS_MCAL_OFFSET, MAG_NONE, OPIHI_FLT); if (!strcasecmp (fieldName, "FLAT")) ESCAPE (MEAS_FLAT, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "CENTER_OFFSET")) ESCAPE (MEAS_CENTER_OFFSET, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX")) ESCAPE (MEAS_FLUX_PSF, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX_ERR")) ESCAPE (MEAS_FLUX_PSF_ERR, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX_PSF")) ESCAPE (MEAS_FLUX_PSF, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX_PSF_ERR")) ESCAPE (MEAS_FLUX_PSF_ERR, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX_KRON")) ESCAPE (MEAS_FLUX_KRON, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "FLUX_KRON_ERR")) ESCAPE (MEAS_FLUX_KRON_ERR, MAG_NONE, OPIHI_FLT); // for words that don't parse, try a photcode @@ -242,4 +260,5 @@ field->magMode = mode; switch (mode) { + case MAG_STACK_ID: case MAG_NCODE: case MAG_NPHOT: @@ -336,4 +355,5 @@ field->magMode = mode; switch (mode) { + case MAG_STACK_ID: case MAG_NCODE: case MAG_NPHOT: Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 34260) @@ -55,4 +55,5 @@ if (!strcasecmp (catformat, "PS1_V2")) return (DVO_FORMAT_PS1_V2); if (!strcasecmp (catformat, "PS1_V3")) return (DVO_FORMAT_PS1_V3); + if (!strcasecmp (catformat, "PS1_V4")) return (DVO_FORMAT_PS1_V4); if (!strcasecmp (catformat, "PS1_REF")) return (DVO_FORMAT_PS1_REF); return (DVO_FORMAT_UNDEF); @@ -66,4 +67,136 @@ if (!strcasecmp (catmode, "SPLIT")) return (DVO_MODE_SPLIT); return (DVO_MODE_UNDEF); +} + +// init all data, or just catalog data +void dvo_average_init (Average *average) { + average->R = 0; + average->D = 0; + average->dR = 0; + average->dD = 0; + + average->uR = 0; + average->uD = 0; + average->duR = 0; + average->duD = 0; + average->P = 0; + average->dP = 0; + + average->ChiSqAve = 0.0; + average->ChiSqPM = 0.0; + average->ChiSqPar = 0.0; + average->Tmean = 0; + average->Trange = 0; + + average->Xp = 0; + average->Npos = 0; + + average->Nmeasure = 0; + average->Nmissing = 0; + average->Nextend = 0; + average->measureOffset = -1; + average->missingOffset = -1; + average->extendOffset = -1; + + average->flags = 0; + average->photFlagsUpper = 0; + average->photFlagsLower = 0; + + average->objID = 0; + average->catID = 0; + average->extID = 0; +} + +// init all data, or just catalog data +void dvo_secfilt_init (SecFilt *secfilt) { + secfilt->M = NAN; + secfilt->Map = NAN; + secfilt->Mkron = NAN; + secfilt->dMkron = NAN; + secfilt->dM = NAN; + secfilt->Xm = NAN_S_SHORT; + + secfilt->FluxPSF = NAN; + secfilt->dFluxPSF = NAN; + secfilt->FluxKron = NAN; + secfilt->dFluxKron = NAN; + + secfilt->flags = 0; + secfilt->Ncode = 0; + secfilt->Nused = 0; + + secfilt->M_20 = NAN_S_SHORT; + secfilt->M_80 = NAN_S_SHORT; + + secfilt->ubercalDist = 1000; + secfilt->Mstdev = NAN_S_SHORT; + secfilt->stackID = 0; + + secfilt->dummy = 0; +} + +// init all data, or just catalog data +void dvo_measure_init (Measure *measure) { + measure->M = NAN; + measure->dR = NAN; + measure->dD = NAN; + measure->M = NAN; + measure->Mcal = NAN; + measure->Map = NAN; + measure->Mkron = NAN; + measure->dMkron = NAN; + measure->dM = NAN; + measure->dMcal = NAN; + measure->dt = NAN; + + measure->FluxPSF = NAN; + measure->dFluxPSF = NAN; + measure->FluxKron = NAN; + measure->dFluxKron = NAN; + + measure->airmass = NAN; + measure->az = NAN; + + measure->Xccd = NAN; + measure->Yccd = NAN; + + measure->Sky = NAN; + measure->dSky = NAN; + + measure->t = 0; + measure->t_msec = 0; + measure->averef = 0; + + measure->detID = 0; + measure->imageID = 0; + measure->objID = 0; + measure->catID = 0; + measure->extID = 0; + + measure->psfQual = NAN; + measure->psfChisq = NAN; + measure->psfNdof = 0; + measure->psfNpix = 0; + measure->crNsigma = NAN; + measure->extNsigma = NAN; + + measure->FWx = 0; + measure->FWy = 0; + measure->theta = 0; + + measure->Mxx = 0; + measure->Mxy = 0; + measure->Myy = 0; + + measure->dXccd = 0; + measure->dYccd = 0; + measure->dRsys = 0; + + measure->posangle = 0; + measure->pltscale = NAN; + + measure->photcode = 0; + measure->dbFlags = 0; + measure->photFlags = 0; } @@ -600,2 +733,3 @@ } + Index: trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 34260) @@ -89,4 +89,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -282,4 +283,5 @@ if (catalog[0].catformat == DVO_FORMAT_PS1_V2) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_V2"); if (catalog[0].catformat == DVO_FORMAT_PS1_V3) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_V3"); + if (catalog[0].catformat == DVO_FORMAT_PS1_V4) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_V4"); if (catalog[0].catformat == DVO_FORMAT_PS1_REF) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_REF"); @@ -384,4 +386,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -434,4 +437,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -488,4 +492,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -538,4 +543,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -592,4 +598,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -642,4 +649,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); Index: trunk/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 34260) @@ -40,4 +40,5 @@ CONVERT_FORMAT ("DVO_AVERAGE_PS1_V2", PS1_V2); CONVERT_FORMAT ("DVO_AVERAGE_PS1_V3", PS1_V3); + CONVERT_FORMAT ("DVO_AVERAGE_PS1_V4", PS1_V4); CONVERT_FORMAT ("DVO_AVERAGE_PS1_REF", PS1_REF); # undef CONVERT_FORMAT @@ -96,4 +97,5 @@ CONVERT_FORMAT ("DVO_AVERAGE_PS1_V2", PS1_V2, PS1_V2); CONVERT_FORMAT ("DVO_AVERAGE_PS1_V3", PS1_V3, PS1_V3); + CONVERT_FORMAT ("DVO_AVERAGE_PS1_V4", PS1_V4, PS1_V4); CONVERT_FORMAT ("DVO_AVERAGE_PS1_REF", PS1_REF, PS1_REF); # undef CONVERT_FORMAT @@ -130,4 +132,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); # undef FORMAT_CASE @@ -185,4 +188,5 @@ CONVERT_FORMAT ("DVO_MEASURE_PS1_V2", PS1_V2, PS1_V2); CONVERT_FORMAT ("DVO_MEASURE_PS1_V3", PS1_V3, PS1_V3); + CONVERT_FORMAT ("DVO_MEASURE_PS1_V4", PS1_V4, PS1_V4); CONVERT_FORMAT ("DVO_MEASURE_PS1_REF", PS1_REF, PS1_REF); # undef CONVERT_FORMAT @@ -219,4 +223,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); # undef FORMAT_CASE @@ -274,4 +279,5 @@ CONVERT_FORMAT ("DVO_SECFILT_PS1_V2", PS1_V2, PS1_V2); CONVERT_FORMAT ("DVO_SECFILT_PS1_V3", PS1_V3, PS1_V3); + CONVERT_FORMAT ("DVO_SECFILT_PS1_V4", PS1_V4, PS1_V4); CONVERT_FORMAT ("DVO_SECFILT_PS1_REF", PS1_REF, PS1_REF); # undef CONVERT_FORMAT @@ -308,4 +314,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); # undef FORMAT_CASE @@ -387,4 +394,5 @@ CONVERT_FORMAT ("DVO_IMAGE_PS1_V2", PS1_V2, PS1_V2); CONVERT_FORMAT ("DVO_IMAGE_PS1_V3", PS1_V3, PS1_V3); + CONVERT_FORMAT ("DVO_IMAGE_PS1_V4", PS1_V4, PS1_V4); CONVERT_FORMAT ("DVO_IMAGE_PS1_REF", PS1_REF, PS1_REF); @@ -425,4 +433,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); @@ -479,4 +488,5 @@ FORMAT_CASE (PS1_V2, PS1_V2); FORMAT_CASE (PS1_V3, PS1_V3); + FORMAT_CASE (PS1_V4, PS1_V4); FORMAT_CASE (PS1_REF, PS1_REF); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c (revision 34260) @@ -57,4 +57,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -215,4 +223,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 34260) @@ -54,4 +54,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -209,4 +217,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c (revision 34260) @@ -52,4 +52,12 @@ out[i].dbFlags = 0; out[i].photFlags = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -169,4 +177,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 34260) @@ -52,4 +52,12 @@ out[i].dbFlags = in[i].dbFlags; out[i].photFlags = in[i].photFlags; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -213,4 +221,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 34260) @@ -53,4 +53,12 @@ out[i].dbFlags = in[i].dbFlags; out[i].photFlags = in[i].photFlags; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -216,4 +224,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 34260) @@ -53,4 +53,12 @@ out[i].dbFlags = in[i].dbFlags; out[i].photFlags = in[i].photFlags; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -214,4 +222,13 @@ out[i].Mstdev = in[i].Mstdev; out[i].ubercalDist = in[i].ubercalDist; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 34260) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 34260) @@ -0,0 +1,463 @@ +# include + +/* convert PS1_V4 formats to internal formats */ + +Measure *Measure_PS1_V4_ToInternal (Measure_PS1_V4 *in, off_t Nvalues) { + + off_t i; + Measure *out; + + ALLOCATE_ZERO (out, Measure, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].M = in[i].M; + out[i].Mcal = in[i].Mcal; + out[i].Map = in[i].Map; + out[i].Mkron = in[i].Mkron; + out[i].dMkron = in[i].dMkron; + out[i].dM = in[i].dM; + out[i].dMcal = in[i].dMcal; + out[i].dt = in[i].dt; + out[i].FluxPSF = in[i].FluxPSF; + out[i].dFluxPSF = in[i].dFluxPSF; + out[i].FluxKron = in[i].FluxKron; + out[i].dFluxKron = in[i].dFluxKron; + out[i].airmass = in[i].airmass; + out[i].az = in[i].az; + out[i].Xccd = in[i].Xccd; + out[i].Yccd = in[i].Yccd; + out[i].Sky = in[i].Sky; + out[i].dSky = in[i].dSky; + out[i].t = in[i].t; + out[i].t_msec = in[i].t_msec; + out[i].averef = in[i].averef; + out[i].detID = in[i].detID; + out[i].imageID = in[i].imageID; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].extID = in[i].extID; + out[i].psfQual = in[i].psfQual; + out[i].psfChisq = in[i].psfChisq; + out[i].psfNdof = in[i].psfNdof; + out[i].psfNpix = in[i].psfNpix; + out[i].crNsigma = in[i].crNsigma; + out[i].extNsigma = in[i].extNsigma; + out[i].FWx = in[i].FWx; + out[i].FWy = in[i].FWy; + out[i].theta = in[i].theta; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].dXccd = in[i].dXccd; + out[i].dYccd = in[i].dYccd; + out[i].dRsys = in[i].dRsys; + out[i].posangle = in[i].posangle; + out[i].pltscale = in[i].pltscale; + out[i].photcode = in[i].photcode; + out[i].dbFlags = in[i].dbFlags; + out[i].photFlags = in[i].photFlags; + } + return (out); +} + +Measure_PS1_V4 *MeasureInternalTo_PS1_V4 (Measure *in, off_t Nvalues) { + + off_t i; + Measure_PS1_V4 *out; + + ALLOCATE_ZERO (out, Measure_PS1_V4, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].M = in[i].M; + out[i].Mcal = in[i].Mcal; + out[i].Map = in[i].Map; + out[i].Mkron = in[i].Mkron; + out[i].dMkron = in[i].dMkron; + out[i].dM = in[i].dM; + out[i].dMcal = in[i].dMcal; + out[i].dt = in[i].dt; + out[i].FluxPSF = in[i].FluxPSF; + out[i].dFluxPSF = in[i].dFluxPSF; + out[i].FluxKron = in[i].FluxKron; + out[i].dFluxKron = in[i].dFluxKron; + out[i].airmass = in[i].airmass; + out[i].az = in[i].az; + out[i].Xccd = in[i].Xccd; + out[i].Yccd = in[i].Yccd; + out[i].Sky = in[i].Sky; + out[i].dSky = in[i].dSky; + out[i].t = in[i].t; + out[i].t_msec = in[i].t_msec; + out[i].averef = in[i].averef; + out[i].detID = in[i].detID; + out[i].imageID = in[i].imageID; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].extID = in[i].extID; + out[i].psfQual = in[i].psfQual; + out[i].psfChisq = in[i].psfChisq; + out[i].psfNdof = in[i].psfNdof; + out[i].psfNpix = in[i].psfNpix; + out[i].crNsigma = in[i].crNsigma; + out[i].extNsigma = in[i].extNsigma; + out[i].FWx = in[i].FWx; + out[i].FWy = in[i].FWy; + out[i].theta = in[i].theta; + out[i].Mxx = in[i].Mxx; + out[i].Mxy = in[i].Mxy; + out[i].Myy = in[i].Myy; + out[i].dXccd = in[i].dXccd; + out[i].dYccd = in[i].dYccd; + out[i].dRsys = in[i].dRsys; + out[i].posangle = in[i].posangle; + out[i].pltscale = in[i].pltscale; + out[i].photcode = in[i].photcode; + out[i].dbFlags = in[i].dbFlags; + out[i].photFlags = in[i].photFlags; + } + return (out); +} + +// 'primary' is needed to conform with the API for Loneos and Elixir, but is not used +Average *Average_PS1_V4_ToInternal (Average_PS1_V4 *in, off_t Nvalues, SecFilt **primary) { + + off_t i; + Average *out; + + ALLOCATE_ZERO (out, Average, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].uR = in[i].uR; + out[i].uD = in[i].uD; + out[i].duR = in[i].duR; + out[i].duD = in[i].duD; + out[i].P = in[i].P; + out[i].dP = in[i].dP; + out[i].Xp = in[i].Xp; + out[i].ChiSqAve = in[i].ChiSqAve; + out[i].ChiSqPM = in[i].ChiSqPM; + out[i].ChiSqPar = in[i].ChiSqPar; + out[i].Tmean = in[i].Tmean; + out[i].Trange = in[i].Trange; + out[i].Npos = in[i].Npos; + out[i].Nmeasure = in[i].Nmeasure; + out[i].Nmissing = in[i].Nmissing; + out[i].Nextend = in[i].Nextend; + out[i].measureOffset = in[i].measureOffset; + out[i].missingOffset = in[i].missingOffset; + out[i].extendOffset = in[i].extendOffset; + out[i].flags = in[i].flags; + out[i].photFlagsUpper = in[i].photFlagsUpper; + out[i].photFlagsLower = in[i].photFlagsLower; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].extID = in[i].extID; + } + return (out); +} + +// 'primary' is needed to conform with the API for Loneos and Elixir, but is not used +Average_PS1_V4 *AverageInternalTo_PS1_V4 (Average *in, off_t Nvalues, SecFilt *primary) { + + off_t i; + Average_PS1_V4 *out; + + ALLOCATE_ZERO (out, Average_PS1_V4, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].uR = in[i].uR; + out[i].uD = in[i].uD; + out[i].duR = in[i].duR; + out[i].duD = in[i].duD; + out[i].P = in[i].P; + out[i].dP = in[i].dP; + out[i].Xp = in[i].Xp; + out[i].ChiSqAve = in[i].ChiSqAve; + out[i].ChiSqPM = in[i].ChiSqPM; + out[i].ChiSqPar = in[i].ChiSqPar; + out[i].Tmean = in[i].Tmean; + out[i].Trange = in[i].Trange; + out[i].Npos = in[i].Npos; + out[i].Nmeasure = in[i].Nmeasure; + out[i].Nmissing = in[i].Nmissing; + out[i].Nextend = in[i].Nextend; + out[i].measureOffset = in[i].measureOffset; + out[i].missingOffset = in[i].missingOffset; + out[i].extendOffset = in[i].extendOffset; + out[i].flags = in[i].flags; + out[i].photFlagsUpper = in[i].photFlagsUpper; + out[i].photFlagsLower = in[i].photFlagsLower; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].extID = in[i].extID; + } + return (out); +} + +SecFilt *SecFilt_PS1_V4_ToInternal (SecFilt_PS1_V4 *in, off_t Nvalues) { + + off_t i; + SecFilt *out; + + ALLOCATE_ZERO (out, SecFilt, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M = in[i].M; + out[i].Map = in[i].Map; + out[i].Mkron = in[i].Mkron; + out[i].dMkron = in[i].dMkron; + out[i].dM = in[i].dM; + out[i].Xm = in[i].Xm; + out[i].FluxPSF = in[i].FluxPSF; + out[i].dFluxPSF = in[i].dFluxPSF; + out[i].FluxKron = in[i].FluxKron; + out[i].dFluxKron = in[i].dFluxKron; + out[i].flags = in[i].flags; + out[i].Ncode = in[i].Ncode; + out[i].Nused = in[i].Nused; + out[i].M_20 = in[i].M_20; + out[i].M_80 = in[i].M_80; + out[i].Mstdev = in[i].Mstdev; + out[i].ubercalDist = in[i].ubercalDist; + out[i].stackID = in[i].stackID; + } + return (out); +} + +SecFilt_PS1_V4 *SecFiltInternalTo_PS1_V4 (SecFilt *in, off_t Nvalues) { + + off_t i; + SecFilt_PS1_V4 *out; + + ALLOCATE_ZERO (out, SecFilt_PS1_V4, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].M = in[i].M; + out[i].Map = in[i].Map; + out[i].Mkron = in[i].Mkron; + out[i].dMkron = in[i].dMkron; + out[i].dM = in[i].dM; + out[i].Xm = in[i].Xm; + out[i].FluxPSF = in[i].FluxPSF; + out[i].dFluxPSF = in[i].dFluxPSF; + out[i].FluxKron = in[i].FluxKron; + out[i].dFluxKron = in[i].dFluxKron; + out[i].flags = in[i].flags; + out[i].Ncode = in[i].Ncode; + out[i].Nused = in[i].Nused; + out[i].M_20 = in[i].M_20; + out[i].M_80 = in[i].M_80; + out[i].Mstdev = in[i].Mstdev; + out[i].ubercalDist = in[i].ubercalDist; + out[i].stackID = in[i].stackID; + } + return (out); +} + +Image *Image_PS1_V4_ToInternal (Image_PS1_V4 *in, off_t Nvalues) { + + off_t i; + Image *out; + + ALLOCATE_ZERO (out, Image, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + + strncpy (out[i].name, in[i].name, 120); // out[121], in[121] + out[i].name[120] = 0; // force termination + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit; + out[i].dapmifit = in[i].dapmifit; + out[i].Mcal = in[i].Mcal; + out[i].dMcal = in[i].dMcal; + out[i].Xm = in[i].Xm; + out[i].photcode = in[i].photcode; + out[i].exptime = in[i].exptime; + out[i].sidtime = in[i].sidtime; + out[i].latitude = in[i].latitude; + + out[i].RAo = in[i].RAo; + out[i].DECo = in[i].DECo; + out[i].Radius = in[i].Radius; + + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].ccdnum = in[i].ccdnum; + out[i].flags = in[i].flags; + out[i].imageID = in[i].imageID; + out[i].parentID = in[i].parentID; + out[i].externID = in[i].externID; + out[i].sourceID = in[i].sourceID; + + // as of 2011.02.03, the old Mx,My,..., Mxxxx,Myyyy have been deprecated and replaced + // with the following. (no real databases used those values -- see + // libdvo/doc/dvo-images.txt) + out[i].nLinkAstrom = in[i].nLinkAstrom; + out[i].nLinkPhotom = in[i].nLinkPhotom; + out[i].ubercalDist = in[i].ubercalDist; + out[i].dXpixSys = in[i].dXpixSys; + out[i].dYpixSys = in[i].dYpixSys; + out[i].dMagSys = in[i].dMagSys; + out[i].nFitAstrom = in[i].nFitAstrom; + out[i].nFitPhotom = in[i].nFitPhotom; + out[i].photom_map_id = in[i].photom_map_id; + out[i].astrom_map_id = in[i].astrom_map_id; + } + return (out); +} + +Image_PS1_V4 *ImageInternalTo_PS1_V4 (Image *in, off_t Nvalues) { + + off_t i; + Image_PS1_V4 *out; + + ALLOCATE_ZERO (out, Image_PS1_V4, Nvalues); + + for (i = 0; i < Nvalues; i++) { + memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + + strncpy (out[i].name, in[i].name, 120); // out[121], in[121] + out[i].name[120] = 0; // force termination + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit; + out[i].dapmifit = in[i].dapmifit; + out[i].Mcal = in[i].Mcal; + out[i].dMcal = in[i].dMcal; + out[i].Xm = in[i].Xm; + out[i].photcode = in[i].photcode; + out[i].exptime = in[i].exptime; + out[i].sidtime = in[i].sidtime; + out[i].latitude = in[i].latitude; + + out[i].RAo = in[i].RAo; + out[i].DECo = in[i].DECo; + out[i].Radius = in[i].Radius; + + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].ccdnum = in[i].ccdnum; + out[i].flags = in[i].flags; + out[i].imageID = in[i].imageID; + out[i].parentID = in[i].parentID; + out[i].externID = in[i].externID; + out[i].sourceID = in[i].sourceID; + + // as of 2011.02.03, the old Mx,My,..., Mxxxx,Myyyy have been deprecated and replaced + // with the following. (no real databases used those values -- see + // libdvo/doc/dvo-images.txt) + out[i].nLinkAstrom = in[i].nLinkAstrom; + out[i].nLinkPhotom = in[i].nLinkPhotom; + out[i].ubercalDist = in[i].ubercalDist; + out[i].dXpixSys = in[i].dXpixSys; + out[i].dYpixSys = in[i].dYpixSys; + out[i].dMagSys = in[i].dMagSys; + out[i].nFitAstrom = in[i].nFitAstrom; + out[i].nFitPhotom = in[i].nFitPhotom; + out[i].photom_map_id = in[i].photom_map_id; + out[i].astrom_map_id = in[i].astrom_map_id; + } + return (out); +} + +PhotCode *PhotCode_PS1_V4_To_Internal (PhotCode_PS1_V4 *in, off_t Nvalues) { + + off_t i; + PhotCode *out; + + ALLOCATE_ZERO (out, PhotCode, Nvalues); + + for (i = 0; i < Nvalues; i++) { + strncpy (out[i].name, in[i].name, 31); // out[32], in[32] + out[i].name[31] = 0; // force termination + + out[i].code = in[i].code; + out[i].type = in[i].type; + out[i].C = in[i].C; + out[i].dC = in[i].dC; + out[i].dX = in[i].dX; + out[i].K = in[i].K; + out[i].c1 = in[i].c1; + out[i].c2 = in[i].c2; + out[i].equiv = in[i].equiv; + out[i].Nc = in[i].Nc; + memcpy (out[i].X, in[i].X, 4*sizeof(float)); + + out[i].astromErrSys = in[i].astromErrSys; + out[i].astromErrScale = in[i].astromErrScale; + out[i].astromErrMagScale = in[i].astromErrMagScale; + out[i].photomErrSys = in[i].photomErrSys; + + out[i].photomPoorMask = in[i].photomPoorMask; + out[i].photomBadMask = in[i].photomBadMask; + out[i].astromPoorMask = in[i].astromPoorMask; + out[i].astromBadMask = in[i].astromBadMask; + } + return (out); +} + +PhotCode_PS1_V4 *PhotCode_Internal_To_PS1_V4 (PhotCode *in, off_t Nvalues) { + + off_t i; + PhotCode_PS1_V4 *out; + + ALLOCATE_ZERO (out, PhotCode_PS1_V4, Nvalues); + + for (i = 0; i < Nvalues; i++) { + strncpy (out[i].name, in[i].name, 31); // out[32], in[32] + out[i].name[31] = 0; // force termination + + out[i].code = in[i].code; + out[i].type = in[i].type; + out[i].C = in[i].C; + out[i].dC = in[i].dC; + out[i].dX = in[i].dX; + out[i].K = in[i].K; + out[i].c1 = in[i].c1; + out[i].c2 = in[i].c2; + out[i].equiv = in[i].equiv; + out[i].Nc = in[i].Nc; + memcpy (out[i].X, in[i].X, 4*sizeof(float)); + + out[i].astromErrSys = in[i].astromErrSys; + out[i].astromErrScale = in[i].astromErrScale; + out[i].astromErrMagScale = in[i].astromErrMagScale; + out[i].photomErrSys = in[i].photomErrSys; + + out[i].photomPoorMask = in[i].photomPoorMask; + out[i].photomBadMask = in[i].photomBadMask; + out[i].astromPoorMask = in[i].astromPoorMask; + out[i].astromBadMask = in[i].astromBadMask; + } + return (out); +} Index: trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 34260) @@ -64,4 +64,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -224,4 +232,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 34260) @@ -66,4 +66,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -217,4 +225,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c (revision 34260) @@ -61,4 +61,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -224,4 +232,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c (revision 34260) @@ -61,4 +61,12 @@ out[i].psfNdof = 0; out[i].psfNpix = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; } return (out); @@ -224,4 +232,13 @@ out[i].Mstdev = 0; out[i].ubercalDist = 0; + + // added for PS1_V4 + out[i].Mkron = NAN; + out[i].dMkron = NAN; + out[i].FluxPSF = NAN; + out[i].dFluxPSF = NAN; + out[i].FluxKron = NAN; + out[i].dFluxKron = NAN; + out[i].stackID = 0; } return (out); Index: trunk/Ohana/src/libdvo/src/dvo_image.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_image.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_image.c (revision 34260) @@ -212,4 +212,5 @@ if (db[0].format == DVO_FORMAT_PS1_V2) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_V2"); if (db[0].format == DVO_FORMAT_PS1_V3) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_V3"); + if (db[0].format == DVO_FORMAT_PS1_V4) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_V4"); if (db[0].format == DVO_FORMAT_PS1_REF) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_REF"); Index: trunk/Ohana/src/libdvo/src/dvo_image_raw.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_image_raw.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_image_raw.c (revision 34260) @@ -58,4 +58,5 @@ if (db[0].format == DVO_FORMAT_PS1_V2) ImageSize = sizeof(Image_PS1_V2); if (db[0].format == DVO_FORMAT_PS1_V3) ImageSize = sizeof(Image_PS1_V3); + if (db[0].format == DVO_FORMAT_PS1_V4) ImageSize = sizeof(Image_PS1_V4); if (db[0].format == DVO_FORMAT_PS1_REF) ImageSize = sizeof(Image_PS1_REF); @@ -87,4 +88,5 @@ if (db[0].format == DVO_FORMAT_PS1_V2) gfits_table_mkheader_Image_PS1_V2 (&db[0].theader); if (db[0].format == DVO_FORMAT_PS1_V3) gfits_table_mkheader_Image_PS1_V3 (&db[0].theader); + if (db[0].format == DVO_FORMAT_PS1_V4) gfits_table_mkheader_Image_PS1_V4 (&db[0].theader); if (db[0].format == DVO_FORMAT_PS1_REF) gfits_table_mkheader_Image_PS1_REF (&db[0].theader); Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 34259) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 34260) @@ -278,4 +278,57 @@ } +float PhotAperInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Map; + return (Minst); + } + Minst = measure[0].Map - measure[0].dt - ZERO_POINT; + + return (Minst); +} + +float PhotKron (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Mkron; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Mkron - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotKronInst (Measure *measure) { + + int Np; + float Minst; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Minst = measure[0].Mkron; + return (Minst); + } + Minst = measure[0].Mkron - measure[0].dt - ZERO_POINT; + + return (Minst); +} + float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) { @@ -307,4 +360,23 @@ Msys = Mcat + Mcol; return (Msys); +} + +float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt) { + + int Np; + float ZP; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + ZP = 0.0; + return (ZP); + } + code = &photcodes[0].code[Np]; + ZP = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (ZP); } @@ -488,4 +560,64 @@ } +float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Mkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mkron = (Ns == -1) ? NAN : secfilt[Ns].Mkron; + return (Mkron); +} + +float PhotAveFluxPSF (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Fpsf; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Fpsf = (Ns == -1) ? NAN : secfilt[Ns].FluxPSF; + return (Fpsf); +} + +float PhotAvedFluxPSF (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dFpsf; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dFpsf = (Ns == -1) ? NAN : secfilt[Ns].dFluxPSF; + return (dFpsf); +} + +float PhotAveFluxKron (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float Fkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Fkron = (Ns == -1) ? NAN : secfilt[Ns].FluxKron; + return (Fkron); +} + +float PhotAvedFluxKron (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + float dFkron; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dFkron = (Ns == -1) ? NAN : secfilt[Ns].dFluxKron; + return (dFkron); +} + float PhotM20 (PhotCode *code, Average *average, SecFilt *secfilt) { @@ -522,4 +654,16 @@ Muc = (Ns == -1) ? NAN : secfilt[Ns].ubercalDist; return (Muc); +} + +unsigned int PhotStackID (PhotCode *code, Average *average, SecFilt *secfilt) { + + int Ns; + unsigned int ID; + + if (code == NULL) return 0; + + Ns = photcodes[0].hashNsec[code[0].code]; + ID = (Ns == -1) ? 0 : secfilt[Ns].stackID; + return (ID); } Index: trunk/Ohana/src/libohana/include/ohana.h =================================================================== --- trunk/Ohana/src/libohana/include/ohana.h (revision 34259) +++ trunk/Ohana/src/libohana/include/ohana.h (revision 34260) @@ -155,25 +155,29 @@ /* if your build crashes on OFF_T_MODE, you probably need to add your 64bit hardware to this list */ # ifdef _LARGEFILE_SOURCE -# define OFF_T_FMT "%jd" -# endif +# define OFF_T_FMT "%jd" +# endif + # ifdef lin64 -# define OFF_T_FMT "%jd" -# endif +# define OFF_T_FMT "%jd" +# endif + // mac is annoying # ifdef _DARWIN_C_SOURCE -# define OFF_T_FMT "%lld" +# define OFF_T_FMT "%lld" # endif # ifdef darwin_x86 -# define OFF_T_FMT "%lld" +# define OFF_T_FMT "%lld" # endif # ifdef darwin -# define OFF_T_FMT "%lld" -# endif +# define OFF_T_FMT "%lld" +# endif + // for mac os x 10.7 (lion): they don't use darwin, and i want to make sure it grabs for 10.7 and 64 bit (are they all 64 bit?) # if defined(__APPLE__) && defined(x86_64) -# define OFF_T_FMT "%lld" -# endif +# define OFF_T_FMT "%lld" +# endif + # ifndef OFF_T_FMT -# define OFF_T_FMT "%ld" +# define OFF_T_FMT "%ld" # endif Index: trunk/Ohana/src/opihi/dvo/dvo_host_utils.c =================================================================== --- trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 34259) +++ trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 34260) @@ -15,4 +15,23 @@ snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + // write the dvo comment (host independent) + char commandBase[DVO_MAX_PATH]; + snprintf (commandBase, DVO_MAX_PATH, "dvo.command.%s.txt", uniquer); + char *commandFile = abspath(commandBase, DVO_MAX_PATH); + + FILE *f = fopen (commandFile, "w"); + fprintf (f, "%s\n", basecmd); + if (fflush (f)) DIE("flush", "failed to flush"); + + int fd = fileno (f); + if (fsync (fd)) DIE("fsync", "failed to fsync"); + + if (fclose (f)) DIE("close", "failed to close"); + + // force NFS to write the file to disk + int state; + f = fsetlockfile (commandFile, 0.5, LCK_XCLD, &state); + fclearlockfile (commandFile, f, LCK_XCLD, &state); + int top_status = TRUE; int i; @@ -29,25 +48,6 @@ snprintf (table->hosts[i].results, DVO_MAX_PATH, "%s/dvo.results.%s.fits", table->hosts[i].pathname, uniquer); - char commandBase[DVO_MAX_PATH]; - snprintf (commandBase, DVO_MAX_PATH, "dvo.command.%s.txt", uniquer); - char *commandFile = abspath(commandBase, DVO_MAX_PATH); - - FILE *f = fopen (commandFile, "w"); - fprintf (f, "%s\n", basecmd); - if (fflush (f)) DIE("flush", "failed to flush"); - - int fd = fileno (f); - if (fsync (fd)) DIE("fsync", "failed to fsync"); - - if (fclose (f)) DIE("close", "failed to close"); - - // force NFS to write the file to disk - int state; - f = fsetlockfile (commandFile, 0.5, LCK_XCLD, &state); - fclearlockfile (commandFile, f, LCK_XCLD, &state); - char command[1024]; snprintf (command, 1024, "dvo_client %s -result %s %s -hostID %d -hostdir %s", commandFile, table->hosts[i].results, options, table->hosts[i].hostID, table->hosts[i].pathname); - free (commandFile); if (VERBOSE) gprint (GP_ERR, "command: %s\n", command); @@ -75,4 +75,5 @@ } } + free (commandFile); return top_status; } Index: trunk/Ohana/src/opihi/dvo/gimages.c =================================================================== --- trunk/Ohana/src/opihi/dvo/gimages.c (revision 34259) +++ trunk/Ohana/src/opihi/dvo/gimages.c (revision 34260) @@ -216,10 +216,10 @@ if (PixelCoords) { - gprint (GP_LOG, "%3d %s %6.1f %6.1f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x\n", - Nfound, image[i].name, X, Y, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags); + gprint (GP_LOG, "%3d %s %6.1f %6.1f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x %7d\n", + Nfound, image[i].name, X, Y, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID); } else { XY_to_RD (&ra, &dec, 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); - gprint (GP_LOG, "%3d %s %8.4f %8.4f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x\n", - Nfound, image[i].name, ra, dec, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags); + gprint (GP_LOG, "%3d %s %8.4f %8.4f %20s %5d %2d %4.2f %6.3f %5.3f %5.3f %4x %7d\n", + Nfound, image[i].name, ra, dec, date, image[i].nstar, image[i].photcode, image[i].secz, image[i].Mcal, image[i].dMcal, image[i].exptime, image[i].flags, image[i].imageID); } sprintf (name, "IMAGEx:%d", Nfound); Index: trunk/Ohana/src/opihi/dvo/hosts.c =================================================================== --- trunk/Ohana/src/opihi/dvo/hosts.c (revision 34259) +++ trunk/Ohana/src/opihi/dvo/hosts.c (revision 34260) @@ -2,4 +2,6 @@ # include # define DVO_MAX_PATH 1024 + +enum {TEMP_NONE, TEMP_DVO_RESULTS, TEMP_DVO_LOG, TEMP_RELPHOT_RESULTS, TEMP_RELPHOT_LOG}; // functions to manage the remote hosts @@ -10,4 +12,6 @@ gprint (GP_ERR, " commands:\n"); gprint (GP_ERR, " purge-temp : delete all tempfiles for this shell\n"); + gprint (GP_ERR, " : [-old-pid] [-all-pid] [-v] [-verbose] [-commit] [-type type]\n"); + gprint (GP_ERR, " get-results : determine name of RESULTS file\n"); return FALSE; } @@ -27,5 +31,4 @@ remove_argument (N, &argc, argv); ALL_PID = TRUE; - remove_argument (N, &argc, argv); } @@ -46,4 +49,25 @@ } + // we have the following types of temp files: + // dvo.results.*.fits + // dvo.results.*.fits.log + // log.rlpc.* (relphot logs) + // relphot.catalog.subset.dat [no glob needed] + int TEMP_TYPE = TEMP_DVO_RESULTS; + if ((N = get_argument (argc, argv, "-type"))) { + TEMP_TYPE = TEMP_NONE; + remove_argument (N, &argc, argv); + if (!strcasecmp(argv[N], "dvo.results")) TEMP_TYPE = TEMP_DVO_RESULTS; + if (!strcasecmp(argv[N], "dvo.log")) TEMP_TYPE = TEMP_DVO_LOG; + if (!strcasecmp(argv[N], "relphot.results")) TEMP_TYPE = TEMP_RELPHOT_RESULTS; + if (!strcasecmp(argv[N], "relphot.log")) TEMP_TYPE = TEMP_RELPHOT_LOG; + if (TEMP_TYPE == TEMP_NONE) { + gprint (GP_ERR, "USAGE: hosts purge-temp [-type (type)]\n"); + gprint (GP_ERR, " allowed types dvo.results, dvo.log, relphot.results, relphot.log]\n"); + return FALSE; + } + remove_argument (N, &argc, argv); + } + // XXX wrap this up in a function: char *CATDIR = GetCATDIR(); @@ -63,4 +87,9 @@ } + if (argc != 2) { + gprint (GP_ERR, "USAGE: hosts purge-temp [-all-pid] [-old-pid PID] [-verbose] [-commit]\n"); + return FALSE; + } + int i; for (i = 0; i < table->Nhosts; i++) { @@ -68,7 +97,13 @@ char name[DVO_MAX_PATH]; if (ALL_PID) { - snprintf (name, DVO_MAX_PATH, "%s/dvo.results.*.*.fits", table->hosts[i].pathname); + if (TEMP_TYPE == TEMP_DVO_RESULTS) snprintf (name, DVO_MAX_PATH, "%s/dvo.results.*.fits", table->hosts[i].pathname); + if (TEMP_TYPE == TEMP_DVO_LOG) snprintf (name, DVO_MAX_PATH, "%s/dvo.results.*.fits.log", table->hosts[i].pathname); + if (TEMP_TYPE == TEMP_RELPHOT_RESULTS) snprintf (name, DVO_MAX_PATH, "%s/log.rlpc.*", table->hosts[i].pathname); + if (TEMP_TYPE == TEMP_RELPHOT_LOG) snprintf (name, DVO_MAX_PATH, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname); } else { - snprintf (name, DVO_MAX_PATH, "%s/dvo.results.%05d.*.fits", table->hosts[i].pathname, PID); + if (TEMP_TYPE == TEMP_DVO_RESULTS) snprintf (name, DVO_MAX_PATH, "%s/dvo.results.%05d.*.fits", table->hosts[i].pathname, PID); + if (TEMP_TYPE == TEMP_DVO_LOG) snprintf (name, DVO_MAX_PATH, "%s/dvo.results.%05d.*.fits.log", table->hosts[i].pathname, PID); + if (TEMP_TYPE == TEMP_RELPHOT_RESULTS) snprintf (name, DVO_MAX_PATH, "%s/log.rlpc.*", table->hosts[i].pathname); + if (TEMP_TYPE == TEMP_RELPHOT_LOG) snprintf (name, DVO_MAX_PATH, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname); } if (VERBOSE) gprint (GP_ERR, "checking %s\n", name); @@ -108,6 +143,4 @@ } - - gprint (GP_ERR, "error: unknown hosts command %s\n", argv[1]); return FALSE; Index: trunk/Ohana/src/opihi/dvo/mextract.c =================================================================== --- trunk/Ohana/src/opihi/dvo/mextract.c (revision 34259) +++ trunk/Ohana/src/opihi/dvo/mextract.c (revision 34260) @@ -129,8 +129,9 @@ if (fields[i].ID == MEAS_YCCD) loadImages = TRUE; } - if (fields[i].ID == MEAS_XMOSAIC) loadImages = mosaicMode = TRUE; - if (fields[i].ID == MEAS_YMOSAIC) loadImages = mosaicMode = TRUE; - if (fields[i].ID == MEAS_EXTERN_ID) loadImages = mosaicMode = TRUE; - if (fields[i].ID == MEAS_FLAT) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_XMOSAIC) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_YMOSAIC) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_EXTERN_ID) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_FLAT) loadImages = mosaicMode = TRUE; + if (fields[i].ID == MEAS_CENTER_OFFSET) loadImages = mosaicMode = TRUE; } if (loadImages && !SetImageSelection (mosaicMode, selection)) goto escape; Index: trunk/Ohana/src/opihi/lib.shell/check_stack.c =================================================================== --- trunk/Ohana/src/opihi/lib.shell/check_stack.c (revision 34259) +++ trunk/Ohana/src/opihi/lib.shell/check_stack.c (revision 34260) @@ -1,3 +1,6 @@ # include "opihi.h" +# ifndef MAX_INT +# define MAX_INT 2147483647 +# endif int check_stack (StackVar *stack, int Nstack, int validsize) { @@ -12,7 +15,19 @@ /** if this is a number, put it on the list of scalars and move on. assume value is - * an int unless proven otherwise **/ + * an int unless proven otherwise. + + If we have built libdvo with opihi_int defined as a 32bit signed int, then we have + an overflow for values > MAX_INT (2^31). If the float value is larger than this, + we should treat the value as a float. (NOTE: this means we cannot use 32bit flags, + only 31bit flags. + + **/ + stack[i].FltValue = strtod (stack[i].name, &c1); stack[i].IntValue = strtol (stack[i].name, &c2, 0); + if ((fabs(stack[i].FltValue) > MAX_INT) && (c1 == stack[i].name + strlen (stack[i].name))) { + stack[i].type = 'S'; // 'S' == (float) + continue; + } if (c2 == stack[i].name + strlen (stack[i].name)) { stack[i].type = 's'; // 's' == (int) Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 34259) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 34260) @@ -86,4 +86,10 @@ } + struct timeval now; + gettimeofday (&now, (void *) NULL); + char *moddate = ohana_sec_to_date (now.tv_sec); + gfits_modify (&catalog.header, "RELASTRO", "%s", 1, moddate); + free (moddate); + save_catalogs (&catalog, 1); } Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 34259) +++ trunk/Ohana/src/relphot/Makefile (revision 34260) @@ -49,4 +49,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o @@ -77,4 +78,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 34259) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 34260) @@ -111,4 +111,6 @@ char *BCATALOG; ModeType MODE; + +char *BOUNDARY_TREE; double MAG_LIM; @@ -242,4 +244,5 @@ float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); short getUbercalDist PROTO((off_t meas, int cat)); +float getCenterOffset PROTO((off_t meas, int cat, Measure *measure, unsigned int *myID)); Image *getimage PROTO((off_t N)); Image *getimages PROTO((off_t *N, off_t **LineNumber)); @@ -340,2 +343,6 @@ int client_logger_message (char *format,...); +int MatchImageName (off_t meas, int cat, char *name); + +int load_tree (char *treefile); +int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec); Index: trunk/Ohana/src/relphot/src/BoundaryTreeOps.c =================================================================== --- trunk/Ohana/src/relphot/src/BoundaryTreeOps.c (revision 34260) +++ trunk/Ohana/src/relphot/src/BoundaryTreeOps.c (revision 34260) @@ -0,0 +1,36 @@ +# include "relphot.h" + +// XXX for the moment, only load one boundary tree at a time +// XXX in fact, only allow RINGS.V3... + +static BoundaryTree *tree = NULL; + +int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec) { + + int zone, band; + + if (!primaryCellName) return FALSE; + + primaryCellName[0] = 0; + + if (!tree) return FALSE; + + if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) { + fprintf (stderr, "mismatch!\n"); + return FALSE; + } + + snprintf (primaryCellName, DVO_MAX_PATH, "RINGS.V3.%s", tree->name[zone][band]); + return TRUE; +} + +int load_tree (char *treefile) { + + tree = BoundaryTreeLoad (treefile); + if (!tree) { + fprintf (stderr, "failed to load boundary tree %s\n", treefile); + exit (2); + } + + return TRUE; +} Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 34259) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 34260) @@ -365,4 +365,45 @@ distance = image[i].ubercalDist; // was dummy3 in structure return (distance); +} + +// returns image.Mcal - ff(x,y) +float getCenterOffset (off_t meas, int cat, Measure *measure, unsigned int *myID) { + + off_t i; + float distance; + + if (!MeasureToImage) return -1; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (1000); + + float Xcenter = 0.5*image[i].NX; + float Ycenter = 0.5*image[i].NY; + + *myID = image[i].imageID; + + distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter); + return (distance); +} + +// returns image.Mcal - ff(x,y) +int MatchImageName (off_t meas, int cat, char *name) { + + off_t i; + + if (!name) return FALSE; + if (!name[0]) return FALSE; + + if (!MeasureToImage) return FALSE; + + i = MeasureToImage[cat][meas]; + if (i == -1) return FALSE; + + // this is a bit crude: stack image names are of the form + // RINGS.V3.skycell.1495.027.sky.191211.stk.988232.cmf + // the primaryCell has a name of the form RINGS.V3.skycell.1495 + + if (!strncmp(image[i].name, name, strlen(name))) return TRUE; + return FALSE; } Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 34259) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 34260) @@ -15,5 +15,6 @@ double *wlist; double *aplist; - double *daplist; + double *kronlist; + double *dkronlist; } SetMrelInfo; @@ -164,5 +165,6 @@ SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist ALLOCATE (results.aplist, double, Nmax); - ALLOCATE (results.daplist, double, Nmax); + ALLOCATE (results.kronlist, double, Nmax); + ALLOCATE (results.dkronlist, double, Nmax); for (i = 0; i < Ncatalog; i++) { @@ -174,5 +176,6 @@ SetMrelInfoFree (&results); free (results.aplist); - free (results.daplist); + free (results.kronlist); + free (results.dkronlist); return (TRUE); } @@ -303,21 +306,26 @@ int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { - off_t j, k, m; + off_t j, k, m, ID; int N; float Msys, Mcal, Mmos, Mgrid; - StatType stats, apstats; + StatType stats, apstats, kronstats; liststats_setmode (&stats, STATMODE); liststats_setmode (&apstats, STATMODE); - - double *list = results->list; - double *dlist = results->dlist; - double *wlist = results->wlist; - double *aplist = results->aplist; - double *daplist = results->daplist; + liststats_setmode (&kronstats, STATMODE); + + double *list = results->list; + double *dlist = results->dlist; + double *wlist = results->wlist; + double *aplist = results->aplist; + double *kronlist = results->kronlist; + double *dkronlist = results->dkronlist; SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays int isSetMrelFinal = (pass >= 0); + + char *primaryCell = NULL; + ALLOCATE (primaryCell, char, DVO_MAX_PATH); for (j = 0; j < catalog[Nc].Naverage; j++) { @@ -325,8 +333,10 @@ // option for a test print - if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) { + if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) { fprintf (stderr, "test obj\n"); print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure); } + + BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D); int GoodPS1 = FALSE; @@ -355,4 +365,15 @@ int haveSynth = FALSE; int haveStack = FALSE; + + // need to find the measurement closest to the center of its skycell, as well as the + // closest for the subset of primary projection cells + + float stackCenterOffsetMin = 1e9; + int stackCenterIDmin = -1; + off_t stackCenterMeasureMin = -1; + + float stackPrimaryOffsetMin = 1e9; + int stackPrimaryIDmin = -1; + off_t stackPrimaryMeasureMin = -1; int forceSynth = FALSE; @@ -401,4 +422,8 @@ float Map = PhotAper (&catalog[Nc].measure[m]); aplist[N] = Map - Mcal - Mmos - Mgrid; + + float Mkron = PhotKron (&catalog[Nc].measure[m]); + kronlist[N] = Mkron - Mcal - Mmos - Mgrid; + dkronlist[N] = catalog[Nc].measure[m].dMkron; // special options for PS1 data @@ -416,6 +441,27 @@ // gpc1 stack data if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) { - if (pass < 2) continue; + // if (pass < 2) continue; haveStack = TRUE; + + unsigned int stackImageID; + + // which stack image should we use for the mean value? + // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell + if (MatchImageName (m, Nc, primaryCell)) { + float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID); + if (stackPrimaryOffset < stackPrimaryOffsetMin) { + stackPrimaryOffsetMin = stackPrimaryOffset; + stackPrimaryIDmin = stackImageID; + stackPrimaryMeasureMin = m; + } + } + + // get the center distance for the generic case: + float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID); + if (stackCenterOffset < stackCenterOffsetMin) { + stackCenterOffsetMin = stackCenterOffset; + stackCenterIDmin = stackImageID; + stackCenterMeasureMin = m; + } } @@ -549,7 +595,58 @@ // NOTE : use the modified weight for apmags as well as psf mags - liststats (aplist, daplist, wlist, N, &apstats); - + liststats (aplist, dlist, wlist, N, &apstats); catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean; + + liststats (kronlist, dkronlist, wlist, N, &kronstats); + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mkron = kronstats.mean; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dMkron = kronstats.error; + + if (haveStack) { + m = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin; + ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; + + // get the zero point for the selected image + float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]); + + // flux_cgs : erg sec^1 cm^-2 Hz^-1 + // mag_inst : -2.5 log (cts/sec) + // mag_inst : -2.5 log (flux_inst) + // flux_inst = ten(-0.4*mag_inst) + + // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band] + // flux_cgs = ten(-0.4*(mag_AB + 48.6)) + + // flux_AB : ten(-0.4*mag_AB) + + // flux_cgs = ten(-0.4*48.6) * flux_AB + // flux_AB = ten(+0.4*48.6) * flux_cgs + + // flux_Jy : flux_cgs * 10^23 + + // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy + + // mag_AB = mag_inst + ZP + + // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB + + // flux_AB = flux_inst * ten(-0.4*ZP) + + // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy + + // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23) + // flux_inst = flux_Jy * ten(0.4*ZP - 3.56) + // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56) + + // zpFactor to go from instrumental flux to Janskies + float zpFactor = pow(10.0, -0.4*zp + 3.56); + + // need to put in AB mag factor to get to Janskies (or uJy?) + catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF = zpFactor * catalog[Nc].measure[m].FluxPSF; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF = zpFactor * catalog[Nc].measure[m].dFluxPSF; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron = zpFactor * catalog[Nc].measure[m].FluxKron; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron; + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID = ID; + } // NOTE: for 2MASS measurements, Next should be 1, as should N @@ -607,4 +704,5 @@ } } + if (primaryCell) free (primaryCell); return (TRUE); } Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 34259) +++ trunk/Ohana/src/relphot/src/args.c (revision 34260) @@ -197,4 +197,10 @@ } + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + load_tree (argv[N]); + remove_argument (N, &argc, argv); + } + SHOW_PARAMS = FALSE; if ((N = get_argument (argc, argv, "-params"))) { Index: trunk/Ohana/src/tools/src/ftable.c =================================================================== --- trunk/Ohana/src/tools/src/ftable.c (revision 34259) +++ trunk/Ohana/src/tools/src/ftable.c (revision 34260) @@ -260,4 +260,5 @@ fprintf (stderr, " -ncolumn N: print column number N\n"); fprintf (stderr, " -list: print extension names\n"); + fprintf (stderr, " -layout: describe extension\n"); exit (2); } Index: trunk/Ohana/src/tools/src/glockfile.c =================================================================== --- trunk/Ohana/src/tools/src/glockfile.c (revision 34259) +++ trunk/Ohana/src/tools/src/glockfile.c (revision 34260) @@ -5,6 +5,13 @@ char *filename; double timeout; - int holdtime, state, type; + int N, holdtime, state, type; FILE *f; + + timeout = 30.0; + if ((N = get_argument (argc, argv, "-timeout"))) { + remove_argument (N, &argc, argv); + timeout = atof (argv[N]); + remove_argument (N, &argc, argv); + } if (argc != 4) { @@ -14,6 +21,4 @@ filename = argv[1]; - timeout = 30.0; - holdtime = atoi (argv[3]);