Changeset 39685 for branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c

Timestamp:

Sep 11, 2016, 2:29:33 PM (10 years ago)

Author:

eugene

Message:

add photometry difference weighting for the astrometric fit to the psastro.config recipe; add ZPN projection to the WCS I/O

File:

• branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c (modified) (5 diffs)

branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c

@@ -277 +277 @@
     wcs->crpix2 = psMetadataLookupF64 (&status, header, "CRPIX2");
     wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
+
+    // XXX if type == ZPN, look for PV2_%d elements:
+    if (type == PS_PROJ_ZPN) {
+      psVector *maxRadial = psVectorAlloc (21, PS_TYPE_F64);
+      for (int i = 0; i <= 20; i++) {
+        char name[64];
+        snprintf (name, 64, "PV2_%d", i);
+
+        maxRadial->data.F64[i] = 0.0;
+        double value = psMetadataLookupF64 (&status, header, name);
+
+        if (status) {
+          maxRadial->data.F64[i] = value;
+          maxRadial->n = i;
+        }
+
+        // PV2_1 is implicit if not present
+        if ((i == 1) && !status) {
+          maxRadial->data.F64[i] = 1.0;
+          continue;
+        }
+      }
+      maxRadial->n ++;
+      wcs->toSky->radial = maxRadial;
+    }
 
     // These aren't needed but having them empty is disconcerting
@@ -375 +400 @@
     psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", wcs->crpix1);
     psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", wcs->crpix2);
+
+    if (wcs->toSky->type == PS_PROJ_ZPN) {
+      psAssert (wcs->toSky->radial, "missing radial vector");
+      for (int i = 0; i < wcs->toSky->radial->n; i++) {
+        if (wcs->toSky->radial->data.F64[i] == 0.0) continue;
+        snprintf (name, 16, "PV2_%d", i);
+        psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->toSky->radial->data.F64[i]);
+      }
+    }
 
     // XXX make it optional to write out CDi_j terms, or other versions
@@ -485 +519 @@
     // projection from TPA (linear microns) to SKY (radians)
     psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG*pdelt1, PM_RAD_DEG*pdelt2, wcs->toSky->type);
+    toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
 
     if (fpa->toSky == NULL) {
@@ -633 +668 @@
     // cdelt1,2 has units of degrees/micron
     fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->cdelt1*PM_RAD_DEG, wcs->cdelt2*PM_RAD_DEG, wcs->toSky->type);
+    fpa->toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
 
     // create transformation with 0,0 reference pixel
@@ -701 +737 @@
     // convert projection from FPA to SKY into wcs projection (degrees to radians)
     wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
+    wcs->toSky->radial = psMemIncrRefCounter (fpa->toSky->radial);
+
     wcs->crval1 = fpa->toSky->R*PS_DEG_RAD;
     wcs->crval2 = fpa->toSky->D*PS_DEG_RAD;