# 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, Xo, xo, phi, T;
  double sin_x, sin_y, cos_x, cos_y;
  float *x, *y;
  struct timeval now;
  struct tm *local;
  Vector *xvec, *yvec;
   
  if (argc != 5) {
    gprint (GP_ERR, "USAGE: csystems [C/G/E/H] [C/G/E/H] X Y\n");
    return (FALSE);
  }

  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;
    }
  }
 
  Xo *= RAD_DEG;

  if (SelectScalar (argv[3], &X)) {
      if (!SelectScalar (argv[4], &Y)) return (FALSE);
      
      X *= RAD_DEG;
      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);

      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);
  }
  
  x = xvec[0].elements;
  y = 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;
  }

  return (TRUE);

}
