# include "uniphot.h"

Group *find_image_sgroups (FITS_DB *db, ImageLink **Imlink, int *Nsgroup) {

  int i, j, Nimage, Ngroup, Nentry, NENTRY;
  double r, d, x, y, radius;
  Group *group;
  Coords coords;
  Image *image;
  ImageLink *imlink;

  imlink = *Imlink;
  image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);

  coords.crpix1 = coords.crpix2 = 0.0;
  coords.cdelt1 = coords.cdelt2 = 1.0; /* scale is degrees, radius in degrees */
  coords.pc1_1 = coords.pc2_2 = 1.0;
  coords.pc1_2 = coords.pc2_1 = 0.0;
  strcpy (coords.ctype, "RA---TAN");
  
  Ngroup = 0;
  ALLOCATE (group, Group, Nimage);

  /* set imlink.sgroups = NULL as a marker */
  for (i = 0; i < Nimage; i++) imlink[i].sgroup = NULL;

  for (i = 0; i < Nimage; i++) {
    if (imlink[i].sgroup != NULL) continue;
    if (image[i].flags & ID_IMAGE_NOCAL) continue;

    /* new sgroup, set ref coords */
    XY_to_RD (&r, &d, 0.5*image[i].NX, 0.5*image[i].NX, &image[i].coords);
    coords.crval1 = r;
    coords.crval2 = d;

    /* init sgroup structure */
    Nentry = 0;
    NENTRY = 100;
    ALLOCATE (group[Ngroup].image, Image *, NENTRY);
    ALLOCATE (group[Ngroup].imlink, ImageLink *, NENTRY);
    group[Ngroup].M = 0;
    group[Ngroup].dM = 0;
    snprintf (group[Ngroup].label, 64, "%10.6f - %10.6f", r, d);
    group[Ngroup].v1 = r;
    group[Ngroup].v2 = d;

    /* link this image to sgroup */
    group[Ngroup].image[Nentry] = &image[i];
    group[Ngroup].imlink[Nentry] = &imlink[i];
    imlink[i].sgroup = &group[Ngroup];
    Nentry ++;

    for (j = 0; j < Nimage; j++) {
      if (image[j].flags & ID_IMAGE_NOCAL) continue;
      if (imlink[j].sgroup != NULL) continue;

      /* project image center to local coords, check radius */
      XY_to_RD (&r, &d, 0.5*image[j].NX, 0.5*image[j].NX, &image[j].coords);
      if (!RD_to_XY (&x, &y, r, d, &coords)) continue; 
      /* RD_to_XY returns FALSE if opposite hemispheres */
      radius = hypot (x, y);

      if (radius > RADIUS) continue;

      /* image in sgroup, add to entry */
      group[Ngroup].image[Nentry] = &image[j];
      group[Ngroup].imlink[Nentry] = &imlink[j];
      imlink[j].sgroup = &group[Ngroup];
      Nentry ++;
      if (Nentry == NENTRY) {
	NENTRY += 100;
	REALLOCATE (group[Ngroup].image, Image *, NENTRY);
	REALLOCATE (group[Ngroup].imlink, ImageLink *, NENTRY);
      }
    }
    group[Ngroup].Nimage = Nentry;
    Ngroup ++;
  }
  *Nsgroup = Ngroup;
  return (group);
}

  /* this is a bit weak: since we use pointers, we can't
     reallocate group after the pointers are assigned.
     therefore, we allocate the max possible groups */
