Changeset 5506 for trunk/psastro/src/psastroBuildFPA.c

Timestamp:

Nov 11, 2005, 9:54:59 AM (21 years ago)

Author:

eugene

Message:

adding WCS interpretation code

File:

• trunk/psastro/src/psastroBuildFPA.c (modified) (1 diff)

trunk/psastro/src/psastroBuildFPA.c

@@ -1 +1 @@
 # include "psastro.h"
 
+pmFPA *psastroBuildFPA (psMetadata *header, psMetadata *config) {
+
+    // this function constructs a basic template pmFPA structure based on the information in the header
+    // metadata.  it converts the astrometric information in the header (WCS, etc) to the corresponding
+    // pmFPA/pmChip, etc distortion and projection elements.
+    
+    // this implementation is basic, assuming a single readout,chip,cell
+    // I will not handle the datasec/biassec information to define the readout/cell
+    // these things are already handled by Paul's code.  this function simple looks for the WCS information
+    // and provides those conversions.
+
+    int Nchips    = 1;
+    int Ncells    = 1;
+    int Nreadouts = 1;
+
+    // allocate the structures
+
+    pmFPA *fpa = pmFPAAlloc ();
+
+    fpa->chips = psArrayAlloc (Nchips);
+    for (int i = 0; i < Nchips; i++) {
+
+        pmChip *chip = pmChipAlloc ();
+        chip->fpa = fpa; // assign parent fpa (view only; don't free)
+
+        chip->cells = psArrayAlloc (Ncells);
+        for (int j = 0; j < Ncells; j++) {
+
+            pmCell *cell = pmCellAlloc ();
+            cell->chip = chip; // assign parent chip (view only; don't free)
+            cell->header = header;
+
+            pmCellInterpretWCS (cell, header);
+
+            cell->readouts = psArrayAlloc (Nreadouts);
+            for (int k = 0; k < Nreadouts; k++) {
+
+                pmReadout *readout = pmReadoutAlloc ();
+                cells->readouts->data[j] = readout;
+            }
+            chips->cells->data[j] = cell;
+        }
+        fpa->chips->data[i] = chip;
+    }
+    return (fpa);
+}
+
+// interpret header WCS
+bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) { 
+
+    float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
+    float pc1_1, pc1_2, pc2_1, pc2_2;
+
+    // *** interpret header data, convert to crval(i), etc
+    char *ctype = pmMetadataLookupPtr (&status, header, "CTYPE2");
+    if (!status) {
+        psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n");
+        return false;
+    }
+
+    // determine projection type
+    type = PS_PROJ_NTYPE;
+    if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN;
+    if (!strcmp (&ctype[4], "-TAN")) type = PS_PROJ_TAN;
+    if (!strcmp (&ctype[4], "-AIT")) type = PS_PROJ_AIT;
+    if (!strcmp (&ctype[4], "-PAR")) type = PS_PROJ_PAR;
+    if (type == PS_PROJ_NTYPE) {
+        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
+        return false;
+    }
+
+    crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");
+    crval2 = psMetadataLookupF32 (&status, header, "CRVAL2");
+    crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1");
+    crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2");
+    
+    // test the CDELTi varient
+    cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1");
+    if (status) {
+        cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2");
+
+        // test the CROTAi varient:
+        float rotate = psMetadataLookupF32 (&status, header, "CROTA2");
+        if (status) {
+            Lambda = cdelt2 / cdelt1;
+            pc1_1 =  cos(rotate*RAD_DEG);
+            pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
+            pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
+            pc2_2 =  cos(rotate*RAD_DEG);
+            goto got_matrix;
+        }
+
+        // test the PC00i00j varient:
+        pc1_1 = psMetadataLookupF32 (&status, header, "PC001001");
+        if (status) {
+            pc1_2 = psMetadataLookupF32 (&status, header, "PC001002");
+            pc2_1 = psMetadataLookupF32 (&status, header, "PC002001");
+            pc2_2 = psMetadataLookupF32 (&status, header, "PC002002");
+
+            // XXX EAM : add Elixir polynomial terms here eventually
+            goto got_matrix;
+        }
+        psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");
+        return false;
+    }
+
+    // test the CDi_j varient
+    pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1");
+    if (status) {
+        pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2");
+        pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1");
+        pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2");
+        
+        // renormalize to cdelt1, cdelt2, etc
+        float scale = hypot (pc1_1, pc1_2);
+        cdelt1 = cdelt2 = scale;
+        pc1_1 /= scale;
+        pc1_2 /= scale;
+        pc2_1 /= scale;
+        pc2_2 /= scale;
+        goto got_matrix;
+    }
+    psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");
+    return false;
+
+got_matrix:
+
+    // XXX EAM : these must already be set
+    chip = cell->chip;
+    fpa = chip->fpa;
+
+    // set toChip to identity as default
+    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
+    cell->toChip   = psPlaneTransformAlloc (2, 2);
+    cell->toChip->x->coeff[1][0] = 1;
+    cell->toChip->x->mask[1][1]  = 1;
+
+    cell->toChip->y->coeff[0][1] = 1;
+    cell->toChip->y->mask[1][1]  = 1;
+
+    // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
+    //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
+
+    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
+    chip->toFPA   = psPlaneTransformAlloc (2, 2);
+    
+    chip->toFPA->x->coeff[0][0] = crpix1;
+    chip->toFPA->x->coeff[1][0] = pc1_1;
+    chip->toFPA->x->coeff[0][1] = pc1_2;
+    chip->toFPA->x->mask[1][1]  = 1;
+
+    chip->toFPA->y->coeff[0][0] = crpix2;
+    chip->toFPA->y->coeff[1][0] = pc2_1;
+    chip->toFPA->y->coeff[0][1] = pc2_2;
+    chip->toFPA->y->mask[1][1]  = 1;
+
+    // set toTPA to identity as default
+    // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581)
+    if (fpa->toTPA == NULL) {
+        fpa->toTPA   = psPlaneDistortAlloc (2, 2, 1, 1);
+        fpa->toTPA->x->coeff[1][0][0][0] = 1;
+        fpa->toTPA->x->mask[1][1][0][0]  = 1;
+
+        fpa->toTPA->y->coeff[0][1][0][0] = 1;
+        fpa->toTPA->y->mask[1][1][0][0]  = 1;
+    } else {
+        psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n");
+    }
+
+    // center of projection is (0,0) coordinate of TPA
+    fpa->toSky = psProjectionAlloc (crval1, crval2, cdelt1, cdelt2, type);
+    return true;
+}