# include "dvoshell.h"

int RD_to_XYpic (double *x, double *y, double r, double d, Coords *coords, double Rmin, double Rmax, double Rmid, int *leftside);

int skycat (int argc, char **argv) {
  
  double Radius;
  int i, j, N, Nregions, kapa, ShowAll, NPTS, Npts, leftside, Depth, VERBOSE;
  struct stat filestat;
  Vector Xvec, Yvec;
  Graphdata graphmode;
  double X[4], Y[4], Rmin, Rmax, Rmid;
  SkyTable *sky;
  SkyList *skylist;
  SkyRegion **regions;

  VERBOSE = FALSE;
  if ((N = get_argument (argc, argv, "-v"))) {
    remove_argument (N, &argc, argv);
    VERBOSE = TRUE;
  }
  ShowAll = FALSE;
  if ((N = get_argument (argc, argv, "-all"))) {
    remove_argument (N, &argc, argv);
    ShowAll = TRUE;
  }
  Depth = -1;
  if ((N = get_argument (argc, argv, "-depth"))) {
    remove_argument (N, &argc, argv);
    Depth = atoi (argv[N]);
    remove_argument (N, &argc, argv);    
  }

  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;

  if (argc != 1) {
    gprint (GP_ERR, "USAGE: skycat [-all]\n");
    return (FALSE);
  }

  Radius = MAX (fabs(graphmode.xmax), fabs(graphmode.ymax));

  sky = GetSkyTable ();
  skylist = SkyListByRadius (sky, Depth, graphmode.coords.crval1, graphmode.coords.crval2, Radius);
  
  if (VERBOSE) gprint (GP_ERR, "region: %6.2f - %6.2f, %6.2f - %6.2f\n", 
			graphmode.coords.crval1 - Radius, graphmode.coords.crval1 + Radius, 
			graphmode.coords.crval2 - Radius, graphmode.coords.crval2 + Radius);

  Rmin = graphmode.coords.crval1 - 180.0;
  Rmax = graphmode.coords.crval1 + 180.0;
  Rmid = 0.5*(Rmin + Rmax);

  NPTS = 200;
  ALLOCATE (Xvec.elements, float, NPTS);
  ALLOCATE (Yvec.elements, float, NPTS);
  Npts = 0;
   
  regions = skylist[0].regions;
  Nregions = skylist[0].Nregions;

  for (i = 0; i < Nregions; i++) {
    if (ShowAll || (stat (skylist[0].filename[i], &filestat) != -1)) {
      if (VERBOSE) gprint (GP_ERR, "%3d %s %6.2f - %6.2f, %6.2f - %6.2f\n", i, regions[i][0].name, 
			    regions[i][0].Rmin, regions[i][0].Rmax, regions[i][0].Dmin, regions[i][0].Dmax);

      leftside = -1;
      RD_to_XYpic (&X[0], &Y[0], regions[i][0].Rmin, regions[i][0].Dmin, &graphmode.coords, Rmin, Rmax, Rmid, &leftside);
      RD_to_XYpic (&X[1], &Y[1], regions[i][0].Rmin, regions[i][0].Dmax, &graphmode.coords, Rmin, Rmax, Rmid, &leftside);
      RD_to_XYpic (&X[2], &Y[2], regions[i][0].Rmax, regions[i][0].Dmax, &graphmode.coords, Rmin, Rmax, Rmid, &leftside);
      RD_to_XYpic (&X[3], &Y[3], regions[i][0].Rmax, regions[i][0].Dmin, &graphmode.coords, Rmin, Rmax, Rmid, &leftside);

      Xvec.elements[Npts] = X[0];
      Yvec.elements[Npts] = Y[0];
      for (j = 1; j < 4; j++) {
	Xvec.elements[Npts + j*2 - 0] = X[j];
	Yvec.elements[Npts + j*2 - 0] = Y[j];
	Xvec.elements[Npts + j*2 - 1] = X[j];
	Yvec.elements[Npts + j*2 - 1] = Y[j];
      }
      Xvec.elements[Npts+7] = Xvec.elements[Npts];
      Yvec.elements[Npts+7] = Yvec.elements[Npts];
      Npts += 8;
      if (Npts > NPTS - 1) {  /* this is OK because NPTS is made always a multiple of 8 */
	NPTS += 200;
	REALLOCATE (Xvec.elements, float, NPTS);
	REALLOCATE (Yvec.elements, float, NPTS);
      }
    }
  }

  gprint (GP_ERR, "plotting %d catalogs\n", Npts/8);
  Xvec.Nelements = Xvec.Nelements = Npts;
  if (Npts > 0) {
    graphmode.style = 2; /* points */
    graphmode.ptype = 100; /* connect pairs of points */
    graphmode.etype = 0;
    PlotVectorPair (kapa, Npts, Xvec.elements, Yvec.elements, &graphmode);
  }

  free (Xvec.elements);
  free (Yvec.elements);
  free (regions);

  return (TRUE);

}


int RD_to_XYpic (double *x, double *y, double r, double d, Coords *coords, double Rmin, double Rmax, double Rmid, int *leftside) {

  while (r < Rmin) { r += 360.0; }
  while (r > Rmax) { r -= 360.0; }

  if (*leftside == -1) {
    *leftside = (r < Rmid);
  } else {
    if (  *leftside && (r > Rmid + 90)) { r -= 360.0; }
    if (! *leftside && (r < Rmid - 90)) { r += 360.0; }
  }

  RD_to_XY (x, y, r, d, coords);

  return (TRUE);
}
