Changeset 7332 for trunk/psastro/src/psastroMosaicGetRefstars.c

Timestamp:

Jun 4, 2006, 7:02:16 AM (20 years ago)

Author:

eugene

Message:

substantial updates

File:

• trunk/psastro/src/psastroMosaicGetRefstars.c (modified) (3 diffs)

trunk/psastro/src/psastroMosaicGetRefstars.c

@@ -1 +1 @@
 # include "psastro.h"
 
-bool psastroMosaicGetRefstars (pmFPA *fpa, psArray *refs, psMetadata *recipe) {
+bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs) {
 
+    bool status;
     pmChip *chip = NULL;
     pmCell *cell = NULL;
     pmReadout *readout = NULL;
+
+    // select the current recipe
+    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
+    if (!recipe) {
+        psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // select the input data sources
+    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
+    if (!input) {
+        psErrorStackPrint(stderr, "Can't find input data!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA"); 
+    if (!status) fieldExtra = 0.0;
+
     pmFPAview *view = pmFPAviewAlloc (0);
+    pmFPA *fpa = input->fpa;
 
     // this loop selects the matched stars for all chips
@@ -22 +42 @@
                 if (! readout->data_exists) { continue; }
 
-                // select the raw objects for this readout
-                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
-                if (rawstars == NULL) { continue; }
+                // read WCS data from the corresponding header
+                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
 
+                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1"); 
+                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2"); 
+
+                float minX = -fieldExtra*Nx;
+                float maxX = (1+fieldExtra)*Nx;
+                float minY = -fieldExtra*Ny;
+                float maxY = (1+fieldExtra)*Ny;
+                
                 // the refstars is a subset within range of this chip
                 psArray *refstars = psArrayAlloc (100);
+
+                int Nst = 0;
 
                 // select the reference objects within range of this readout
                 // project the reference objects to this chip
                 for (int i = 0; i < refs->n; i++) {
-                    pmAstromObj *ref = refs->data[i];
+                    pmAstromObj *ref = pmAstromObjCopy(refs->data[i]);
         
                     p_psProject (ref->TP, ref->sky, fpa->projection);
@@ -38 +67 @@
                     psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
 
-                    // XXX what is the X,Y range of the selected chip?
-                    // XXX for the moment, force these to be good for Megacam
-                    if (ref->chip->x <  -10) continue;
-                    if (ref->chip->x > 2060) continue;
-                    if (ref->chip->y <  -10) continue;
-                    if (ref->chip->y > 4610) continue;
+                    // limit the X,Y range of the refs to the selected chip
+                    if (ref->chip->x < minX) continue;
+                    if (ref->chip->x > maxX) continue;
+                    if (ref->chip->y < minY) continue;
+                    if (ref->chip->y > maxY) continue;
                     
+                    if (Nst < 0) {
+                        fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n", 
+                                 ref->chip->x, ref->chip->y, 
+                                 ref->FP->x, ref->FP->y, 
+                                 ref->TP->x, ref->TP->y, 
+                                 ref->sky->r, ref->sky->d);
+
+                        psSphere *sk = psSphereAlloc();
+                        psPlane *fp = psPlaneAlloc();
+                        psPlane *tp = psPlaneAlloc();
+
+                        psPlaneTransformApply (fp, chip->toFPA, ref->chip);
+                        psPlaneDistortApply (tp, fpa->toTangentPlane, fp, 0.0, 0.0);
+                        p_psDeproject (sk, tp, fpa->projection);
+
+                        fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n", 
+                                 ref->chip->x, ref->chip->y, 
+                                 fp->x, fp->y, 
+                                 tp->x, tp->y, 
+                                 sk->r, sk->d);
+                    }
+                    Nst ++;
+
                     psArrayAdd (refstars, 100, ref);
                 }
+                psTrace (__func__, 4, "Added %d refstars\n", refstars->n);
+
                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
+
+                // XXX test point
+                if (0) { 
+                    pmChip *chip1 = fpa->chips->data[10];
+                    pmCell *cell1 = chip1->cells->data[0];
+                    pmReadout *readout1 = cell1->readouts->data[0];
+                    psArray *refstars1 = psMetadataLookupPtr (NULL, readout1->analysis, "PSASTRO.REFSTARS");
+                    pmAstromObj *ref1 = refstars1->data[0];
+                    fprintf (stderr, "test 3: %f %f  %f %f  %f %f\n", ref1->chip->x, ref1->chip->y, ref1->FP->x, ref1->FP->y, ref1->sky->r, ref1->sky->d);
+                }
+
+                psFree (refstars);
             }
         }
     }
+
+    psFree (view);
     return true;
 }