Index: /trunk/Ohana/src/addstar/src/sky_tessalation.c
===================================================================
--- /trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 12750)
+++ /trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 12751)
@@ -32,47 +32,58 @@
 int sky_triangle_to_image (Image *image, SkyTriangle *triangle) {
 
-  int i, peak, b1, b2, NX, NY;
+  int i, n, parity, peak, b1, b2, NX, NY;
   double xv[3], yv[3];	      // coordinates of the vertex in the reference projection 
-  double xo, yo, angle, scale;
+  double xo, yo, xc, yc, 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;
+  for (n = 0; n < 3; n++) {
+    // 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
+    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);
+
+    // find the base and height
+    b1 = (peak + 1) % 3;
+    b2 = (peak + 2) % 3;
+
+    // xo, yo is the center of the baseline
+    xo = 0.5*(xv[b2] + xv[b1]);
+    yo = 0.5*(yv[b2] + yv[b1]);
+
+    if (n != 2) {
+      // xc, yc is the true image center
+      xc = 0.5*(xv[peak] + xo);
+      yc = 0.5*(yv[peak] + yo);
+      XY_to_RD (&triangle[0].r, &triangle[0].d, xc, yc, refcoords);
+    }
+  }  
 
   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));
+  angle = atan2(parity*xv[peak], parity*yv[peak]); // note that this is x/y not y/x (and in radians)
 
   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.pc1_2 = +sin(angle);
+  image[0].coords.pc2_1 = -sin(angle);
   image[0].coords.pc2_2 = +cos(angle);
 
@@ -83,5 +94,4 @@
   }
   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;
@@ -95,4 +105,7 @@
   image[0].NY = NY;
   image[0].code = 1;
+
+  // fprintf (stderr, "%d  %f  %d %d %d  %d %d\n", parity, angle*DEG_RAD, peak, b1, b2, NX, NY);
+
   return (TRUE);
 }
@@ -122,4 +135,6 @@
     triangle[0].rv[i] = DEG_RAD * atan2(triangle[0].vertex[i].y, triangle[0].vertex[i].x);
   }
+
+  return TRUE;
 }
 
@@ -128,5 +143,5 @@
 SkyTriangle *sky_divide_triangles (SkyTriangle *in, int *ntriangles) {
 
-  int i, j, Ntriangles;
+  int i, j, Ntriangles, Nt;
   SkyTriangle *out;
 
@@ -134,4 +149,5 @@
   ALLOCATE (out, SkyTriangle, Ntriangles);
 
+  Nt = 0;
   for (i = 0; i < *ntriangles; i++) {
     for (j = 0; j < 3; j++) {
@@ -145,6 +161,19 @@
       out[4*i + 3].vertex[j] = out[4*i + j].vertex[1];
     }
-  }
-  *ntriangles = Ntriangles;
+# if (0)
+    fprintf (stdout, "%f %f %f  %f %f %f  %f %f %f\n", 
+	     in[i].vertex[0].x, in[i].vertex[0].y, in[i].vertex[0].z, 
+	     in[i].vertex[1].x, in[i].vertex[1].y, in[i].vertex[1].z, 
+	     in[i].vertex[2].x, in[i].vertex[2].y, in[i].vertex[2].z);
+    for (j = 0; j < 4; j++) {
+      fprintf (stdout, "%f %f %f  %f %f %f  %f %f %f\n", 
+	       out[4*i+j].vertex[0].x, out[4*i+j].vertex[0].y, out[4*i+j].vertex[0].z, 
+	       out[4*i+j].vertex[1].x, out[4*i+j].vertex[1].y, out[4*i+j].vertex[1].z, 
+	       out[4*i+j].vertex[2].x, out[4*i+j].vertex[2].y, out[4*i+j].vertex[2].z);
+    }
+# endif  
+    Nt += 4;
+  }
+  *ntriangles = Nt;
   return (out);
 }
@@ -239,10 +268,10 @@
     tri[i].vertex[1].x = +ctht*cos(psi);
     tri[i].vertex[1].y = +ctht*sin(psi);
-    tri[i].vertex[1].z = +stht;
+    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[2].z = -stht;
 
     tri[i].vertex[0].x = +0;