# include "dvoshell.h"
# define BETA 0.41421

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

  int i, j, Nimage, status, InPic, leftside, *plist, TimeSelect, ByName;
  time_t tzero, tend;
  int N, NPTS, n, npts, Npts;
  Vector Xvec, Yvec;
  double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius;
  Image *image;
  Graphdata graphmode;
  char name[256];

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

  ByName = FALSE;
  if ((N = get_argument (argc, argv, "-name"))) {
    remove_argument (N, &argc, argv);
    strcpy (name, argv[N]);
    remove_argument (N, &argc, argv);
    ByName = TRUE;
  }

  Radius = 45;
  if ((N = get_argument (argc, argv, "-radius"))) {
    remove_argument (N, &argc, argv);
    Radius = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  TimeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-time"))) {
    remove_argument (N, &argc, argv);
    if (!str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!str_to_dtime (argv[N], &trange)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (trange < 0) {
      trange = fabs (trange);
      tzero -= trange;
    }
    TimeSelect = TRUE;
  }
  if ((N = get_argument (argc, argv, "-trange"))) {
    remove_argument (N, &argc, argv);
    if (!str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!str_to_time (argv[N], &tend)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    trange = tend - tzero;
    if (trange < 0) {
      trange = fabs (trange);
      tzero -= trange;
    }
    TimeSelect = TRUE;
  }
 
  if (argc != 1) {
    gprint (GP_ERR, "USAGE: image [-time start range]\n");
    return (FALSE);
  }
  SetGraph (graphmode);
  
  /* it is not an error for the database not to have any images */
  if ((image = LoadImages (&Nimage)) == NULL) return (TRUE);
  BuildChipMatch (image, Nimage);

  Rmin = graphmode.coords.crval1 - 180.0;
  Rmax = graphmode.coords.crval1 + 180.0;
  Rmid = 0.5*(Rmin + Rmax);
  
  npts = NPTS = 200;
  ALLOCATE (Xvec.elements, float, NPTS);
  ALLOCATE (Yvec.elements, float, NPTS);
  ALLOCATE (plist, int, NPTS);
  n = N = 0;
  for (i = 0; i < Nimage; i++) {
    if (ByName && strcmp (name, image[i].name)) continue;
    if (TimeSelect && ((image[i].tzero < tzero) || (image[i].tzero+image[i].trate*image[i].NY > tzero + trange))) continue;
    if (!FindMosaicForImage (image, Nimage, i)) continue;

    /* project this image to screen display coords */
    /* DIS images represent a field, not a chip */
    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) {
      x[0] = -0.5*image[i].NX; y[0] = -0.5*image[i].NY;
      x[1] = +0.5*image[i].NX; y[1] = -0.5*image[i].NY;
      x[2] = +0.5*image[i].NX; y[2] = +0.5*image[i].NY;
      x[3] = -0.5*image[i].NX; y[3] = +0.5*image[i].NY;
    } else {
      x[0] = 0;           y[0] = 0;
      x[1] = image[i].NX; y[1] = 0;
      x[2] = image[i].NX; y[2] = image[i].NY;
      x[3] = 0;           y[3] = image[i].NY;
    }

    Npts = 4;
    status = TRUE;
    leftside = FALSE;
    for (j = 0; j < 4; j++) {
      status = XY_to_RD (&r[j], &d[j], x[j], y[j], &image[i].coords);
      if (!status) break;
      while (r[j] < Rmin) { r[j] += 360.0; }
      while (r[j] > Rmax) { r[j] -= 360.0; }
      if (j == 0) {
	leftside = (r[j] < Rmid);
      } 
      if (j > 0) { 
	if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; }
	if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; }
      }
    }

    // extremely large-scale images with certain projections will have odd boundaries.
    // eg, the ASCA images are essentially circles of radius ~60 degrees.  plot these as 
    // octagons with some dummy size.
    if (!status) {
      int jp, xo, yo;
      double rc, dc;
      xo = 0.5*image[i].NX;
      yo = 0.5*image[i].NY;
      // is the image center on the screen?
      status = XY_to_RD (&rc, &dc, 0.5*x[2], 0.5*y[2], &image[i].coords);
      if (status && (image[i].NX * image[i].coords.cdelt2 > 90)) {
	// draw an octagon with radius 45 degrees
	double dX = Radius / image[i].coords.cdelt2;
	x[0] = xo+dX;      y[0] = yo+BETA*dX;
	x[1] = xo+BETA*dX; y[1] = yo+dX;
	x[2] = xo-BETA*dX; y[2] = yo+dX;
	x[3] = xo-dX;      y[3] = yo+BETA*dX;
	x[4] = xo-dX;      y[4] = yo-BETA*dX;
	x[5] = xo-BETA*dX; y[5] = yo-dX;
	x[6] = xo+BETA*dX; y[6] = yo-dX;
	x[7] = xo+dX;      y[7] = yo-BETA*dX;
	Npts = 8;
	j = 0;
	for (jp = 0; jp < 8; jp++) {
	  status = XY_to_RD (&r[j], &d[j], x[jp], y[jp], &image[i].coords);
	  if (!status) continue;
	  while (r[j] < Rmin) { r[j] += 360.0; }
	  while (r[j] > Rmax) { r[j] -= 360.0; }
	  if (j == 0) {
	    leftside = (r[j] < Rmid);
	  } 
	  if (j > 0) { 
	    if ( leftside && (r[j] > Rmid + 90)) { r[j] -= 360.0; }
	    if (!leftside && (r[j] < Rmid - 90)) { r[j] += 360.0; }
	  }
	  j++;
	}
	Npts = j;
      } else {
	continue;
      }
    }
    if (Npts == 0) continue;

    // if any of the points is on the screen, plot the image
    status = FALSE;
    for (j = 0; j < Npts; j++) {
      status |= fRD_to_XY (&Xvec.elements[N+2*j], &Yvec.elements[N+2*j], r[j], d[j], &graphmode.coords);
      if (j > 0) {
	Xvec.elements[N+2*j - 1] = Xvec.elements[N+2*j];
	Yvec.elements[N+2*j - 1] = Yvec.elements[N+2*j];
      }
    }
    Xvec.elements[N+2*Npts-1] = Xvec.elements[N];
    Yvec.elements[N+2*Npts-1] = Yvec.elements[N];
    InPic = FALSE;
    for (j = 0; j < 2*Npts; j+=2) {
      if ((Xvec.elements[N+j] >= graphmode.xmin) && 
	  (Xvec.elements[N+j] <= graphmode.xmax) && 
	  (Yvec.elements[N+j] >= graphmode.ymin) && 
	  (Yvec.elements[N+j] <= graphmode.ymax))
	InPic = TRUE;
    }
    if (InPic && status) {
      plist[n] = i;
      n++;
      if (n > npts - 1) {
	npts += 200;
	REALLOCATE (plist, int, npts);
      }
      N+=2*Npts;
      if (N + 16 >= NPTS) {  /* need to leave room for 8 point image */
	NPTS += 400;
	REALLOCATE (Xvec.elements, float, NPTS);
	REALLOCATE (Yvec.elements, float, NPTS);
      }
    }
  }

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

  free (Xvec.elements);
  free (Yvec.elements);
  free (image);
  return (TRUE);

}


