# include "astro.h"

int csystem (int argc, char **argv) {

  /* USAGE: csystem [C/G/E/H] [C/G/E/H] [epoch] */
  int i;
  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");
    return (FALSE);
  }

  switch (argv[1][0]) {
    case 'C': input = COORD_CELESTIAL; break;
    case 'G': input = COORD_GALACTIC; break;
    case 'E': input = COORD_ECLIPTIC; break;
    default: abort();
  }

  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);
      
    ApplyTransform (&x, &y, X, Y, transform);

    gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
    return (TRUE);
  }

  /* find vectors */
  if ((xvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((yvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);

  if (xvec[0].Nelements != yvec[0].Nelements) {
    gprint (GP_ERR, "vectors %s and %s not the same length\n", argv[3], argv[4]);
    return (FALSE);
  }
  
  xptr = xvec[0].elements;
  yptr = yvec[0].elements;

  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;
  }

  return (TRUE);

}
