# include "gastro2.h"

void greference (CmpCatalog *Target, RefCatalog *Ref) {

  CatStats catstats;

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

  define_region (&catstats, Target);

  Ref[0].N = 0;
  /* get stars from the USNO catalog for the given region */
  if (!strcmp (REFCAT, "USNO")) {
    getusno (&catstats, Ref);
    /* calculate Ref[0].Moff from Target & Ref dMdN, Mo */
  }

  /* get stars from the HST GSC catalog for the given region */
  if (!strcmp (REFCAT, "GSC")) {
    getgsc (&catstats, Ref);
  }
  
  /* get stars from the HST GSC catalog for the given region */
  if (!strcmp (REFCAT, "2MASS")) {
    get2mass (&catstats, Ref);
  }
  
  /* get stars from the HST GSC catalog for the given region */
  if (!strcmp (REFCAT, "PTOLEMY")) {
    getptolemy (&catstats, Ref);
  }
  
  if (Ref[0].N == 0) {
    fprintf (stderr, "no ref objs: %s\n", REFCAT);
    exit (1);
  }

  Ref[0].lum.Mz = (2 + log (Ref[0].Area) / log (10.) - Ref[0].lum.Mo) / Ref[0].lum.dNdM;
  Target[0].lum.Mz = (2 + log (Target[0].Area) / log (10.) - Target[0].lum.Mo) / Target[0].lum.dNdM;
  Ref[0].Moff = Target[0].lum.Mz - Ref[0].lum.Mz;
  fprintf (stderr, "mag offset: %f (Areas: %f vs %f)\n", Ref[0].Moff, Target[0].Area, Ref[0].Area);

  if (PLOTSTUFF) plot_lumfunc (Target, Ref);

}

/* return RA, DEC bounds of the reigon of interest */  
void define_region (CatStats *catstats, CmpCatalog *Target) {
   
  int NX, NY;
  double x, y, X, Y, R, D, dX, dY, Xo, Yo;

  NX = Target[0].header.Naxis[0];
  NY = Target[0].header.Naxis[1];

  dX = NX + NFIELD*NX;
  dY = NY + NFIELD*NY;

  Xo = -0.5*NFIELD*NX;
  Yo = -0.5*NFIELD*NY;

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

  for (x = 0; x <= 1.0; x += 0.5) {
    for (y = 0; y <= 1.0; y += 0.5) {

      X = x*dX + Xo;
      Y = y*dY + Yo;
      XY_to_RD (&R, &D, X, Y, &Target[0].coords);

      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);
    }
  }
  if (VERBOSE) 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]);
}
