Changeset 38183

Timestamp:

Apr 24, 2015, 3:40:47 PM (11 years ago)

Author:

eugene

Message:

add PROJ_ZPN (currently only valid for UKIRT data with P1 & P3 non-zero)

Location:

branches/eam_branches/ipp-20150419/Ohana/src/libdvo

Files:

• include/libdvo_astro.h (modified) (1 diff)
• src/coordops.c (modified) (12 diffs)

branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/libdvo_astro.h

@@ -26 +26 @@
   PROJ_GLS, // pseudocyl
   PROJ_PAR, // pseudocyl
+  PROJ_ZPN, // zenithal
 } OhanaProjection;
 

branches/eam_branches/ipp-20150419/Ohana/src/libdvo/src/coordops.c

@@ -1 +1 @@
 # include <dvo.h>
+# define UKIRT_ONLY 1
 
 /* note that Coords.ctype carries the DEC (ctype2) value */
@@ -81 +82 @@
 
   /** extra polynomial terms **/
-  if (coords[0].Npolyterms > 1) {
+  if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
     X2 = X*X;
     Y2 = Y*Y;
@@ -171 +172 @@
         stht = cos (RAD_DEG * R);
         break;
+
+      case PROJ_ZPN:
+        // R = 90 - theta (degrees)
+        // this is wrong because we are ignoring the distortion
+        // XXX For now, just solve for terms up to n = 3
+
+        // Ro = (pi/180)(90 - theta)
+        // R = (180/pi)sum (P_i R^i)
+
+        if (UKIRT_ONLY) {
+          double Ro = RAD_DEG * R;
+          double P1 = coords[0].polyterms[0][1];
+          double P3 = coords[0].polyterms[0][3];
+
+          // find the roots of f(gamma) = P1 gamma + P3 gamma^3 - Ro
+          // starting guess for gamma is Ro / P1
+          double gamma = Ro / P1;
+
+          int i;
+          for (i = 0; i < 5; i++) {
+            double F = P1*gamma + P3*gamma*gamma*gamma - Ro;
+            double dFdgamma = P1 + 3.0*P3*gamma*gamma;
+
+            double gamma_new = gamma - F / dFdgamma;
+            gamma = gamma_new;
+          }
+          
+          double theta = 90.0 - gamma * DEG_RAD ;
+          ctht = cos (RAD_DEG * theta);
+          stht = sin (RAD_DEG * theta);
+          break;
+        }
+
       case PROJ_ZEA:
       case PROJ_ZPL:
@@ -238 +272 @@
 int RD_to_LM (double *L, double *M, double ra, double dec, Coords *coords) {
 
+  int i;
   double phi, theta;
   double sphi, cphi, stht, ctht;
@@ -303 +338 @@
         *M = -DEG_RAD * cphi;
         return (stht > 0);
+
       case PROJ_ARC:
         // R = 90 - theta
@@ -311 +347 @@
         *M = (ctht == 0.0) ? 0.0 : -Rc * cphi / ctht ;
         return (TRUE);
+
+      case PROJ_ZPN:
+        // Ro = (pi/180)(90 - theta)
+        // R = (180/pi)sum (P_i R^i)
+        ctht = hypot(sphi, cphi);
+        theta = atan2 (stht, ctht);
+
+        double Rc;
+        if (UKIRT_ONLY) {
+          double P1 = coords[0].polyterms[0][1];
+          double P3 = coords[0].polyterms[0][3];
+          double gamma = RAD_DEG * (90 - DEG_RAD * theta);
+          Rc = P1*gamma + P3*gamma*gamma*gamma;
+        } else {
+          double Ro = RAD_DEG * (90 - DEG_RAD * theta);
+
+          // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
+          i = coords[0].Npolyterms - 1;
+          Rc = 0.0;
+          while (i > 0) {
+            if (i < 7) {
+              Rc = (Rc + coords[0].polyterms[0][i])*Ro;
+            } else {
+              Rc = (Rc + coords[0].polyterms[1][i-7])*Ro;
+            }
+            i --;
+          }
+          Rc += coords[0].polyterms[0][i];
+        }
+        Rc = DEG_RAD * Rc;
+
+        *L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ;
+        *M = (ctht == 0.0) ? 0.0 : -Rc * cphi / ctht ;
+        return (TRUE);
+
       case PROJ_ZEA:
       case PROJ_ZPL:
@@ -386 +457 @@
 
   /** extra polynomial terms **/
-  if (coords[0].Npolyterms > 1) {
+  // polyterm values are only valid for some projections.  for PROJ_ZPN, they are applied to R, not X,Y
+  if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
     for (i = 0; i < 10; i++) {
       // find derivatives at Xo, Yo
@@ -531 +603 @@
 int GetCoords (Coords *coords, Header *header) {
   
-  int status, status1, status2, itmp, Polynomial, Polyterm;
+  int i, status, status1, status2, itmp, Polynomial, Polyterm;
   double Lambda, rotate, rotate1, rotate2, scale;
   double equinox;
@@ -663 +735 @@
       coords[0].pc2_1 /= scale;
       coords[0].pc2_2 /= scale;
+
+      if (!strcmp (&coords[0].ctype[4], "-ZPN")) {
+        int found;
+        for (i = 0; i < 14; i++) {
+          char name[64];
+          sprintf (name, "PV2_%d", i);
+          if (i < 7) {
+            found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[0][i]);
+          } else {
+            found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[1][i-7]);
+          }
+          if ((i == 0) && !found) {
+            coords[0].polyterms[0][0] = 0.0;
+            continue;
+          }
+          if ((i == 1) && !found) {
+            coords[0].polyterms[0][1] = 1.0;
+            continue;
+          }
+          if (!found) {
+            coords[0].Npolyterms = i;
+            break;
+          }
+        }
+      }
       break;
 
@@ -895 +992 @@
   if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA;
   if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL;
+  if (!strcmp(&ctype[4], "-ZPN")) return PROJ_ZPN;
   if (!strcmp(&ctype[4], "-ARC")) return PROJ_ARC;
   if (!strcmp(&ctype[4], "-STG")) return PROJ_STG;
@@ -916 +1014 @@
     case PROJ_ZEA: strcpy(&ctype[4], "-ZEA"); return TRUE;
     case PROJ_ZPL: strcpy(&ctype[4], "-ZPL"); return TRUE;
+    case PROJ_ZPN: strcpy(&ctype[4], "-ZPN"); return TRUE;
     case PROJ_ARC: strcpy(&ctype[4], "-ARC"); return TRUE;
     case PROJ_STG: strcpy(&ctype[4], "-STG"); return TRUE;
@@ -937 +1036 @@
     case PROJ_ZEA:
     case PROJ_ZPL:
+    case PROJ_ZPN:
     case PROJ_ARC:
     case PROJ_STG: