Index: trunk/ppSim/src/ppSimLoop.c
===================================================================
--- trunk/ppSim/src/ppSimLoop.c	(revision 12907)
+++ trunk/ppSim/src/ppSimLoop.c	(revision 12917)
@@ -6,4 +6,5 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psastro.h>
 
 #include "ppSim.h"
@@ -268,28 +269,61 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
 
-    psVector *ra = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.RA"); // Star RAs
-    psVector *dec = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.DEC"); // Star Declinations
-    psVector *mag = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.MAG"); // Star magnitudes
+    psVector *ra = NULL;                // Star RAs
+    psVector *dec = NULL;               // Star Declinations
+    psVector *mag = NULL;               // Star magnitudes
     float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
-    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // Plate scale (arcsec/pixel)
-    scale *= M_PI / 3600.0 / 180.0; // Convert to radians/pixel
+    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") *
+        M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
     float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
 
     // Add catalogue stars
-    if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) {
-        // Convert star ra,dec to pixels from FPA centre
+    if (type == PPSIM_TYPE_OBJECT) {
         float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
         float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
         float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
 
+        // Read catalogue stars using psastro
+        psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE);
+        if (!astroRecipe) {
+            psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE);
+            psFree(rng);
+            psFree(bounds);
+            return PS_EXIT_CONFIG_ERROR;
+        }
+
+        // Size of FPA
+        float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;
+
+        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, "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);
+        psArray *refStars = psastroLoadRefstars(config);
+        if (!refStars || refStars->n == 0) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars.");
+            psFree(refStars);
+            psFree(rng);
+            psFree(bounds);
+            return PS_EXIT_CONFIG_ERROR;
+        }
+        psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);
+
+        ra = psVectorAlloc(refStars->n, PS_TYPE_F32);
+        dec = psVectorAlloc(refStars->n, PS_TYPE_F32);
+        mag = psVectorAlloc(refStars->n, PS_TYPE_F32);
+
         // Conversion loop
-        for (long i = 0; i < ra->n; i++) {
-            float div = sin(dec->data.F32[i]) * sin(dec0) +
-                cos(dec->data.F32[i]) * cos(dec0) * cos(ra->data.F32[i] - ra0); // Divisor
+        for (long i = 0; i < refStars->n; i++) {
+            pmAstromObj *ref = refStars->data[i]; // Reference star
+            float raStar = ref->sky->r; // RA of star
+            float decStar = ref->sky->d; // Dec of star
+            float magStar = ref->Mag;       // Magnitude of star
+
 
             // Convert to x,y position on tangent plane, in pixels
-            float xi = cos(dec->data.F32[i]) * sin(ra->data.F32[i] - ra0) / div / scale;
-            float eta = (sin(dec->data.F32[i]) * cos(dec0) -
-                cos(dec->data.F32[i]) * sin(dec0) * cos(ra->data.F32[i] - ra0)) / div / scale;
+            float div = sin(decStar) * sin(dec0) + cos(decStar) * cos(dec0) * cos(raStar - ra0); // Divisor
+            float xi = cos(decStar) * sin(raStar - ra0) / div / scale;
+            float eta = (sin(decStar) * cos(dec0) - cos(decStar) * sin(dec0) * cos(raStar - ra0)) /
+                div / scale;
 
             // Apply rotation
@@ -298,6 +332,7 @@
 
             // Convert magnitude to peak flux
-            mag->data.F32[i] = powf(10.0, -0.4 * (mag->data.F32[i] - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
-        }
+            mag->data.F32[i] = powf(10.0, -0.4 * (magStar - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
+        }
+        psFree(refStars);
 
         psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);