# include "dvoshell.h"

int imbox (int argc, char **argv) {
  
  int j, kapa, Nskip, status, InPic, flipped, N, haveNx, haveNy, Nx, Ny, SOLO_PHU;
  Vector Xvec, Yvec;
  double r, d, x[4], y[4], Rmin, Rmax, Rmid;
  Header header;
  Coords coords;
  Coords mosaic;
  Graphdata graphmode;
  FILE *f;

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

  SOLO_PHU = FALSE;
  if ((N = get_argument (argc, argv, "-phu"))) {
    SOLO_PHU = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc != 2) {
    gprint (GP_ERR, "USAGE: imbox (filename)\n");
    return (FALSE);
  }

  f = fopen (argv[1], "r");
  if (f == NULL) {
    gprint (GP_ERR, "file not found\n");
    return (FALSE);
  }
  
  Rmin = graphmode.coords.crval1 - 182.0;
  Rmax = graphmode.coords.crval1 + 182.0;
  Rmid = 0.5*(Rmin + Rmax);
  
  /* project this image to screen display coords */
  SetVector (&Xvec, OPIHI_FLT, 8);
  SetVector (&Yvec, OPIHI_FLT, 8);

  while (gfits_fread_header (f, &header)) {
    if (!GetCoords (&coords, &header)) goto skip;
    if (!SOLO_PHU && !strcmp (&coords.ctype[4], "-DIS")) {
      mosaic = coords;
      RegisterMosaic (&mosaic);
      goto skip;
    }

    // XXX currently, image uses an unsigned short for NX,XY. this is rather restrictive
    // and needs to be at least checked.
    haveNx = gfits_scan (&header, "IMNAXIS1",   "%d", 1, &Nx);
    haveNy = gfits_scan (&header, "IMNAXIS2",   "%d", 1, &Ny);

    if (!haveNx && !haveNy) {
	haveNx = gfits_scan (&header, "ZNAXIS1",   "%d", 1, &Nx);
	haveNy = gfits_scan (&header, "ZNAXIS2",   "%d", 1, &Ny);
    }

    if (!haveNx && !haveNy) {
	haveNx = gfits_scan (&header, "NAXIS1",   "%d", 1, &Nx);
	haveNy = gfits_scan (&header, "NAXIS2",   "%d", 1, &Ny);
    }

    if (!haveNx || !haveNy) {
	fprintf (stderr, "missing image dimensions in header\n");
	goto skip;
    }

    x[0] = 0;  y[0] = 0;
    x[1] = Nx; y[1] = 0;
    x[2] = Nx; y[2] = Ny;
    x[3] = 0;  y[3] = Ny;
    status = FALSE;
    flipped = FALSE;
    for (j = 0; j < 4; j++) {
      XY_to_RD (&r, &d, x[j], y[j], &coords);
      r = ohana_normalize_angle (r);
      while ((j == 0) && (r < Rmin)) { flipped = TRUE; r += 360.0; }
      while ((j == 0) && (r > Rmax)) { flipped = TRUE; r -= 360.0; }
      if ((j > 0) && flipped) {
	while (r < Rmid) r+= 360.0;
	while (r > Rmid) r-= 360.0;
      }
      status |= RD_to_XY (&Xvec.elements.Flt[2*j], &Yvec.elements.Flt[2*j], r, d, &graphmode.coords);
      if (j > 0) {
	Xvec.elements.Flt[2*j - 1] = Xvec.elements.Flt[2*j];
	Yvec.elements.Flt[2*j - 1] = Yvec.elements.Flt[2*j];
      }
    }
    Xvec.elements.Flt[7] = Xvec.elements.Flt[0];
    Yvec.elements.Flt[7] = Yvec.elements.Flt[0];
    InPic = FALSE;
    for (j = 0; j < 8; j+=2) {
      if ((Xvec.elements.Flt[j] >= graphmode.xmin) && 
	  (Xvec.elements.Flt[j] <= graphmode.xmax) && 
	  (Yvec.elements.Flt[j] >= graphmode.ymin) && 
	  (Yvec.elements.Flt[j] <= graphmode.ymax))
	InPic = TRUE;
    }

    Xvec.Nelements = Yvec.Nelements = 8;
    if (InPic) {
      graphmode.style = 2; /* points */
      graphmode.ptype = 100; /* connect pairs of points */
      graphmode.etype = 0;
      PlotVectorPair (kapa, &Xvec, &Yvec, &graphmode);
    }
  skip:
    Nskip = gfits_data_size (&header);
    fseek (f, Nskip, SEEK_CUR); 
    gfits_free_header (&header);
  }
  fclose (f);
  free (Xvec.elements.Flt);
  free (Yvec.elements.Flt);
  return (TRUE);

}


