Index: /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/include/libdvo_astro.h
===================================================================
--- /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/include/libdvo_astro.h	(revision 41965)
+++ /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/include/libdvo_astro.h	(revision 41966)
@@ -27,4 +27,5 @@
   PROJ_PAR, // pseudocyl
   PROJ_ZPN, // zenithal
+  PROJ_MOL, // pseudocyl 
 } OhanaProjection;
 
Index: /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/src/coordops.c
===================================================================
--- /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/src/coordops.c	(revision 41965)
+++ /branches/eam_branches/ipp-20211108/Ohana/src/libdvo/src/coordops.c	(revision 41966)
@@ -266,11 +266,10 @@
     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,4 +282,9 @@
 	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,4 +452,20 @@
 	*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,4 +1079,5 @@
   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,4 +1100,5 @@
     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,4 +1125,5 @@
     case PROJ_GLS:
     case PROJ_PAR:
+    case PROJ_MOL:
       return PROJ_MODE_PSEUDOCYL;
     default: