Changeset 27609

Timestamp:

Apr 5, 2010, 3:30:27 PM (16 years ago)

Author:

Paul Price

Message:

Allow reading of streaks file.

Location:

trunk/ppViz/src/ppCoord

Files:

• ppCoord.h (modified) (1 diff)
• ppCoordArguments.c (modified) (3 diffs)
• ppCoordData.c (modified) (2 diffs)
• ppCoordLoop.c (modified) (4 diffs)

trunk/ppViz/src/ppCoord/ppCoord.h

@@ -16 +16 @@
     psString chipName;                  // Name of chip of interest
     psString radecName;                 // Filename with sky coordinates
+    psString streaksName;               // Filename with streaks (sky coordinates)
     pmConfig *config;                   // Configuration
     bool radians;                       // RA,Dec are in radians?

trunk/ppViz/src/ppCoord/ppCoordArguments.c

@@ -51 +51 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-chip", 0, "Chip for pixel coordinates", NULL);
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-radec", 0, "Filename with RA,Dec (default decimal degrees)", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-streaks", 0, "Filename with streaks", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-radians", 0, "RA,Dec in radians?", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-all", 0, "Output all coordinates?", NULL);
@@ -65 +66 @@
     data->chipName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-chip"));
     data->radecName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-radec"));
+    data->streaksName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-streaks"));
     data->radians = psMetadataLookupBool(NULL, arguments, "-radians");
     data->all = psMetadataLookupBool(NULL, arguments, "-all");
@@ -88 +90 @@
     }
 
-    if (!data->pixelsName && !data->radecName) {
+    if (!data->pixelsName && !data->radecName && !data->streaksName) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Neither -pixels nor -radec provided.");
         return false;

trunk/ppViz/src/ppCoord/ppCoordData.c

@@ -18 +18 @@
     psFree(data->chipName);
     psFree(data->radecName);
+    psFree(data->streaksName);
     psFree(data->config);
     psFree(data->ds9name);
@@ -38 +39 @@
     data->chipName = NULL;
     data->radecName = NULL;
+    data->streaksName = NULL;
     data->config = NULL;
     data->radians = false;

trunk/ppViz/src/ppCoord/ppCoordLoop.c

@@ -23 +23 @@
         NULL; // File with raw image
 
-    psArray *pixels = NULL, *radec = NULL; // Array of vectors with coordinates
-    psArray *radecOut = NULL;              // Output for sky coordinates
+    psArray *pixels = NULL, *radec = NULL, *streaks = NULL; // Array of vectors with coordinates
+    psArray *radecOut = NULL, *streaksOut = NULL;           // Output for sky coordinates
     if (data->pixelsName) {
         pixels = psVectorsReadFromFile(data->pixelsName, "%f %f");
@@ -47 +47 @@
         psVectorInit(radecOut->data[1], NAN);
         psVectorInit(radecOut->data[2], NAN);
+    }
+
+    if (data->streaksName) {
+        FILE *streaksFile = fopen(data->streaksName, "r"); // File handle for streaks
+        if (!streaksFile) {
+            psError(PS_ERR_IO, true, "Unable to open streaks file %s", data->streaksName);
+            return false;
+        }
+        int numStreaks = 0;             // Number of streaks
+        if (fscanf(streaksFile, "%d", &numStreaks) != 1) {
+            psError(PS_ERR_IO, true, "Unable to read number of streaks from %s", data->streaksName);
+            return false;
+        }
+
+        streaks = psArrayAlloc(4);
+        psVector *ra1 = streaks->data[0] = psVectorAlloc(numStreaks, PS_TYPE_F64);
+        psVector *dec1 = streaks->data[1] = psVectorAlloc(numStreaks, PS_TYPE_F64);
+        psVector *ra2 = streaks->data[2] = psVectorAlloc(numStreaks, PS_TYPE_F64);
+        psVector *dec2 = streaks->data[3] = psVectorAlloc(numStreaks, PS_TYPE_F64);
+
+        for (int i = 0; i < numStreaks; i++) {
+            if (fscanf(streaksFile, "%lf %lf %lf %lf %*d",
+                       &ra1->data.F64[i], &dec1->data.F64[i], &ra2->data.F64[i], &dec2->data.F64[i]) != 4) {
+                psError(PS_ERR_IO, true, "Unable to read streak %d of %d from %s",
+                        i, numStreaks, data->streaksName);
+                return false;
+            }
+        }
+        streaksOut = psArrayAlloc(4);
+        streaksOut->data[0] = psVectorAlloc(numStreaks, PS_TYPE_F32);
+        streaksOut->data[1] = psVectorAlloc(numStreaks, PS_TYPE_F32);
+        streaksOut->data[2] = psVectorAlloc(numStreaks, PS_TYPE_F32);
+        streaksOut->data[3] = psVectorAlloc(numStreaks, PS_TYPE_F32);
+        psVectorInit(streaksOut->data[0], NAN);
+        psVectorInit(streaksOut->data[1], NAN);
+        psVectorInit(streaksOut->data[2], NAN);
+        psVectorInit(streaksOut->data[3], NAN);
     }
 
@@ -282 +319 @@
         }
 
+        if (streaks && data->chipName && !rawFile) {
+            psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1
+            psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1
+            psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2
+            psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2
+
+            psVector *ra1 = streaks->data[0]; // RA coordinate for point 1
+            psVector *dec1 = streaks->data[1]; // Dec coordinate for point 1
+            psVector *ra2 = streaks->data[2]; // RA coordinate for point 2
+            psVector *dec2 = streaks->data[3]; // Dec coordinate for point 2
+
+            psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
+            psPlane *fp = psPlaneAlloc();    // Focal plane coordinates
+            psPlane *tp = psPlaneAlloc();    // Tangent plane coordinates
+            psSphere *sky = psSphereAlloc(); // Sky coordinates
+            for (long i = 0; i < xPix1->n; i++) {
+                sky->r = ra1->data.F64[i];
+                sky->d = dec1->data.F64[i];
+                psProject(tp, sky, astromFile->fpa->toSky);
+                psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp);
+                psPlaneTransformApply(pix, chip->fromFPA, fp);
+                xPix1->data.F32[i] = pix->x;
+                yPix1->data.F32[i] = pix->y;
+
+                sky->r = ra2->data.F64[i];
+                sky->d = dec2->data.F64[i];
+                psProject(tp, sky, astromFile->fpa->toSky);
+                psPlaneTransformApply(fp, astromFile->fpa->fromTPA, tp);
+                psPlaneTransformApply(pix, chip->fromFPA, fp);
+                xPix2->data.F32[i] = pix->x;
+                yPix2->data.F32[i] = pix->y;
+            }
+        }
+
         // Chip
         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
@@ -323 +394 @@
     }
 
+    if (streaksOut) {
+        psVector *xPix1 = streaksOut->data[0]; // x coordinate for point 1
+        psVector *yPix1 = streaksOut->data[1]; // y coordinate for point 1
+        psVector *xPix2 = streaksOut->data[2]; // x coordinate for point 2
+        psVector *yPix2 = streaksOut->data[3]; // y coordinate for point 2
+
+        psVector *ra1 = streaks->data[0]; // RA coordinate for point 1
+        psVector *dec1 = streaks->data[1]; // Dec coordinate for point 1
+        psVector *ra2 = streaks->data[2]; // RA coordinate for point 2
+        psVector *dec2 = streaks->data[3]; // Dec coordinate for point 2
+
+        for (long i = 0; i < xPix1->n; i++) {
+            if (data->ds9) {
+                fprintf(data->ds9, "image;line(%f,%f,%f,%f) # line= 0 0 color=%s\n",
+                        xPix1->data.F32[i], yPix1->data.F32[i], xPix2->data.F32[i], yPix2->data.F32[i],
+                        data->ds9color);
+            } else {
+                fprintf(stdout, "%.10lf %.10lf %.10lf %.10lf --> %.3f %.3f %.3f %.3f\n",
+                        ra1->data.F64[i], dec1->data.F64[i],
+                        ra2->data.F64[i], dec2->data.F64[i],
+                        xPix1->data.F32[i], yPix1->data.F32[i],
+                        xPix2->data.F32[i], yPix2->data.F32[i]);
+            }
+        }
+    }
+
     psFree(pixels);
     psFree(radec);
     psFree(radecOut);
+    psFree(streaks);
+    psFree(streaksOut);
 
     return true;