# include "dvoshell.h"

int wordhash (char *word);

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

  int WITH_MOSAIC, SOLO_MOSAIC;
  off_t i, Nimage;
  int N, status, TimeSelect, ByName, xs, ys;
  time_t tzero, tend;
  double pixscale, dX, dY, Npts, r, d, Xi, Yi, Xs, Ys, x[2], y[2], trange, RaCenter, DecCenter;
  Image *image;
  char name[256], projection[16];
  float *V;
  int Nx, Ny;
  Buffer *buf;
  Coords coords;
  int typehash;
  int PhotcodeSelect;
  PhotCode *PhotcodeValue;

  WITH_MOSAIC = FALSE;
  if ((N = get_argument (argc, argv, "+mosaic"))) {
    remove_argument (N, &argc, argv);
    WITH_MOSAIC = TRUE;
  }

  SOLO_MOSAIC = FALSE;
  if ((N = get_argument (argc, argv, "-mosaic"))) {
    remove_argument (N, &argc, argv);
    SOLO_MOSAIC = TRUE;
    WITH_MOSAIC = TRUE;
  }

  RaCenter = 0.0;
  DecCenter = 0.0;
  if ((N = get_argument (argc, argv, "-center"))) {
    remove_argument (N, &argc, argv);
    RaCenter = atof (argv[N]);
    remove_argument (N, &argc, argv);
    DecCenter = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  Nx = 0;
  Ny = 0;
  if ((N = get_argument (argc, argv, "-size"))) {
    remove_argument (N, &argc, argv);
    Nx = atof (argv[N]);
    remove_argument (N, &argc, argv);
    Ny = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  pixscale = 1.0;
  if ((N = get_argument (argc, argv, "-scale"))) {
    remove_argument (N, &argc, argv);
    pixscale = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

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

  strcpy (projection, "DEC--AIT");
  if ((N = get_argument (argc, argv, "-proj"))) {
    remove_argument (N, &argc, argv);
    if (!strcasecmp(argv[N], "TAN")) {
	strcpy (projection, "DEC--TAN");
    }	
    if (!strcasecmp(argv[N], "SIN")) {
	strcpy (projection, "DEC--SIN");
    }	
    if (!strcasecmp(argv[N], "GLS")) {
	strcpy (projection, "DEC--GLS");
    }	
    if (!strcasecmp(argv[N], "PAR")) {
	strcpy (projection, "DEC--PAR");
    }	
    remove_argument (N, &argc, argv);
  }

  PhotcodeValue = NULL;
  PhotcodeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-photcode"))) {
    if (!InitPhotcodes ()) return (FALSE);
    PhotcodeSelect = TRUE;
    remove_argument (N, &argc, argv);
    PhotcodeValue = GetPhotcodebyName (argv[N]);
    if (PhotcodeValue == NULL) {
      gprint (GP_ERR, "photcode not found in photcode table\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
  }

  TimeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-time"))) {
    remove_argument (N, &argc, argv);
    if (!ohana_str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!ohana_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 (!ohana_str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!ohana_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 != 3) {
    gprint (GP_ERR, "USAGE: skycoverage (buffer) (Npts)\n");
    gprint (GP_ERR, "  options: [-scale pixscale] [-center ra dec] [-size Nx Nx] [-proj projection] [-time start range] [-trange start stop] [-name name] [-photcode name] [+mosaic] [-mosaic]\n");
    gprint (GP_ERR, "       (buffer) saves bitmapped image\n");
    gprint (GP_ERR, "       (Npts) gives the number of test points per image in each dimension\n");
    gprint (GP_ERR, "       -scale (pixscale)  : specifies the pixel size in degrees [1.0]\n");
    gprint (GP_ERR, "       -center (ra) (dec) : specifies the center of the field [0.0, 0.0]\n");
    gprint (GP_ERR, "       -size (Nx) (Ny)    : specifies the size of the image [360/scale, 180/scale]\n");
    gprint (GP_ERR, "       -proj (projection) : specifies the projection choice [AIT]\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);
  Npts = atof(argv[2]);

  if (!Nx || !Ny) {
    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 = RaCenter;
  coords.crval2 = DecCenter;
  coords.crpix1 = 0.5*Nx;
  coords.crpix2 = 0.5*Ny;
  strcpy (coords.ctype, projection);
  coords.pc1_1 = -1.0;
  coords.pc2_2 = +1.0;
  coords.pc1_2 =  0.0;
  coords.pc2_1 =  0.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;
      }
    }
  }

  int DistortImage = wordhash ("-DIS");

  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 (PhotcodeSelect) {
      if (PhotcodeValue[0].type == PHOT_DEP) {
	if (PhotcodeValue[0].code != image[i].photcode) continue;
      } else {
	if (PhotcodeValue[0].code != GetPhotcodeEquivCodebyCode (image[i].photcode)) continue;
      }
    }

    if (!FindMosaicForImage (image, Nimage, i)) continue;

    typehash = wordhash (&image[i].coords.ctype[4]);

    /* DIS images represent a field, not a chip */
    if ((typehash == DistortImage) && !WITH_MOSAIC) continue;
    if ((typehash != DistortImage) &&  SOLO_MOSAIC) 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;
    
    dX = (x[1] - x[0]) / Npts;
    dY = (y[1] - y[0]) / Npts;
    
    for (Yi = y[0] + 0.5*dY; Yi < y[1]; Yi += dY) {
      for (Xi = x[0] + 0.5*dX; Xi < x[1]; Xi += dX) {
	XY_to_RD (&r, &d, Xi, Yi, &image[i].coords);
	r = ohana_normalize_angle (r);
	if (r - RaCenter > +180.0) r -= 360.0;
	if (r - RaCenter < -180.0) r += 360.0;
	status = RD_to_XY (&Xs, &Ys, r, d, &coords);
	if (Xs < 0) continue;
	if (Ys < 0) continue;
	if (Xs >= Nx) continue;
	if (Ys >= Ny) continue;
	if (status) {
	  xs = (int)Xs;
	  ys = (int)Ys;
	  V[ys*Nx + xs] = 1;
	}
      }
    }
  }

  free (image);
  return (TRUE);
}


