Index: trunk/Ohana/src/addstar/Makefile =================================================================== --- trunk/Ohana/src/addstar/Makefile (revision 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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 34088) +++ 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