Changeset 41966 for branches/eam_branches/ipp-20211108

Timestamp:

Dec 25, 2021, 2:19:38 PM (5 years ago)

Author:

eugene

Message:

add MOL projection

Location:

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

Files:

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

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

@@ -27 +27 @@
   PROJ_PAR, // pseudocyl
   PROJ_ZPN, // zenithal
+  PROJ_MOL, // pseudocyl 
 } OhanaProjection;
 

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

@@ -266 +266 @@
     switch (proj) {
       case PROJ_AIT:
-      Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
-      if (Z2 < 0) return (FALSE);
-      Z = sqrt (Z2);
-      alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0);
-      delta = DEG_RAD * asin (RAD_DEG*M*Z);
-      break;
-      
+        Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
+        if (Z2 < 0) return (FALSE);
+        Z = sqrt (Z2);
+        alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0);
+        delta = DEG_RAD * asin (RAD_DEG*M*Z);
+        break;
       case PROJ_GLS:
         /* L,M in degrees, alpha,delta in degrees */
@@ -283 +282 @@
         delta = 3 * DEG_RAD * asin (M/180.0);
         break;
+      case PROJ_MOL:
+        Z = sqrt(2.0 - SQ(RAD_DEG*M));
+        alpha = M_PI*L/(2*Z);
+        double T1 = asin(RAD_DEG*M*M_SQRT1_2)/90.0 + M*Z/180.0;
+        delta = DEG_RAD*asin(T1);
       default:
         return (FALSE);
@@ -448 +452 @@
         *M = 180.0 * sin (RAD_DEG*theta/3.0);
         return (TRUE);
+      case PROJ_MOL: 
+        phi = ra - coords[0].crval1;
+        theta = dec - coords[0].crval2;
+        // given theta, solve for gamma:
+        double So = sin(theta*RAD_DEG);
+        double Go = theta;
+        if (theta > +80) { Go = theta - 10.0;}
+        if (theta < -80) { Go = theta + 10.0;}
+        for (int iter = 0; iter < 5; iter ++) {
+          double Fo = (Go / 90.0) + sin(Go*M_PI/90.0) / M_PI - So;
+          double dFdG = (1 / 90.0) + cos(Go*M_PI/90.0) / 90.0;
+          Go -= Fo / dFdG;
+        }
+        *L = 2*phi * cos(RAD_DEG*Go) * M_SQRT2 / M_PI;
+        *M = 180.0 * sin(RAD_DEG*Go) * M_SQRT2 / M_PI;
+        return (TRUE);
       default:
         return (FALSE);
@@ -1059 +1079 @@
   if (!strncmp(&ctype[4], "-GLS", 4)) return PROJ_GLS;
   if (!strncmp(&ctype[4], "-PAR", 4)) return PROJ_PAR;
+  if (!strncmp(&ctype[4], "-MOL", 4)) return PROJ_MOL;
   return PROJ_NONE;
 }
@@ -1079 +1100 @@
     case PROJ_GLS: strcpy(&ctype[4], "-GLS"); return TRUE;
     case PROJ_PAR: strcpy(&ctype[4], "-PAR"); return TRUE;
+    case PROJ_MOL: strcpy(&ctype[4], "-MOL"); return TRUE;
     case PROJ_NONE: return FALSE;
   }
@@ -1103 +1125 @@
     case PROJ_GLS:
     case PROJ_PAR:
+    case PROJ_MOL:
       return PROJ_MODE_PSEUDOCYL;
     default: