# include "addstar.h"
double opening_angle ();

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

  /* check if image datafile exists, get header */
  if (!fits_read_header (ImageCat, &header)) {
    if (VERBOSE) fprintf (stderr, "New image catalog %s (1)\n", ImageCat);
    *Npimage = 0;
    return ((Image *) NULL);
  }

  /* 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);
  }
  
  /* get ready to read data on images */ 
  f = fopen (ImageCat, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "New image catalog: %s (2)\n", ImageCat);
    free (header.buffer);
    *Npimage = 0;
    return ((Image *) NULL);
  }
  fseek (f, header.size, SEEK_SET); 

  Nimage = 0;
  fits_scan (&header, "NIMAGES", "%d", 1, &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; j++) {
	for (k = 0; k < 4; k++) {
	  addtolist |= edge_check (&Xi[j], &Yi[j], &Xo[k], &Yo[k]);
	}
      }
      
      /*
      fprintf (stderr, "%f %f %f %d %f %f %d %d %d %f %f %s: %d\n", 0.1*timage[i].detection_limit, 0.1*timage[i].saturation_limit, 0.02*timage[i].cerror, 
	       timage[i].tzero, 1e-4*timage[i].trate, 
	       0.001*timage[i].secz, timage[i].nstar, timage[i].NX, timage[i].NY, 0.001*timage[i].Mcal, 0.1*timage[i].Xm, timage[i].name, addtolist);
      
      fprintf (stderr, "%s\n", timage[i].coords.ctype);
      fprintf (stderr, "%f %f\n", timage[i].coords.crval1, timage[i].coords.crval2);
      fprintf (stderr, "%f %f\n", timage[i].coords.crpix1, timage[i].coords.crpix2);
      fprintf (stderr, "%f %f\n", timage[i].coords.pc1_1, timage[i].coords.pc1_2);
      fprintf (stderr, "%f %f\n", timage[i].coords.pc2_1, timage[i].coords.pc2_2);
      fprintf (stderr, "%f %f\n", timage[i].coords.cdelt1, timage[i].coords.cdelt2);      
      */

      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);
      clear_lockfile ();
      exit (0);
    }
  }

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

  fclose (f);
  *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);

}




/* probably not needed now 
char *_parse_nextword();


int fparse (X, NX, line)
float *X;
int NX;
char *line;
{

  int i;
  char *word;
  char *ptr;

  word = line;
  for (i = 0; i < NX - 1; i++)
    word = _parse_nextword (word);

  *X = strtod (word, &ptr);
  if (ptr == word)
    return (FALSE);
  else {
    if (word[0] == '-') 
      return (-1);
    else
      return (1);
  }
}

*/
