Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 12748) +++ /trunk/Ohana/src/addstar/Makefile (revision 12749) @@ -215,8 +215,17 @@ $(SRC)/SetSignals.$(ARCH).o +SKYCELLS = \ +$(SRC)/skycells.$(ARCH).o \ +$(SRC)/args_skycells.$(ARCH).o \ +$(SRC)/ConfigInit_skycells.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/sky_tessalation.$(ARCH).o + $(ADDSTARC) : $(INC)/addstar.h $(ADDSTARD) : $(INC)/addstar.h $(ADDSTART) : $(INC)/addstar.h $(ADDSTAR) : $(INC)/addstar.h +$(SKYCELLS) : $(INC)/addstar.h $(BIN)/addstar.$(ARCH) : $(ADDSTAR) @@ -227,6 +236,7 @@ $(BIN)/sedstar.$(ARCH) : $(SEDSTAR) $(BIN)/load2mass.$(ARCH) : $(LOAD-2MASS) - -INSTALL = addstar addstarc addstard addstart sedstar load2mass +$(BIN)/skycells.$(ARCH) : $(SKYCELLS) + +INSTALL = addstar addstarc addstard addstart sedstar load2mass skycells # dependancy rules for binary code ######################### Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 12748) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 12749) @@ -226,4 +226,7 @@ int NewReflist_Thread PROTO((int BindSocket)); +int args_skycells (int argc, char **argv); +int ConfigInit_skycells (int *argc, char **argv); + // this is a gnu extension?? caution! void *memrchr(const void *s, int c, size_t n); Index: /trunk/Ohana/src/addstar/include/skycells.h =================================================================== --- /trunk/Ohana/src/addstar/include/skycells.h (revision 12749) +++ /trunk/Ohana/src/addstar/include/skycells.h (revision 12749) @@ -0,0 +1,65 @@ +# include +# include +# include +# include +# include +# include + +/* linux is happy with this, not solaris */ +# include +# include +# include +# include + +typedef struct { + double x, y, z; +} Point; + +typedef struct { + Point vertex[3]; // triangle vertices (3d) + Point center; // triangle center (3d) + double r, d; // triangle center (2d) + double rv[3], dv[3]; // triangle center (2d) +} SkyTriangle; + +/* globals which define database info / data sources (KEEP) */ +char ImageCat[256]; +char GSCFILE[256]; +char CATDIR[256]; +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char PASSWORD[80]; +char HOSTNAME[80]; +int NVALID, *VALID_IP; +char SKY_TABLE[256]; +int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ +char CameraLayout[256]; +SkyTable *ServerSky; + +/* used to select entries from header (gstars or parse_time) (KEEP) */ +char DateKeyword[64]; +char DateMode[64]; +char UTKeyword[64]; +char MJDKeyword[64]; +char JDKeyword[64]; +char ExptimeKeyword[64]; +char AirmassKeyword[64]; +char CCDNumKeyword[64]; +char STKeyword[64]; + +int VERBOSE; + +void SetProtect PROTO((int mode)); +int SetSignals PROTO(()); +int Shutdown PROTO((char *message, ...); ) +void TrapSignal PROTO((int sig)); +int args_skycells (int argc, char **argv); +int ConfigInit_skycells (int *argc, char **argv); + +Image *sky_tessalation (int level, int *nimages); +int sky_triangle_to_image (Image *image, SkyTriangle *triangle); +int sky_triangle_coords (SkyTriangle *triangle); +SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles); +Point sky_divide_edge (Point v1, Point v2); +SkyTriangle *sky_base_triangles (int *ntriangles); +int sky_tessalation_init (); Index: /trunk/Ohana/src/addstar/src/ConfigInit_skycells.c =================================================================== --- /trunk/Ohana/src/addstar/src/ConfigInit_skycells.c (revision 12749) +++ /trunk/Ohana/src/addstar/src/ConfigInit_skycells.c (revision 12749) @@ -0,0 +1,78 @@ +# include "addstar.h" + +void GetConfig (char *config, char *field, char *format, int N, void *ptr); + +int ConfigInit_skycells (int *argc, char **argv) { + + char *config, *file; + char RadiusWord[80]; + char CatdirPhotcodeFile[256]; + char MasterPhotcodeFile[256]; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (1); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + /* used by parse_time to find time-related keywords */ + strcpy (DateKeyword, "NONE"); + strcpy (DateMode, "NONE"); + strcpy (UTKeyword, "NONE"); + strcpy (JDKeyword, "NONE"); + strcpy (MJDKeyword, "NONE"); + ScanConfig (config, "DATE-KEYWORD", "%s", 0, DateKeyword); + ScanConfig (config, "DATE-MODE", "%s", 0, DateMode); + ScanConfig (config, "UT-KEYWORD", "%s", 0, UTKeyword); + ScanConfig (config, "MJD-KEYWORD", "%s", 0, MJDKeyword); + ScanConfig (config, "JD-KEYWORD", "%s", 0, JDKeyword); + + ScanConfig (config, "EXPTIME-KEYWORD", "%s", 0, ExptimeKeyword); + ScanConfig (config, "AIRMASS-KEYWORD", "%s", 0, AirmassKeyword); + ScanConfig (config, "CCDNUM-KEYWORD", "%s", 0, CCDNumKeyword); + ScanConfig (config, "ST-KEYWORD", "%s", 0, STKeyword); + + GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + GetConfig (config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); + ScanConfig (config, "CATMODE", "%s", 0, CATMODE); + ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { + SKY_DEPTH = SKY_DEPTH_HST; + } + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { + SKY_TABLE[0] = 0; + } + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + /* default mode, format, if not specified */ + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + /* XXX this does not yet write out the master photcode table */ + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) { + fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile); + exit (1); + } + + free (config); + free (file); + return (TRUE); +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /trunk/Ohana/src/addstar/src/args_skycells.c =================================================================== --- /trunk/Ohana/src/addstar/src/args_skycells.c (revision 12749) +++ /trunk/Ohana/src/addstar/src/args_skycells.c (revision 12749) @@ -0,0 +1,36 @@ +# include "addstar.h" +static void help (void); + +int args_skycells (int argc, char **argv) { + + int i, N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc == 2) return (TRUE); + + fprintf (stderr, "USAGE: skycells (level)\n"); + exit (2); +} + +static void help () { + + fprintf (stderr, "USAGE\n"); + fprintf (stderr, " skycells\n\n"); + + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} Index: /trunk/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 12749) +++ /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 12749) @@ -0,0 +1,270 @@ +# include "skycells.h" +# include "assert.h" + +static Coords *refcoords = NULL; + +Image *sky_tessalation (int level, int *nimages) { + + int i, Ntriangles; + SkyTriangle *tri, *new; + Image *image; + + sky_tessalation_init (); + + tri = sky_base_triangles (&Ntriangles); + + for (i = 0; i < level; i++) { + new = sky_divide_triangles (tri, &Ntriangles); + free (tri); + tri = new; + } + + ALLOCATE (image, Image, Ntriangles); + for (i = 0; i < Ntriangles; i++) { + sky_triangle_to_image (&image[i], &tri[i]); + } + + *nimages = Ntriangles; + return (image); +} + +// an allocated image is supplied, we fill in the values +int sky_triangle_to_image (Image *image, SkyTriangle *triangle) { + + int i, peak, b1, b2, NX, NY; + double xv[3], yv[3]; // coordinates of the vertex in the reference projection + double xo, yo, angle, scale; + + // calculate the triangle coordinates in r,d + sky_triangle_coords (triangle); + + // we will project to the triangle center position + refcoords[0].crval1 = triangle[0].r; + refcoords[0].crval2 = triangle[0].d; + + // find the size, rotation, and parity of the image + // project the vertices and find the image parity + int parity = 1; + for (i = 0; i < 3; i++) { + RD_to_XY (&xv[i], &yv[i], triangle[0].rv[i], triangle[0].dv[i], refcoords); + parity *= SIGN(yv[i]); + } + + // choose the peak vertex + peak = -1; + for (i = 0; (peak == -1) && (i < 3); i++) { + if (parity == SIGN(yv[i])) { + peak = i; + } + } + assert (peak != -1); + angle = atan2(xv[peak], yv[peak]); // note that this is x/y not y/x (and in radians) + + // find the base and height + b1 = (peak + 1) % 3; + b2 = (peak + 2) % 3; + + NX = hypot((xv[b2]-xv[b1]),(yv[b2]-yv[b1])); + xo = 0.5*(xv[b2] + xv[b1]); + yo = 0.5*(yv[b2] + yv[b1]); + NY = hypot((xv[peak]-xo),(yv[peak]-yo)); + + memset (image, 0, sizeof(Image)); + image[0].coords = *refcoords; + image[0].coords.pc1_1 = +cos(angle); + image[0].coords.pc1_2 = -sin(angle); + image[0].coords.pc2_1 = +sin(angle); + image[0].coords.pc2_2 = +cos(angle); + + if (parity == 1) { + strcpy (image[0].coords.ctype, "TRP--TAN"); + } else { + strcpy (image[0].coords.ctype, "TRM--TAN"); + } + sprintf (image[0].name, "test.%03d", i); + fprintf (stderr, "NX, NY: %d, %d\n", NX, NY); + if ((NX > 60000) || (NY > 60000)) { + scale = NX / 60000.0; + scale = MAX(scale, NY / 60000.0); + NX /= scale; + NY /= scale; + image[0].coords.cdelt1 *= scale; + image[0].coords.cdelt2 *= scale; + } + image[0].NX = NX; + image[0].NY = NY; + image[0].code = 1; + return (TRUE); +} + +int sky_triangle_coords (SkyTriangle *triangle) { + + int i; + double r; + + // calculate the triangle center + triangle[0].center.x = (triangle[0].vertex[0].x + triangle[0].vertex[1].x + triangle[0].vertex[2].x)/3.0; + triangle[0].center.y = (triangle[0].vertex[0].y + triangle[0].vertex[1].y + triangle[0].vertex[2].y)/3.0; + triangle[0].center.z = (triangle[0].vertex[0].z + triangle[0].vertex[1].z + triangle[0].vertex[2].z)/3.0; + + // renormalize + r = 1.0 / sqrt (SQ(triangle[0].center.x) + SQ(triangle[0].center.y) + SQ(triangle[0].center.z)); + + triangle[0].center.x *= r; + triangle[0].center.y *= r; + triangle[0].center.z *= r; + + triangle[0].d = DEG_RAD * asin(triangle[0].center.z); + triangle[0].r = DEG_RAD * atan2(triangle[0].center.y, triangle[0].center.x); + + for (i = 0; i < 3; i++) { + triangle[0].dv[i] = DEG_RAD * asin(triangle[0].vertex[i].z); + triangle[0].rv[i] = DEG_RAD * atan2(triangle[0].vertex[i].y, triangle[0].vertex[i].x); + } +} + +// take a list of triangles from one level and return a list of triangles in the next level +// we are doing basic edge division, always yielding 4x as many new triangles as old; +SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles) { + + int i, j, Ntriangles; + SkyTriangle *out; + + Ntriangles = *ntriangles * 4; + ALLOCATE (out, SkyTriangle, Ntriangles); + + for (i = 0; i < *ntriangles; i++) { + for (j = 0; j < 3; j++) { + out[4*i + j].vertex[0] = in[i].vertex[j]; + out[4*i + j].vertex[1] = sky_divide_edge (in[i].vertex[j], in[i].vertex[(j+1)%3]); + } + for (j = 0; j < 3; j++) { + out[4*i + j].vertex[2] = out[4*i + (j+2)%3].vertex[1]; + } + for (j = 0; j < 3; j++) { + out[4*i + 3].vertex[j] = out[4*i + j].vertex[1]; + } + } + *ntriangles = Ntriangles; + return (out); +} + +// take a list of triangles from one level and return a list of triangles in the next level +// we are doing basic edge division, always yielding 4x as many new triangles as old; +Point sky_divide_edge (Point v1, Point v2) { + + double r; + Point out; + + out.x = v1.x + v2.x; + out.y = v1.y + v2.y; + out.z = v1.z + v2.z; + + r = 1.0 / sqrt (SQ(out.x) + SQ(out.y) + SQ(out.z)); + + out.x *= r; + out.y *= r; + out.z *= r; + + return (out); +} + +# define THETA RAD_DEG*26.565 +# define D_PSI RAD_DEG*360.0/5.0 + +SkyTriangle *sky_base_triangles (int *ntriangles) { + + int i; + double ctht, stht, psi; + SkyTriangle *tri; + + // generate 0-level triangles + ALLOCATE (tri, SkyTriangle, 20); + + ctht = cos(THETA); + stht = sin(THETA); + + for (i = 0; i < 5; i++) { + tri[i].vertex[0].x = +0; + tri[i].vertex[0].y = +0; + tri[i].vertex[0].z = +1; + + psi = (i + 0.0)*D_PSI; + tri[i].vertex[1].x = +ctht*cos(psi); + tri[i].vertex[1].y = +ctht*sin(psi); + tri[i].vertex[1].z = +stht; + + psi = (i + 1.0)*D_PSI; + tri[i].vertex[2].x = +ctht*cos(psi); + tri[i].vertex[2].y = +ctht*sin(psi); + tri[i].vertex[2].z = +stht; + } + + for (i = 5; i < 10; i++) { + psi = (i + 0.0)*D_PSI; + tri[i].vertex[0].x = +ctht*cos(psi); + tri[i].vertex[0].y = +ctht*sin(psi); + tri[i].vertex[0].z = +stht; + + psi = (i + 0.5)*D_PSI; + tri[i].vertex[1].x = +ctht*cos(psi); + tri[i].vertex[1].y = +ctht*sin(psi); + tri[i].vertex[1].z = -stht; + + psi = (i + 1.0)*D_PSI; + tri[i].vertex[2].x = +ctht*cos(psi); + tri[i].vertex[2].y = +ctht*sin(psi); + tri[i].vertex[2].z = +stht; + } + + for (i = 10; i < 15; i++) { + psi = (i + 0.5)*D_PSI; + tri[i].vertex[0].x = +ctht*cos(psi); + tri[i].vertex[0].y = +ctht*sin(psi); + tri[i].vertex[0].z = -stht; + + psi = (i + 1.0)*D_PSI; + tri[i].vertex[1].x = +ctht*cos(psi); + tri[i].vertex[1].y = +ctht*sin(psi); + tri[i].vertex[1].z = +stht; + + psi = (i + 1.5)*D_PSI; + tri[i].vertex[2].x = +ctht*cos(psi); + tri[i].vertex[2].y = +ctht*sin(psi); + tri[i].vertex[2].z = -stht; + } + + for (i = 15; i < 20; i++) { + psi = (i + 0.5)*D_PSI; + tri[i].vertex[1].x = +ctht*cos(psi); + tri[i].vertex[1].y = +ctht*sin(psi); + tri[i].vertex[1].z = +stht; + + psi = (i + 1.5)*D_PSI; + tri[i].vertex[2].x = +ctht*cos(psi); + tri[i].vertex[2].y = +ctht*sin(psi); + tri[i].vertex[2].z = +stht; + + tri[i].vertex[0].x = +0; + tri[i].vertex[0].y = +0; + tri[i].vertex[0].z = -1; + } + + *ntriangles = 20; + return tri; +} + +int sky_tessalation_init () { + + ALLOCATE (refcoords, Coords, 1); + refcoords[0].crval1 = refcoords[0].crval2 = 0.0; + refcoords[0].crpix1 = refcoords[0].crpix2 = 0.0; + refcoords[0].cdelt1 = refcoords[0].cdelt2 = 1.0 / 3600; + refcoords[0].pc1_1 = refcoords[0].pc2_2 = 1.0; + refcoords[0].pc1_2 = refcoords[0].pc2_1 = 0.0; + refcoords[0].Npolyterms = 0; + memset (refcoords[0].polyterms, 0, 14*sizeof(float)); + strcpy (refcoords[0].ctype, "RA---TAN"); + return (TRUE); +} + Index: /trunk/Ohana/src/addstar/src/skycells.c =================================================================== --- /trunk/Ohana/src/addstar/src/skycells.c (revision 12749) +++ /trunk/Ohana/src/addstar/src/skycells.c (revision 12749) @@ -0,0 +1,75 @@ +# include "addstar.h" + +int main (int argc, char **argv) { + + int i, Nimages, status, level; + Image *images; + Catalog catalog; + FITS_DB db; + + // XXX for test data + Coords coords; + double r, d; + + SetSignals (); + ConfigInit_skycells (&argc, argv); + args_skycells (argc, argv); + level = atoi (argv[1]); + + images = sky_tessalation (level, &Nimages); + + /*** update the image table ***/ + /* setup image table format and lock */ + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + /* load or create the image table */ + if (db.dbstate == LCK_EMPTY) { + if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); + dvo_image_create (&db, ZeroPt); + } else { + if (!dvo_image_load (&db, VERBOSE, FORCE_READ)) { + Shutdown ("can't read image catalog %s", db.filename); + } + } + + /* add the new images and save */ + dvo_image_addrows (&db, images, Nimages); + SetProtect (TRUE); + dvo_image_update (&db, VERBOSE); + SetProtect (FALSE); + + dvo_image_unlock (&db); + exit (0); +} + +# if (0) + coords.crval1 = coords.crval2 = 0.0; + coords.crpix1 = coords.crpix2 = 0.0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 0; + memset (coords.polyterms, 0, 14*sizeof(float)); + strcpy (coords.ctype, "TRP--TAN"); + + /** quick demo to create a set of triangular images **/ + Nimages = 100; + ALLOCATE (images, Image, Nimages); + for (r = 0, i = 0; r < 10; r++) { + for (d = 0; d < 10; d++, i++) { + memset (&images[i], 0, sizeof(Image)); + images[i].coords = coords; + images[i].coords.crval1 = r; + images[i].coords.crval2 = d; + if (i%2) strcpy (images[i].coords.ctype, "TRM--TAN"); + sprintf (images[i].name, "test.%03d", i); + images[i].NX = 1000; + images[i].NY = 1000; + images[i].code = 1; + } + } +# endif +