# include "dvoshell.h"

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

  int i, N, Nimage, status, TimeSelect, ByName, xs, ys;
  time_t tzero, tend;
  double pixscale, dPix, r, d, Xi, Yi, Xs, Ys, x[2], y[2], trange;
  Image *image;
  char name[256];
  float *V;
  int Nx, Ny;
  Buffer *buf;
  Coords coords;

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

  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 != 4) {
    gprint (GP_ERR, "USAGE: skycoverage (buffer) (pixscale) (dPix) [-time start range] [-name name]\n");
    gprint (GP_ERR, "       (buffer) saves bitmapped AIT plot\n");
    gprint (GP_ERR, "       (pixscale) specifies the pixel size in degrees\n");
    gprint (GP_ERR, "       note: we need 64800 / (pixscale)^2 pixels to represent the sky\n");
    return (FALSE);
  }
  
  if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
  pixscale = atof(argv[2]);
  dPix     = atof(argv[3]);

  Nx = 360/pixscale;
  Ny = 180/pixscale;

  gfits_free_matrix (&buf[0].matrix);
  gfits_free_header (&buf[0].header);
  CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0);
  strcpy (buf[0].file, "(empty)");

  coords.crval1 = 180;
  coords.crval2 = 0;
  coords.crpix1 = 0.5*Nx;
  coords.crpix2 = 0.5*Ny;
  strcpy (coords.ctype, "DEC--AIT");
  coords.pc1_1 = coords.pc2_2 = 1;
  coords.pc1_2 = coords.pc2_1 = 0;
  coords.cdelt1 = coords.cdelt2 = pixscale;
  coords.Npolyterms = 0;

  PutCoords (&coords, &buf[0].header);

  if ((image = LoadImages (&Nimage)) == NULL) return (FALSE);
  BuildChipMatch (image, Nimage);

  V = (float *)buf[0].matrix.buffer;
  bzero (V, Nx*Ny*sizeof(float));

  for (ys = 0; ys < Ny; ys++) {
    for (xs = 0; xs < Nx; xs++) {
      status = XY_to_RD (&r, &d, (double)(xs), (double)(ys), &coords);
      status &= (r > 0);
      status &= (r < 360);
      if (status) {
	V[ys*Nx + xs] = 2;
      }
    }
  }

  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;
    } else {
      x[0] = 0;                y[0] = 0;
      x[1] = image[i].NX;      y[1] = image[i].NY;
    }
    status = FALSE;
    
    for (Yi = y[0] + 0.5*dPix; Yi < y[1]; Yi += dPix) {
      for (Xi = x[0] + 0.5*dPix; Xi < x[1]; Xi += dPix) {
	XY_to_RD (&r, &d, Xi, Yi, &image[i].coords);
	while (r <   0.0) { r += 360.0; }
	while (r > 360.0) { r -= 360.0; }
	status = RD_to_XY (&Xs, &Ys, r, d, &coords);
	if (status) {
	  xs = (int)Xs;
	  ys = (int)Ys;
	  V[ys*Nx + xs] = 1;
	}
      }
    }
  }

  free (image);
  return (TRUE);
}


