# include "addstar.h"

Image *gimages (FILE *f, Image *image, int *Npimage) {
  
  int i, j, k, Nimage, addtolist, size;
  int NTIMAGE, Ntimage, ntimage;
  int NPIMAGE, npimage;
  Image *timage, *pimage;
  Coords tcoords;
  Header header;
  double r, d;
  double Xi[5], Yi[5], Xo[5], Yo[5];  /* image and original corners */
  double zeropt;
  struct stat filestatus;
  
  fseek (f, 0, SEEK_SET);

  /* read header */
  if (!fits_load_header (f, &header)) {
    if (VERBOSE) fprintf (stderr, "can't read image catalog %s\n", ImageCat);
    Shutdown ();
  }

  fits_scan (&header, "ZERO_PT", "%lf", 1, &zeropt);
  if (fabs (ZeroPt - zeropt) > 1e-4) {
    fprintf (stderr, "zero point in image (%f:%s) inconsistent with zero point in catalog (%f)\n",
	     ZeroPt, image[0].name, zeropt);
    Shutdown ();
  }

  /* project onto rectilinear grid with 1 arcsec pixels */
  /* we keep the original crpix1,2 and crref1,2 */
  tcoords = image[0].coords;
  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
  strcpy (tcoords.ctype, "RA---TAN");

  /* define original corners */
  Xo[0] = 0;           Yo[0] = 0;
  Xo[1] = image[0].NX; Yo[1] = 0;
  Xo[2] = image[0].NX; Yo[2] = image[0].NY;
  Xo[3] = 0;           Yo[3] = image[0].NY;
  Xo[4] = 0;           Yo[4] = 0;  /* so we can make a loop easily */
  for (j = 0; j < 5; j++) {
    XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords);
    RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords);
  }
  
  /* advance to data segment */ 
  fseek (f, header.size, SEEK_SET); 

  Nimage = 0;
  fits_scan (&header, "NIMAGES", "%d", 1, &Nimage);
  /* check that file size makes sense */
  if (stat (ImageCat, &filestatus) == -1) {
    fprintf (stderr, "failed to get status of image catalog\n");
    Shutdown ();
  }
  size = Nimage*sizeof(Image) + header.size;
  if (size != filestatus.st_size) {
    int Ndata;

    Ndata = (filestatus.st_size - header.size) / sizeof (Image);
    fprintf (stderr, "image catalog has inconsistent size\n");
    fprintf (stderr, "header: %d images, data: %d images\n", Nimage, Ndata);
    if (!FORCE) Shutdown ();
    
    Nimage = MIN (Nimage, Ndata);
    fprintf (stderr, "using %d images\n", Nimage);
  } 

  NTIMAGE = 100;
  ALLOCATE (timage, Image, NTIMAGE);

  npimage = 0;
  NPIMAGE = 20;
  ALLOCATE (pimage, Image, NPIMAGE);

  for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) {
    ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image");
    for (i = 0; i < ntimage; i++) {
      addtolist = FALSE;
      /* define image corners */
      Xi[0] = 0;            Yi[0] = 0;
      Xi[1] = timage[i].NX; Yi[1] = 0;
      Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
      Xi[3] = 0;            Yi[3] = timage[i].NY;
      Xi[4] = 0;            Yi[4] = 0;  /* so we can make a loop easily */
      /* transform to tcoords */
      for (j = 0; j < 5; j++) {
	XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords);
	RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords);
      }
      /* check if any edges cross */
      for (j = 0; (j < 4) && !addtolist; j++) {
	for (k = 0; (k < 4) && !addtolist; k++) {
	  addtolist |= edge_check (&Xi[j], &Yi[j], &Xo[k], &Yo[k]);
	}
      }
      if (addtolist) {
	pimage[npimage] = timage[i];
	npimage ++;
	if (npimage == NPIMAGE) {
	  NPIMAGE += 20;
	  REALLOCATE (pimage, Image, NPIMAGE);
	}
      }
    }

    if (ntimage == 0) {
      fprintf (stderr, "ERROR: expected %d images, only found %d\n", Nimage, Ntimage);
      Shutdown ();
    }
  }

  if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage);

  *Npimage = npimage;
  return (pimage);

}
  
int edge_check (x1, y1, x2, y2)
double *x1, *y1, *x2, *y2;
{

  double theta1, theta2;
  double Theta1, Theta2;

  theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); 
  theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); 

  if (theta1*theta2 < 0.0) {
    return (FALSE);
  }

  if (fabs(theta1) < fabs(theta2)) {
    return (FALSE);
  }

  Theta1 = theta1;
  Theta2 = theta2;
  theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); 
  theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); 
  
 
  if (theta1*theta2 < 0.0) {
    return (FALSE);
  }

  if (fabs(theta1) < fabs(theta2)) {
    return (FALSE);
  }

  return (TRUE);

}

double opening_angle (x1, y1, x2, y2, x3, y3)
double x1, y1, x2, y2, x3, y3;
{

  double dx1, dy1, dx2, dy2, ct, st, theta;

  dx1 = x1 - x2;
  dy1 = y1 - y2;
  
  dx2 = x3 - x2;
  dy2 = y3 - y2;
  
  ct = (dx1*dx2 + dy1*dy2);
  st = (dx1*dy2 - dx2*dy1);

  theta = atan2 (st, ct);

  return (theta);

}

