Index: trunk/ppSim/src/ppSimLoadStars.c
===================================================================
--- trunk/ppSim/src/ppSimLoadStars.c	(revision 18047)
+++ trunk/ppSim/src/ppSimLoadStars.c	(revision 18464)
@@ -35,6 +35,8 @@
     float y0fpa = 0.5*(bounds->y0 + bounds->y1);
 
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
-    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius);
+    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "",
+                  ra0 - radius / cos(dec0));
+    psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "",
+                  ra0 + radius / cos(dec0));
     psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius);
     psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius);
@@ -50,4 +52,6 @@
     stars = psArrayRealloc (stars, refStars->n);
 
+    psProjection *proj = psProjectionAlloc(ra0, dec0, scale, scale, PS_PROJ_TAN); // Projection
+
     // Conversion loop
     for (long i = 0; i < refStars->n; i++) {
@@ -60,11 +64,13 @@
 
         // Convert to x,y position on tangent plane, in pixels
-        float div = (sin(star->ra) * sin(dec0) + cos(star->dec) * cos(dec0) * cos(star->ra - ra0)) * scale;
-        float xi = cos(star->dec) * sin(star->ra - ra0) / div;
-        float eta = (sin(star->dec) * cos(dec0) - cos(star->dec) * sin(dec0) * cos(star->ra - ra0)) / div;
+        psPlane plane;                  // Plane (detector) coordinates
+        psSphere sphere;                // Sphere (sky) coordinates
+        sphere.r = star->ra;
+        sphere.d = star->dec;
+        psProject(&plane, &sphere, proj);
 
         // Apply rotation, make FPA center of boresite
-        star->x = cos(pa) * xi - sin(pa) * eta + x0fpa;
-        star->y = sin(pa) * xi + cos(pa) * eta + y0fpa;
+        star->x = cos(pa) * plane.x - sin(pa) * plane.y + x0fpa;
+        star->y = sin(pa) * plane.x + cos(pa) * plane.y + y0fpa;
 
         // Convert magnitude to peak flux
@@ -77,4 +83,5 @@
     }
     stars->n = oldSize + refStars->n;
+    psFree(proj);
 
     pmLumFunc *lumfunc = psastroLuminosityFunction (refStars);
@@ -83,14 +90,14 @@
     psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);
     psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0);
-    
+
     if (lumfunc) {
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MIN",   PS_META_REPLACE, "min valid magnitude",      	  lumfunc->mMin);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MAX",   PS_META_REPLACE, "max valid magnitude",      	  lumfunc->mMax);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.SLOPE",  PS_META_REPLACE, "log-mag histogram slope",  	  lumfunc->slope);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.OFFSET", PS_META_REPLACE, "log-mag histogram offset", 	  lumfunc->offset);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.PEAK",  PS_META_REPLACE, "magnitude of peak bin",    	  lumfunc->mPeak);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.NUM.PEAK",  PS_META_REPLACE, "number of stars in peak bin", lumfunc->nPeak);
-	psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.SUM.PEAK",  PS_META_REPLACE, "sum of stars up to peak bin", lumfunc->sPeak);
-	psFree (lumfunc);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MIN",   PS_META_REPLACE, "min valid magnitude",             lumfunc->mMin);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MAX",   PS_META_REPLACE, "max valid magnitude",             lumfunc->mMax);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.SLOPE",  PS_META_REPLACE, "log-mag histogram slope",         lumfunc->slope);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.OFFSET", PS_META_REPLACE, "log-mag histogram offset",        lumfunc->offset);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.PEAK",  PS_META_REPLACE, "magnitude of peak bin",           lumfunc->mPeak);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.NUM.PEAK",  PS_META_REPLACE, "number of stars in peak bin", lumfunc->nPeak);
+        psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.SUM.PEAK",  PS_META_REPLACE, "sum of stars up to peak bin", lumfunc->sPeak);
+        psFree (lumfunc);
     }