Index: /branches/eam_branches/eam_branch_20090303/Ohana/src/addstar/src/sky_tessalation.c
===================================================================
--- /branches/eam_branches/eam_branch_20090303/Ohana/src/addstar/src/sky_tessalation.c	(revision 23219)
+++ /branches/eam_branches/eam_branch_20090303/Ohana/src/addstar/src/sky_tessalation.c	(revision 23220)
@@ -240,11 +240,18 @@
   // integer number of dec zones between -90 and +90
 
-  nDEC = 180.0 / CELLSIZE;
-  dDEC = 180.0 / nDEC;
+  // 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;
+
+  // a test
+  // for (dec = 0.0 + 0.5*dDEC; dec < +90.0; dec += dDEC) {
 
   // generate the a collection of rectangles for each ring
-  for (dec = -90.0 + 0.5*dDEC; dec < +90.0; dec += dDEC) {
+  for (dec = -90.0; dec < +90.0 + 0.5*dDEC; dec += dDEC) {
 
     ring = sky_rectangle_ring (dec, dDEC, &Nring);
+    if (!ring) continue;
 
     // subdivide each image (Nx x Ny subcells)
@@ -529,7 +536,7 @@
 SkyRectangle *sky_rectangle_ring (float dec, float dDEC, int *nring) {
 
-  int i, nRA, NX, NY;
-  float dRA, decLower;
+  int i, NX, NY, nRA;
   SkyRectangle *ring;
+  float theta, dRA;
 
   // 'dec' is a guess at the center of the cell; in fact, we need to choose decLower and
@@ -537,9 +544,62 @@
 
   // we can determine the 'lower' bound (bound closest to the equator):
-  decLower = (dec > 0.0) ? dec - 0.5*dDEC : dec + 0.5*dDEC;
-
-  // Subdivide the 'lower' bound into an integer number of segments:
-  nRA = cos(dec*RAD_DEG) * 360.0 / CELLSIZE; // CELLSIZE is a projection size
-  dRA = 360.0 / nRA;                         // dRA is a size in RA degrees
+  float decLower = (dec > 0.0) ? dec - 0.5*dDEC : dec + 0.5*dDEC;
+
+  // solve for actual cellsize (\theta):  tan(\delta_{n+1} - \theta/2) = tan(\delta_n + \theta/2)cos(\alpha_n / 2)
+  float decUpper = (dec > 0.0) ? dec + dDEC : dec - dDEC;
+
+  if (fabs(dec) + 0.5*dDEC > 90.0) {
+    // onPole = TRUE;
+    theta = dDEC;
+    nRA = 1;
+    dRA = theta / cos(decLower*RAD_DEG); // make a square at the pole
+  } else {
+    // onPole = FALSE;
+    // Subdivide the 'lower' bound into an integer number of segments:
+    nRA = cos(RAD_DEG*decLower) * 360.0 / CELLSIZE; // CELLSIZE is a projection size
+    dRA = 360.0 / nRA;                         // dRA is a size in RA degrees == \alpha_n
+
+    // tan(decUpper - theta/2) = tan(dec + theta/2) cos(dRA / 2);
+
+    // we solve this equation for theta (fairly ugly: expand the tangents into sin/cos, expand the 
+    // sum-of-angle sine and cosine, multiply through, convert via half-angle formulae and write 
+    // as a quadratic expression in sine(theta/2)
+  
+    float sd1 = sin(RAD_DEG*decUpper);
+    float cd1 = cos(RAD_DEG*decUpper);
+    float sd2 = sin(RAD_DEG*dec);
+    float cd2 = cos(RAD_DEG*dec);
+    float   k = cos(RAD_DEG*dRA/2.0);
+
+    float c1 =  (sd1*cd2 + sd2*cd1)*(1.0 - k);
+    float c2 =  (sd1*cd2 - sd2*cd1)*(1.0 + k);
+    float c3 = -(sd1*sd2 + cd1*cd2)*(1.0 + k); 
+
+    float A = SQ(c3) + SQ(c2);
+    float B = 2*c1*c3;
+    float C = SQ(c1) - SQ(c2);
+
+    float arg = SQ(B) - 4.0*A*C;
+
+    float root;
+
+    if (dec >= 0.0) {
+      root = (-B + sqrt (arg)) / (2.0*A);
+      theta = +DEG_RAD*asin(root);
+    } else {
+      root = (-B - sqrt (arg)) / (2.0*A);
+      theta = -DEG_RAD*asin(root);
+    }
+
+    // the negative solution yields a negative cellsize 
+    // float root2 = (-B - sqrt (arg)) / (2.0*A);
+    // float theta2 = DEG_RAD*asin(root2);
+
+    // test lines:
+    // float r1 = tan(RAD_DEG*(decUpper - 0.5*theta1));
+    // float r2 = tan(RAD_DEG*(dec + 0.5*theta1));
+    // fprintf (stdout, "%f %f  %f  %f  %f %f  %f %f  %f %f %f\n", dec, decUpper, dRA, arg, root1, root2, theta1, theta2, r1, r2, k*r2);
+  }
+  fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
 
   // I think we need to return the value of dec for the next ring, but I am not sure...
@@ -559,6 +619,6 @@
   
     // range values are in projected degrees
-    NX = cos(dec*RAD_DEG) * dRA  * 3600.0 / SCALE;
-    NY =                    dDEC * 3600.0 / SCALE;
+    NX = cos(decLower*RAD_DEG) * dRA   * 3600.0 / SCALE;
+    NY =                         theta * 3600.0 / SCALE;
 
     // crpix1,crpix2 is the projection center
@@ -576,7 +636,7 @@
 
 
-    fprintf (stderr, "%f %f  : %f %f\n", 
-	     ring[i].coords.crval1, ring[i].coords.crval2, 
-	     ring[i].coords.crpix1, ring[i].coords.crpix2);
+    // fprintf (stderr, "%f %f  : %f %f\n", 
+    // ring[i].coords.crval1, ring[i].coords.crval2, 
+    // ring[i].coords.crpix1, ring[i].coords.crpix2);
   }