Index: trunk/Ohana/src/opihi/cmd.astro/csystem.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/csystem.c	(revision 8427)
+++ 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;
   }