Index: /trunk/Ohana/src/addstar/include/skycells.h
===================================================================
--- /trunk/Ohana/src/addstar/include/skycells.h	(revision 13185)
+++ /trunk/Ohana/src/addstar/include/skycells.h	(revision 13186)
@@ -53,4 +53,5 @@
 int    MODE;
 int    NMAX;
+int    NX_SUB, NY_SUB;
 double SCALE;
 
Index: /trunk/Ohana/src/addstar/src/args_skycells.c
===================================================================
--- /trunk/Ohana/src/addstar/src/args_skycells.c	(revision 13185)
+++ /trunk/Ohana/src/addstar/src/args_skycells.c	(revision 13186)
@@ -43,4 +43,17 @@
   }
 
+  /* max number of skycells kept in memory */
+  NX_SUB = NY_SUB = 1;
+  if ((N = get_argument (argc, argv, "-nx"))) {
+    remove_argument (N, &argc, argv);
+    NX_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-ny"))) {
+    remove_argument (N, &argc, argv);
+    NY_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc == 2) {
     strtol (argv[1], &ptr, 10);
@@ -48,5 +61,5 @@
   }
 
-  fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix] [-nmax Max Cells in memory]\n");
+  fprintf (stderr, "USAGE: skycells (level) [-triangles] [-scale arcsec/pix] [-nmax Max Cells in memory] [-nx (Nx cells)] [-ny (Ny cells)]\n");
   exit (2);
 }
Index: /trunk/Ohana/src/addstar/src/sky_tessalation.c
===================================================================
--- /trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 13185)
+++ /trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 13186)
@@ -93,5 +93,5 @@
 int sky_tessalation_squares (FITS_DB *db, int level, int Nmax) {
 
-  int i, j, Ndigit, Ntriangles, Nbase, Ntotal, Ntop, Ltop, Nout;
+  int i, j, Ndigit, Ntriangles, Nbase, Ntotal, Ntop, Ltop, Nout, Nx, Ny;
   double fLtop;
   SkyTriangle *base, *tri, *new;
@@ -99,4 +99,7 @@
   char format[16];
 
+  Nx = NX_SUB;
+  Ny = NY_SUB;
+
   // generate the initial base set
   base = sky_base_triangles (&Nbase);
@@ -108,5 +111,5 @@
 
   // to what depth do we need to go to have only Nmax foreach subcell?
-  Ntop = Ntotal / Nmax;
+  Ntop = Ntotal / (Nmax*Nx*Ny) ;
   if (Ntop > Nbase) {
     fLtop = log10(Ntotal / (double)(Ntop * Nbase)) / log10(4.0);
@@ -154,4 +157,15 @@
     free (image);
 
+    // subdivide each image (Nx x Ny subcells)
+    if (Nx*Ny > 1) {
+      ALLOCATE (image, Image, Nx*Ny*Nout);
+      for (j = 0; j < Nout; j++) {
+	sky_subdivide_image (&image[j*Nx*Ny], &out[j], Nx, Ny);
+      }
+      free (out);
+      out = image;
+      Nout *= Nx*Ny;
+    }
+
     /* add the new images and save */
     dvo_image_addrows (db, out, Nout);
@@ -229,4 +243,32 @@
   image[0].Mxxy = xv[2];  image[0].Myyy = yv[2];
 
+  return (TRUE);
+}
+
+// an allocated image set is supplied, we fill in the values
+int sky_subdivide_image (Image *output, Image *input, int Nx, int Ny) {
+
+  int i, j, N, NX, NY;
+
+  NX = input[0].NX/(double)Nx + 0.5;
+  NY = input[0].NY/(double)Ny + 0.5;
+
+  N = 0;
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++) {
+
+      memcpy (&output[N], &input[0], sizeof(Image));
+
+      // output[0].coords = input[0].coords;
+      // strcpy (output[0].coords.ctype, input[0].coords.ctype);
+
+      output[N].NX = NX;
+      output[N].NY = NY;
+
+      output[N].coords.crpix1 = input[0].coords.crpix1 - i*NX;
+      output[N].coords.crpix2 = input[0].coords.crpix2 - j*NY;
+      N++;
+    }
+  }
   return (TRUE);
 }