# include "gastro.h"

# define USNO 1
# define GSC  0

int greference (SStars **cat, int *Ncat, Coords *coords, int NX, int NY) {

  char **regions;
  int  i, Nregions, Nusno;
  CatStats catstats;
  USNOstats usnostats;
  USNOdata *usno;

  if (VERBOSE) fprintf (stderr, "\nloading astrometric reference data from %s\n", REFCAT); 

  *Ncat = 0;

  define_region (&catstats, coords, NX, NY);

  /* get stars from the USNO catalog for the given region */
  if (!strcmp (REFCAT, "USNO") || !strcmp (REFCAT, "BOTH")) {
    while (catstats.RA[0] > 360.0) catstats.RA[0] -= 360.0;
    while (catstats.RA[0] <   0.0) catstats.RA[0] += 360.0;
    while (catstats.RA[1] > 360.0) catstats.RA[1] -= 360.0;
    while (catstats.RA[1] <   0.0) catstats.RA[1] += 360.0;

    /* if RA crosses 0,360 boundary, do 2 passes */
    if (catstats.RA[0] > catstats.RA[1]) {
        CatStats substats;
	USNOdata *usno1, *usno2;
	int Nusno1, Nusno2;

	substats = catstats;
	substats.RA[0] = 0.0;
	usno1 = getusno (&usnostats, &substats, &Nusno1);

	substats = catstats;
	substats.RA[1] = 360.0;
	usno2 = getusno (&usnostats, &substats, &Nusno2);

	ALLOCATE (*cat, SStars, MAX (Nusno1 + Nusno2, 1));
	for (i = 0; i < Nusno1; i++) {
	  cat[0][i].X = usno1[i].R;
	  cat[0][i].Y = usno1[i].D;
	  cat[0][i].mag = fabs(usno1[i].r);
	}      
	for (i = Nusno1; i < Nusno2; i++) {
	  cat[0][i].X = usno2[i].R;
	  cat[0][i].Y = usno2[i].D;
	  cat[0][i].mag = fabs(usno2[i].r);
	}      
	*Ncat = Nusno1 + Nusno2;
	free (usno1);
	free (usno2);
    } else {
	usno = getusno (&usnostats, &catstats, &Nusno);

	ALLOCATE (*cat, SStars, MAX (Nusno, 1));
	for (i = 0; i < Nusno; i++) {
	  cat[0][i].X = usno[i].R;
	  cat[0][i].Y = usno[i].D;
	  cat[0][i].mag = fabs(usno[i].r);
	}      
	*Ncat = Nusno;
	free (usno);
    }
    if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", *Ncat);
  }

  if (!strcmp (REFCAT, "GSC") || !strcmp (REFCAT, "BOTH")) {
    regions = gregions (&catstats, &Nregions);
    for (i = 0; i < Nregions; i++) {
      if (!gcatalog (regions[i], cat, Ncat))       
	return (FALSE);
    }
    if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", *Ncat);
  }

  if (!strcmp (REFCAT, "PTOLEMY")) {
    regions = gregions (&catstats, &Nregions);
    for (i = 0; i < Nregions; i++) {
      if (!gptolemy (regions[i], cat, Ncat)) continue;
    }
    if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", *Ncat);
  }

  if (!strcmp (REFCAT, "2MASS")) {
    /* this needs to be written to handle the complete 2mass catalog.
       at the moment, I only have a small subset */
    if (!g2mass ("/data/milo/eugene/elixir/refs/2mass/2mass.dat", cat, Ncat)) exit (1);
    if (VERBOSE) fprintf (stderr, "%d stars from 2MASS\n", *Ncat);
  }

  return (TRUE);

}

void define_region (CatStats *catstats, Coords *coords, int NX, int NY) {
   
  int i, j;
  double X, Y, R, D;

  catstats[0].RA[0] = catstats[0].DEC[0] =  360.0;
  catstats[0].RA[1] = catstats[0].DEC[1] = -360.0;

  for (i = -1; i < 2; i++) {
    for (j = -1; j < 2; j++) {
      X = 0.5*(1.0 + i*NFIELD)*NX;
      Y = 0.5*(1.0 + j*NFIELD)*NY;
      XY_to_RD (&R, &D, X, Y, coords);
      /* coords returns a region all in same phase 
      while (R < 0.0)    R += 360.0;
      while (R >= 360.0) R -= 360.0;
      */
      catstats[0].RA[0]  = MIN (catstats[0].RA[0], R);
      catstats[0].RA[1]  = MAX (catstats[0].RA[1], R);
      catstats[0].DEC[0] = MIN (catstats[0].DEC[0], D);
      catstats[0].DEC[1] = MAX (catstats[0].DEC[1], D);
    }
  }
  fprintf (stderr, "full region: %f - %f, %f - %f\n", 
	   catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]);
}
