Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 16934)
+++ /trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 16935)
@@ -22,4 +22,5 @@
 $(SRC)/cplot.$(ARCH).o		   \
 $(SRC)/csystem.$(ARCH).o	   \
+$(SRC)/coord_systems.$(ARCH).o	   \
 $(SRC)/ctimes.$(ARCH).o	   \
 $(SRC)/cval.$(ARCH).o		   \
Index: /trunk/Ohana/src/opihi/cmd.astro/coord_systems.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/coord_systems.c	(revision 16935)
+++ /trunk/Ohana/src/opihi/cmd.astro/coord_systems.c	(revision 16935)
@@ -0,0 +1,135 @@
+# include "astro.h"
+    
+CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output) {
+
+  CoordTransform *transform;
+  struct timeval now;
+  struct tm *local;
+  double T;
+
+  ALLOCATE (transform, CoordTransform, 1);
+  transform->isIdentity = FALSE;
+  
+  switch (input) {
+    case COORD_CELESTIAL:
+      switch (output) {
+	case COORD_CELESTIAL:
+	  transform->phi = 0.0;
+	  transform->Xo  = 0.0;
+	  transform->xo  = 0.0;
+	  transform->isIdentity = TRUE;
+	  break;
+	case COORD_GALACTIC:
+	  transform->phi = -62.60*RAD_DEG;
+	  transform->Xo  = 282.25*RAD_DEG;
+	  transform->xo  =  33.00;
+	  break;
+	case COORD_ECLIPTIC:
+	  gettimeofday (&now, NULL);
+	  local = localtime (&now.tv_sec);
+	  T = local[0].tm_year / 100.0;
+	  transform->phi = -1.0*RAD_DEG*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
+	  transform->Xo  = 0.0;
+	  transform->xo  = 0.0;
+	  break;
+	default:
+	  abort();
+      }
+      break;
+    case COORD_ECLIPTIC:
+      switch (output) {
+	case COORD_CELESTIAL:
+	  gettimeofday (&now, (struct timezone *) NULL);
+	  local = localtime (&now.tv_sec);
+	  T = local[0].tm_year / 100.0;
+	  transform->phi = RAD_DEG*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
+	  transform->Xo = 0.0;
+	  transform->xo = 0.0;
+	  break;
+	case COORD_GALACTIC:
+	  return NULL;
+	  break;
+	case COORD_ECLIPTIC:
+	  transform->phi = 0.0;
+	  transform->Xo  = 0.0;
+	  transform->xo  = 0.0;
+	  transform->isIdentity = TRUE;
+	  break;
+      }
+    case COORD_GALACTIC:
+      switch (output) {
+	case COORD_CELESTIAL:
+	  transform->phi =  62.60*RAD_DEG;
+	  transform->Xo  =  33.00*RAD_DEG;
+	  transform->xo  = 282.25;
+	  break;
+	case COORD_GALACTIC:
+	  transform->phi = 0.0;
+	  transform->Xo  = 0.0;
+	  transform->xo  = 0.0;
+	  transform->isIdentity = TRUE;
+	  break;
+	case COORD_ECLIPTIC:
+	  return NULL;
+      }
+  }
+ 
+  // pre-calculated constants:
+  transform->sin_phi_cos_Xo = sin(transform->phi)*cos(transform->Xo);
+  transform->sin_phi_sin_Xo = sin(transform->phi)*sin(transform->Xo);
+  transform->cos_phi        = cos(transform->phi);
+  
+  transform->cos_phi_cos_Xo = cos(transform->phi)*cos(transform->Xo);
+  transform->cos_phi_sin_Xo = cos(transform->phi)*sin(transform->Xo);
+  transform->sin_phi        = sin(transform->phi);
+  transform->cos_Xo 	    = cos(transform->Xo);
+  transform->sin_Xo 	    = sin(transform->Xo);
+
+  return transform;
+}
+
+// input and output coordinates are in degrees
+int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform) {
+
+  double sin_x, sin_y, cos_x, cos_y;
+
+
+  if (transform == NULL) return (FALSE);
+
+  if (transform->isIdentity) {
+    *x = X;
+    *y = Y;
+    return (TRUE);
+  }
+
+  X *= RAD_DEG;
+  Y *= RAD_DEG;
+
+  // recast with constants extracted:
+  sin_y = cos(Y)*sin(X)*transform->sin_phi_cos_Xo - cos(Y)*cos(X)*transform->sin_phi_sin_Xo + sin(Y)*transform->cos_phi;
+  cos_y = sqrt (1 - sin_y*sin_y);
+
+  sin_x = cos(Y)*sin(X)*transform->cos_phi_cos_Xo - cos(Y)*cos(X)*transform->cos_phi_sin_Xo - sin(Y)*transform->sin_phi;
+  cos_x = cos(Y)*cos(X)*transform->cos_Xo + cos(Y)*sin(X)*transform->sin_Xo;
+      
+  // atan2 returns -pi : +pi
+  *x = DEG_RAD * atan2 (sin_x, cos_x) + transform->xo;
+  if ((*x) < 0.0) (*x) += 360;
+
+  // should be in range -pi/2 : +pi/2
+  *y = DEG_RAD * atan2 (sin_y, cos_y);
+
+  return (TRUE);
+}
+
+// sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);
+
+// sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
+// cos_x = cos(Y)*cos(X - Xo) / cos_y;
+
+// multiplying both sides by cos_y:
+// sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi));
+// cos_x = cos(Y)*cos(X - Xo);
+      
+// sin(a - b) = sin(a)*cos(b) - sin(b)*cos(a);
+// cos(a - b) = sin(a)*sin(b) + cos(a)*cos(b);
Index: /trunk/Ohana/src/opihi/cmd.astro/csystem.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 16934)
+++ /trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 16935)
@@ -5,11 +5,10 @@
   /* USAGE: csystem [C/G/E/H] [C/G/E/H] [epoch] */
   int i;
-  double X, Y, Xo, xo, phi, T;
-  double sin_x, sin_y, cos_x, cos_y;
-  float *x, *y;
-  struct timeval now;
-  struct tm *local;
+  double X, Y, x, y;
+  float *xptr, *yptr;
   Vector *xvec, *yvec;
-   
+  CoordTransformSystem input, output;
+  CoordTransform *transform;
+
   if (argc != 5) {
     gprint (GP_ERR, "USAGE: csystems [C/G/E/H] [C/G/E/H] X Y\n");
@@ -17,92 +16,31 @@
   }
 
-  Xo = xo = phi = 0;
   switch (argv[1][0]) {
-  case 'C':
-    switch (argv[2][0]) {
-    case 'C': 
-      gprint (GP_ERR, "same coordinate system\n");
-      return (TRUE);
-      break;
-    case 'G':
-      phi = -62.6*RAD_DEG;
-      Xo = 282.25;
-      xo = 33;
-      break;
-    case 'E':
-      gettimeofday (&now, (struct timezone *) NULL);
-      local = localtime (&now.tv_sec);
-      T = local[0].tm_year / 100.0;
-      phi = -1*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
-      phi *= RAD_DEG;
-      Xo = xo = 0.0;
-      break;
-    }
-    break;
-  case 'E':
-    switch (argv[2][0]) {
-    case 'C': 
-      gettimeofday (&now, (struct timezone *) NULL);
-      local = localtime (&now.tv_sec);
-      T = local[0].tm_year / 100.0;
-      phi = 23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T;
-      phi *= RAD_DEG;
-      Xo = xo = 0.0;
-      break;
-    case 'G':
-      gprint (GP_ERR, "error: conversions between galactic and ecliptic not implemented\n");
-      return (FALSE);
-      phi = -62.6*RAD_DEG;
-      Xo = 282.25;
-      xo = 33;
-      break;
-    case 'E':
-      phi = Xo = xo = 0.0;
-      break;
-    }
-    break;
-  case 'G':
-    switch (argv[2][0]) {
-    case 'C': 
-      phi = 62.6*RAD_DEG;
-      Xo = 33;
-      xo = 282.25;
-      break;
-    case 'G':
-      phi = Xo = xo = 0.0;
-      break;
-    case 'E':
-      gprint (GP_ERR, "error: conversions between galactic and ecliptic not implemented\n");
-      return (FALSE);
-      gettimeofday (&now, (struct timezone *) NULL);
-      local = localtime (&now.tv_sec);
-      T = local[0].tm_year / 100.0;
-      phi = -1*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
-      Xo = xo = 0.0;
-      break;
-    }
+    case 'C': input = COORD_CELESTIAL; break;
+    case 'G': input = COORD_GALACTIC; break;
+    case 'E': input = COORD_ECLIPTIC; break;
+    default: abort();
   }
- 
-  Xo *= RAD_DEG;
 
+  switch (argv[2][0]) {
+    case 'C': output = COORD_CELESTIAL; break;
+    case 'G': output = COORD_GALACTIC; break;
+    case 'E': output = COORD_ECLIPTIC; break;
+    default: abort();
+  }
+
+  transform = InitTransform (input, output);
+  if (transform == NULL) {
+    gprint (GP_ERR, "transform %c to %c is not yet defined\n", argv[1][0], argv[2][0]);
+    return (FALSE);
+  }
+    
   if (SelectScalar (argv[3], &X)) {
-      if (!SelectScalar (argv[4], &Y)) return (FALSE);
+    if (!SelectScalar (argv[4], &Y)) return (FALSE);
       
-      X *= RAD_DEG;
-      Y *= RAD_DEG;
+    ApplyTransform (&x, &y, X, Y, transform);
 
-      sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);
-      cos_y = sqrt (1 - sin_y*sin_y);
-      sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
-      cos_x = cos(Y)*cos(X - Xo) / cos_y;
-      
-      X = (DEG_RAD * atan2 (sin_x, cos_x) + xo + 360);
-      
-      while (X >= 360.0)
-	  X -= 360;
-      Y = DEG_RAD * atan2 (sin_y, cos_y);
-
-      gprint (GP_LOG, "%10.6f %10.6f\n", X, Y);
-      return (TRUE);
+    gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
+    return (TRUE);
   }
 
@@ -116,23 +54,12 @@
   }
   
-  x = xvec[0].elements;
-  y = yvec[0].elements;
+  xptr = xvec[0].elements;
+  yptr = yvec[0].elements;
 
-  for (i = 0; i < xvec[0].Nelements; i++, x++, y++) {
-    X = *x*RAD_DEG;
-    Y = *y*RAD_DEG;
-
-    sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);
-    cos_y = sqrt (1 - sin_y*sin_y);
-    sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
-    cos_x = cos(Y)*cos(X - Xo) / cos_y;
-    
-    X = (DEG_RAD * atan2 (sin_x, cos_x) + xo + 360);
-    
-    while (X >= 360.0)
-      X -= 360;
-    Y = DEG_RAD * atan2 (sin_y, cos_y);
-    *x = X;
-    *y = Y;
+  for (i = 0; i < xvec[0].Nelements; i++, xptr++, yptr++) {
+    // ApplyTransform takes (double *), but xptr, yptr are (float *)
+    ApplyTransform (&x, &y, *xptr, *yptr, transform);
+    *xptr = x;
+    *yptr = y;
   }
 
Index: /trunk/Ohana/src/opihi/include/astro.h
===================================================================
--- /trunk/Ohana/src/opihi/include/astro.h	(revision 16934)
+++ /trunk/Ohana/src/opihi/include/astro.h	(revision 16935)
@@ -11,3 +11,23 @@
 void InitAstro ();
 
+typedef struct {
+  int    isIdentity;	      // identity transformation 
+  double phi;		      // saved in radians
+  double Xo;		      // saved in radians
+  double xo;		      // saved in degrees
+  double sin_phi_cos_Xo;      // pre-computed values
+  double sin_phi_sin_Xo;      // pre-computed values
+  double cos_phi_cos_Xo;      // pre-computed values
+  double cos_phi_sin_Xo;      // pre-computed values
+  double cos_phi;	      // pre-computed values
+  double sin_phi;	      // pre-computed values
+  double cos_Xo;	      // pre-computed values
+  double sin_Xo;	      // pre-computed values
+} CoordTransform;
+
+typedef enum {COORD_NONE, COORD_CELESTIAL, COORD_GALACTIC, COORD_ECLIPTIC} CoordTransformSystem;
+
+CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output);
+int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform);
+
 # endif