# include "addstar.h"
# include "supercos.h"

# define NBYTE 228
# define NRECORDS 1000000
# define DEBUG 0

int loadsupercos_rawdata (Image *image, int *imlist, int Nimage, SkyList *skytable, char *filename, AddstarClientOptions options) {
  
  int i, j, verbose;
  int Nstars, NSTARS, Nrecords;
  
  FILE *f;
  char *buffer;

  Detection *sstars;
  Stars **stars, *tstars;
  SkyList *skylist;
  SkyRegion *region;
  Catalog catalog;

  // the data files are written as an integer number of 228 byte records
  ALLOCATE (buffer, char, NRECORDS*NBYTE);

  f = fopen (filename, "r");
  if (f == NULL) Shutdown ("can't read the Supercosmos data file: %s", filename);

  ALLOCATE (tstars, Stars, NRECORDS);
  ALLOCATE (sstars, Detection, NRECORDS); 

  double Rmin = 360.0;
  double Rmax =   0.0;
  double Dmin = +90.0;
  double Dmax = -90.0;

  while ((Nrecords = fread (buffer, NBYTE, NRECORDS, f)) != 0) {
    if (Nrecords == -1) Shutdown ("error reading from raw file %s", filename);
    if (DEBUG) fprintf (stderr, "read %d bytes (%d bytes x %d rows)", NBYTE * Nrecords, NBYTE, Nrecords);

    if (VERBOSE) fprintf (stderr, "read .. ");

    // extract the basic values from the buffer
    for (i = 0; i < Nrecords; i++) {
	sstars[i].objID    = *(int64_t *) &buffer[i*NBYTE + 0];
	sstars[i].surveyID = *(int8_t *)  &buffer[i*NBYTE + 8];
	sstars[i].plateID  = *(int32_t *) &buffer[i*NBYTE + 9];
	sstars[i].ra       = *(double *)  &buffer[i*NBYTE + 33];
	sstars[i].dec      = *(double *)  &buffer[i*NBYTE + 41];
	sstars[i].xCen     = *(double *)  &buffer[i*NBYTE + 129];
	sstars[i].yCen     = *(double *)  &buffer[i*NBYTE + 137];
	sstars[i].aU       = *(float *)   &buffer[i*NBYTE + 145];
	sstars[i].bU       = *(float *)   &buffer[i*NBYTE + 149];
	sstars[i].thetaU   = *(int16_t *) &buffer[i*NBYTE + 154];
	sstars[i].class    = *(int8_t *)  &buffer[i*NBYTE + 165];
	sstars[i].quality  = *(int32_t *) &buffer[i*NBYTE + 204];
	sstars[i].prfStat  = *(float *)   &buffer[i*NBYTE + 208];
	sstars[i].prfMag   = *(float *)   &buffer[i*NBYTE + 212];
	sstars[i].gMag     = *(float *)   &buffer[i*NBYTE + 216];
	sstars[i].sMag     = *(float *)   &buffer[i*NBYTE + 220];
    }

    // convert to the Stars format
    for (i = 0; i < Nrecords; i++) {
	InitStar (&tstars[i]);
	tstars[i].measure.Xccd = sstars[i].xCen;
	tstars[i].measure.Yccd = sstars[i].yCen;

	tstars[i].measure.M = sstars[i].sMag;
	tstars[i].measure.Map = sstars[i].gMag;
	tstars[i].measure.photFlags = sstars[i].class;

	tstars[i].average.R = sstars[i].ra;
	tstars[i].average.D = sstars[i].dec;

	int Ni = imlist[sstars[i].plateID];
	if (Ni == -1) abort();

	// XXX fix these
	tstars[i].measure.Mcal = 0.0;
	tstars[i].measure.dt = image[Ni].exptime;

	double sidtime = image[Ni].sidtime;
	double latitude = image[Ni].latitude;
	
	double alt, az;
	altaz (&alt, &az, 15.0*sidtime - tstars[i].average.R, tstars[i].average.D, latitude);

	tstars[i].measure.airmass = 1.0 / cos(RAD_DEG*alt);
	tstars[i].measure.az = az;
	tstars[i].measure.t = image[Ni].tzero;
	tstars[i].measure.imageID = image[Ni].imageID;
	tstars[i].measure.photcode = image[Ni].photcode;

	tstars[i].measure.FWx = ToShortPixels(sstars[i].aU);
	tstars[i].measure.FWy = ToShortPixels(sstars[i].bU);
	tstars[i].measure.theta = ToShortDegrees(sstars[i].thetaU);

	Rmin = MIN (Rmin, tstars[i].average.R);
	Rmax = MAX (Rmax, tstars[i].average.R);
	Dmin = MIN (Dmin, tstars[i].average.D);
	Dmax = MAX (Dmax, tstars[i].average.D);
    }

    if (VERBOSE) fprintf (stderr, "scan %d stars (%10.6f - %10.6f, %10.6f - %10.6f) .. ", Nrecords, Rmin, Rmax, Dmin, Dmax);

    // sort the tstars by RA
    loadsupercos_sortStars (tstars, Nrecords);

    // scan through the stars, loading the containing catalogs
    // skip through table for unsaved stars
    for (i = 0; i < Nrecords; i++) {
      if (tstars[i].found != -1) continue;

      // scan forward until we reach the UserPatch (stars sorted by RA)
      if (tstars[i].average.R < UserPatch.Rmin) continue;
      if (tstars[i].average.R > UserPatch.Rmax) break;
      if (tstars[i].average.D < UserPatch.Dmin) continue;
      if (tstars[i].average.D > UserPatch.Dmax) continue;

      // identify the relevant catalog
      skylist = SkyRegionByPoint_List (skytable, -1, tstars[i].average.R, tstars[i].average.D);
      if (skylist[0].Nregions == 0) {
	  SkyListFree (skylist);
	  continue;
      }
      region = skylist[0].regions[0];
      if (DEBUG) fprintf (stderr, "writing to %s\n", skylist[0].filename[0]);

      int this_image_id = tstars[i].measure.imageID;

      // loop over stars in this WISE region that are also in this output region
      Nstars = 0;
      NSTARS = 1000;
      ALLOCATE (stars, Stars *, NSTARS);

      for (j = i; j < Nrecords; j++) {
	if (tstars[j].found != -1) continue;

	if (tstars[j].measure.imageID != this_image_id) continue;

	// check if in skyregion
	if (tstars[j].average.R < region[0].Rmin) continue;
	if (tstars[j].average.R > region[0].Rmax) break;
	if (tstars[j].average.D < region[0].Dmin) continue;
	if (tstars[j].average.D > region[0].Dmax) continue;
	  
	// check if in UserPatch
	if (tstars[j].average.R < UserPatch.Rmin) continue;
	if (tstars[j].average.R > UserPatch.Rmax) break;
	if (tstars[j].average.D < UserPatch.Dmin) continue;
	if (tstars[j].average.D > UserPatch.Dmax) continue;

	stars[Nstars] = &tstars[j];
	Nstars ++;

	CHECK_REALLOCATE (stars, Stars *, NSTARS, Nstars, 1000);
      }
	  
      if (DEBUG) fprintf (stderr, "selected %d stars (%10.6f - %10.6f, %10.6f - %10.6f)\n", Nstars, 
			  region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax);

      if (1) {
	verbose = VERBOSE;
	VERBOSE = FALSE;

	// now we have all of the loaded stars in this catalog
	catalog.filename = skylist[0].filename[0];
	catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
	catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
	catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
	catalog.Nsecfilt  = GetPhotcodeNsecfilt ();

	// an error exit status here is a significant error
	if (!dvo_catalog_open (&catalog, skylist[0].regions[0], VERBOSE, "w")) {
	  fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename);
	  exit (2);
	}

	// assume no input star matches an existing star 
	// simply add to the existing table
	// loadWISE_catalog (&catalog, stars, Nstars);
	find_matches_refstars (skylist[0].regions[0], stars, Nstars, &catalog, options);
	// loadWISE_catalog (&catalog, stars, Nstars);

	dvo_catalog_save (&catalog, VERBOSE);
	dvo_catalog_unlock (&catalog);
	dvo_catalog_free (&catalog);
	// free (catalog.filename);
	// XXX don't free this! it points to an element of the skytable
      }

      SkyListFree (skylist);

      free (stars);
      VERBOSE = verbose;
    }
    if (VERBOSE) fprintf (stderr, "done\n");
  }

  if (VERBOSE) fprintf (stderr, "\n");
  
  free (tstars);
  free (sstars);

  fclose (f);
  free (buffer);
  return (TRUE);
}

/*
  for each WISE file:
  for each data block
  generate a table of: R, D, byte, flag
  for each unsaved star
  find containing catalog
  load catalog
  find all contained stars
  add to catalog
  save catalog 
  mark all contained stars
*/

