Changeset 10864 for trunk/psastro/src

Timestamp:

Jan 1, 2007, 11:03:03 AM (20 years ago)

Author:

eugene

Message:

added various plots, more work on the mosaic astrometry

Location:

trunk/psastro/src

Files:

• Makefile.am (modified) (2 diffs)
• psastro.h (modified) (1 diff)
• psastroAstromGuess.c (modified) (3 diffs)
• psastroDemoDump.c (added)
• psastroDemoPlot.c (added)
• psastroDemoPlots.c (deleted)
• psastroLoadRefstars.c (modified) (4 diffs)
• psastroMosaicAstrom.c (modified) (2 diffs)
• psastroMosaicSetAstrom.c (modified) (1 diff)
• psastroOneChip.c (modified) (3 diffs)

trunk/psastro/src/Makefile.am

@@ -5 +5 @@
 bin_PROGRAMS = psastro
 psastro_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
-psastro_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS)
+psastro_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
 psastro_LDADD = libpsastro.la
 
@@ -34 +34 @@
 psastroMosaicOneChip.c      \
 psastroMosaicSetAstrom.c    \
-psastroMosaicSetMatch.c
+psastroMosaicSetMatch.c     \
+psastroDemoDump.c           \
+psastroDemoPlot.c
 
 include_HEADERS = \

trunk/psastro/src/psastro.h

@@ -70 +70 @@
 psString psastroVersionLong(void);
 
+// demo plots
+bool psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip);
+bool psastroPlotRefstars (psArray *refstars);
+bool psastroPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, pmAstromFitResults *results);
+
+bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip);
+bool psastroDumpRefstars (psArray *refstars);
+
+bool psastroMosaicSetAstrom_tmp (pmFPA *fpa);

trunk/psastro/src/psastroAstromGuess.c

@@ -75 +75 @@
                 if (rawstars == NULL) { continue; }
 
-                if (psTraceGetLevel("psastro.dump.rawstars") > 2) {
-                  dump_stars (rawstars, fpa, chip);
-                }
-
-                if (psTraceGetLevel("psastro.plot.rawstars") > 2) {
-                  plot_stars (rawstars, fpa, chip);
-                }
-
                 for (int i = 0; i < rawstars->n; i++) {
                     pmAstromObj *raw = rawstars->data[i];
@@ -100 +92 @@
                     DECmax = PS_MAX (raw->sky->d, DECmax);
                 }
+
+                // dump or plot the resulting projected positions
+                if (psTraceGetLevel("psastro.dump") > 0) {
+                    psastroDumpRawstars (rawstars, fpa, chip);
+                }
+
+                if (psTraceGetLevel("psastro.plot") > 0) {
+                    psastroPlotRawstars (rawstars, fpa, chip);
+                }
             }
         }
@@ -125 +126 @@
    sky (ra, dec)
 */
-
-dump_stars (psArray *rawstars, pmFPA *fpa, pmChip *chip) {
-
-  FILE *f1 = fopen ("rawstars.up.dat", "w");
-  FILE *f2 = fopen ("rawstars.dn.dat", "w");
-  for (int i = 0; i < rawstars->n; i++) {
-    pmAstromObj *raw = rawstars->data[i];
-
-    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
-    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
-    psDeproject (raw->sky, raw->TP, fpa->toSky);
-
-    fprintf (f1, "%d  %f %f  %f %f  %f %f  %f %f\n", i,
-             raw->sky->r, raw->sky->d, 
-             raw->TP->x, raw->TP->y, 
-             raw->FP->x, raw->FP->y, 
-             raw->chip->x, raw->chip->y);
-                
-    psPlane *fp = psPlaneAlloc();
-    psPlane *tp = psPlaneAlloc();
-    psPlane *ch = psPlaneAlloc();
-                        
-    psProject (tp, raw->sky, fpa->toSky);
-    psPlaneTransformApply (fp, fpa->fromTPA, tp);
-    psPlaneTransformApply (ch, chip->fromFPA, fp);
-                        
-    fprintf (f2, "%d  %f %f  %f %f  %f %f  %f %f\n", i,
-             raw->sky->r, raw->sky->d, 
-             tp->x, tp->y, 
-             fp->x, fp->y, 
-             ch->x, ch->y);
-                
-    psFree (fp);
-    psFree (tp);
-    psFree (ch);
-  }
-
-  fclose (f1);
-  fclose (f2);
-}

trunk/psastro/src/psastroLoadRefstars.c

@@ -1 +1 @@
 # include "psastro.h"
-// int mkstemp(char *template);
 # define ELIXIR_MODE 1
 
@@ -84 +83 @@
     // convert the Average table to the pmAstromObj entries
     psTimerStart ("psastro");
-    psArray *refs = psArrayAllocEmpty (table->n);
+    psArray *refstars = psArrayAllocEmpty (table->n);
     for (int i = 0; i < table->n; i++) {
         pmAstromObj *ref = pmAstromObjAlloc ();
@@ -101 +100 @@
         }
 
-        psArrayAdd (refs, 100, ref);
+        psArrayAdd (refstars, 100, ref);
         psFree (ref);
     }
@@ -109 +108 @@
 
     psTrace ("psastro", 3, "loaded %ld reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n",
-             refs->n, RAmin, DECmin, RAmax, DECmax);
+             refstars->n, RAmin, DECmin, RAmax, DECmax);
 
-    return refs;
+    // dump or plot the available refstars
+    if (psTraceGetLevel("psastro.dump") > 0) {
+        psastroDumpRefstars (refstars);
+    }
+
+    if (psTraceGetLevel("psastro.plot") > 0) {
+        psastroPlotRefstars (refstars);
+    }
+
+    return refstars;
 }

trunk/psastro/src/psastroMosaicAstrom.c

@@ -32 +32 @@
 
     // fitted chips will follow the local plate-scale, hiding the distortion
-    // modify the chip->toFPA scaling to match knowledge about pixel scale
+    // modify the chip->toFPA scaling to match knowledge about pixel scale,
+    // then recalculate raw and ref positions
     psastroMosaicCommonScale (fpa, recipe);
+    psastroMosaicSetAstrom (fpa);
 
     gradients = psastroMosaicGradients (fpa, recipe);
@@ -62 +64 @@
     // XXX modify match radius
     // XXX set chip.order to 1
-    psastroMosaicSetAstrom (fpa);
+    psastroMosaicSetAstrom_tmp (fpa);
     psastroMosaicSetMatch (fpa, recipe, 1);
     psastroMosaicChipAstrom (fpa, recipe, 1);

trunk/psastro/src/psastroMosaicSetAstrom.c

@@ -49 +49 @@
     return true;
 }
+
+bool psastroMosaicSetAstrom_tmp (pmFPA *fpa) {
+
+    pmChip *chip = NULL;
+    pmCell *cell = NULL;
+    pmReadout *readout = NULL;
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    FILE *f1 = fopen ("raw.tmp.dat", "w");
+    FILE *f2 = fopen ("ref.tmp.dat", "w");
+
+    // this loop selects the matched stars for all chips
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+        
+        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
+            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+            // process each of the readouts
+            // XXX there can only be one readout per chip, right?
+            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
+                if (! readout->data_exists) { continue; }
+
+                // select the raw objects for this readout
+                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
+                if (rawstars == NULL) { continue; }
+
+                for (int i = 0; i < rawstars->n; i++) {
+                    pmAstromObj *raw = rawstars->data[i];
+        
+                    fprintf (f1, "%f %f  %f %f  %f %f  ", raw->sky->r, raw->sky->d, raw->TP->x, raw->TP->y, raw->FP->x, raw->FP->y);
+
+                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
+                    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
+                    psDeproject (raw->sky, raw->TP, fpa->toSky);
+
+                    fprintf (f1, "|  %f %f  %f %f  %f %f   %f\n", raw->sky->r, raw->sky->d, raw->TP->x, raw->TP->y, raw->FP->x, raw->FP->y, raw->Mag);
+                }
+
+                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
+                if (refstars == NULL) { continue; }
+
+                for (int i = 0; i < refstars->n; i++) {
+                    pmAstromObj *ref = refstars->data[i];
+        
+                    fprintf (f2, "%f %f  %f %f  %f %f  ", ref->sky->r, ref->sky->d, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y);
+
+                    psProject (ref->TP, ref->sky, fpa->toSky);
+                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
+                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+
+                    fprintf (f2, "|  %f %f  %f %f  %f %f   %f\n", ref->sky->r, ref->sky->d, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y, ref->Mag);
+                }
+            }
+        }
+    }
+    fclose (f1);
+    fclose (f2);
+    return true;
+}

trunk/psastro/src/psastroOneChip.c

@@ -90 +90 @@
     }
     
-    // toTPA converts from FPA units (microns) to TPA units (linear degrees)
-    float plateScale1 = hypot (fpa->toTPA->x->coeff[1][0], fpa->toTPA->x->coeff[0][1]);
-    float plateScale2 = hypot (fpa->toTPA->y->coeff[1][0], fpa->toTPA->y->coeff[0][1]);
-    float plateScale = 0.5*(plateScale1 + plateScale2)*3600.0;
+    // toSky converts from FPA & TPA units (microns) to sky units (radians)
+    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
 
     // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
@@ -105 +103 @@
 
     // XXX should these result in errors or be handled another way?
-    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
     if (astError > maxError) {
         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
@@ -127 +125 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
 
-// write results
+    // write results
     psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
+    
+    if (psTraceGetLevel("psastro.plot") > 0) {
+        psastroPlotOneChipFit (rawstars, refstars, match, results);
+    }
 
     psFree (match);