Changeset 12917

Timestamp:

Apr 18, 2007, 4:11:34 PM (19 years ago)

Author:

Paul Price

Message:

Using psastro to get the reference stars, instead of a text file.

Location:

trunk/ppSim

Files:

• configure.ac (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/ppSimArguments.c (modified) (3 diffs)
• src/ppSimLoop.c (modified) (3 diffs)

trunk/ppSim/configure.ac

@@ -19 +19 @@
 PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0])
 PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0])
+PKG_CHECK_MODULES([PSASTRO], [psastro >= 0.9.0])
 
 IPP_STDOPTS

trunk/ppSim/src/Makefile.am

@@ -1 +1 @@
 bin_PROGRAMS = ppSim
 
-ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(ppSim_CFLAGS)
-ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS)
+ppSim_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS) $(ppSim_CFLAGS)
+ppSim_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
 ppSim_SOURCES = \
         ppSim.c                 \

trunk/ppSim/src/ppSimArguments.c

@@ -170 +170 @@
                      config->argv[1]);
 
-
-    // Values common for stars
-    float scale = psMetadataLookupF32(NULL, arguments, "-scale"); // Plate scale
-    float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point
-    float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // Zero point
-    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale);
-    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp);
-    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SEEING", 0, "Seeing simga (pix)",
-                     seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);
-
-    const char *starName = psMetadataLookupStr(NULL, arguments, "-stars"); // Name of file with stars
-    if (starName) {
+    if (type == PPSIM_TYPE_OBJECT) {
+        // Load values required for adding stars
+        float scale = psMetadataLookupF32(NULL, arguments, "-scale"); // Plate scale
+        float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point
+        float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // Zero point
         float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight
         float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight
@@ -188 +181 @@
         if (isnan(ra0) || isnan(dec0) || isnan(pa) || isnan(scale) || isnan(zp) || isnan(seeing)) {
             psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                    "-ra, -dec, -pa -scale, -zp, -seeing must be specified with -stars");
+                    "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type");
             usage(arguments, config);
         }
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale);
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp);
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SEEING", 0, "Seeing simga (pix)",
+                         seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);
         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)",
                          ra0 * M_PI / 180.0);
@@ -197 +194 @@
         psMetadataAddF32(config->arguments, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",
                          pa * M_PI / 180.0);
-
-        psArray *raDecMag = psVectorsReadFromFile(starName, "%f %f %f"); // Vectors for RA Dec Mag
-        if (!raDecMag || raDecMag->n != 3) {
-            psError(PS_ERR_IO, false, "Unable to read RA Dec Mag from %s", starName);
-            psFree(raDecMag);
-            psFree(arguments);
-            psFree(config);
-            exit(PS_EXIT_CONFIG_ERROR);
-        }
-
-        psVector *ra = raDecMag->data[0]; // Star RAs
-        psVector *dec = raDecMag->data[1]; // Star Declinations
-        psVector *mag = raDecMag->data[2]; // Star magnitudes
-        for (long i = 0; i < ra->n; i++) {
-            ra->data.F32[i] *= M_PI / 180.0;
-            dec->data.F32[i] *= M_PI / 180.0;
-        }
-
-        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.RA", 0, "Star RAs (radians)", ra);
-        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.DEC", 0, "Star Decs (radians)", dec);
-        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.MAG", 0, "Star magnitudes", mag);
-        psFree(raDecMag);
-    }
-
-    if (psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY") > 0 &&
-        (isnan(scale) || isnan(zp) || isnan(seeing))) {
-        psWarning("-scale, -zp, -seeing must be specified if fake stars are to be generated --- turned off.");
-        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STASR.NUM", PS_META_REPLACE,
-                         "Fake stars turned off", 0);
-    }
-
+    }
 
     psFree(arguments);

trunk/ppSim/src/ppSimLoop.c

@@ -6 +6 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psastro.h>
 
 #include "ppSim.h"
@@ -268 +269 @@
     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 +332 @@
 
             // 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);