# include "addstar.h"
# define NZONE 180
# define NBYTE 207
# define NITEM 100

static SkyTable *tychoTable = NULL;

Stars *gettycho (SkyRegion *catstats, int photcode, unsigned int *nstars) {

  int i, j, nitems;
  char *buffer;
  char filename[128], *line;
  FILE *f;
  int NTYCHO, Ntycho;
  off_t firstRow, firstByte;
  short int TYCHO_B, TYCHO_V;
  e_time TychoEpoch;
  SkyRegion *region;
  SkyList  *skylist;
  Stars *stars;

  /* require photcode */
  NAMED_PHOTCODE (TYCHO_B, "TYCHO_B");
  NAMED_PHOTCODE (TYCHO_V, "TYCHO_V");
  if (photcode == TYCHO_B) goto good_code;
  if (photcode == TYCHO_V) goto good_code;
  Shutdown ("TYCHO photcode not specified");
good_code:

  fprintf (stderr, "loading TYCHO catalog\n");

  /* load tycho index file into sky table */
  if (tychoTable == NULL) {
    sprintf (filename, "%s/index.dat", TYCHO_DIR);
    tychoTable = SkyTableFromTychoIndex (filename, VERBOSE);
    if (tychoTable == NULL) {
      fprintf (stderr, "ERROR: problem loading tycho data\n");
      exit (2);
    }
  }

  /* identify tycho region files overlapping the requested region */
  skylist = SkyListByPatch (tychoTable, 2, catstats);

  /* open the Tycho catalog file */
  sprintf (filename, "%s/tycho.dat", TYCHO_DIR);
  f = fopen (filename, "r");

  Ntycho = 0;
  NTYCHO = 5000;
  ALLOCATE (stars, Stars, NTYCHO);

  ALLOCATE (buffer, char, NITEM*NBYTE);

  TychoEpoch = ohana_date_to_sec ("1991/04/02,07:30:00");

  for (i = 0; i < skylist[0].Nregions; i++) {
    region = skylist[0].regions[i];
    if (VERBOSE) fprintf (stderr, "section %d (%f - %f, %f - %f)...", i, 
			  region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax);

    firstRow  = strtol (region[0].name, NULL, 0);
    firstByte = firstRow * NBYTE;
    fseeko (f, firstByte, SEEK_SET);

    while (1) {
      nitems = fread (buffer, NBYTE, NITEM, f);
      if (nitems == 0) break;

      for (j = 0; j < nitems; j++) {

	line = &buffer[j*NBYTE];

	InitStar (&stars[Ntycho]);

	stars[Ntycho].average.R = atof (&line[15]);
	stars[Ntycho].average.D = atof (&line[28]);

	if (stars[Ntycho].average.R < region[0].Rmin) goto next_section;
	if (stars[Ntycho].average.R > region[0].Rmax) goto next_section;
	if (stars[Ntycho].average.D < region[0].Dmin) goto next_section;
	if (stars[Ntycho].average.D > region[0].Dmax) goto next_section;

	if (stars[Ntycho].average.R < UserPatch.Rmin) continue;
	if (stars[Ntycho].average.R > UserPatch.Rmax) continue;
	if (stars[Ntycho].average.D < UserPatch.Dmin) continue;
	if (stars[Ntycho].average.D > UserPatch.Dmax) continue;

	stars[Ntycho].average.dR  = 1000 * atof (&line[57]);
	stars[Ntycho].average.dD  = 1000 * atof (&line[64]);

	/* XXX : we need to apply uR,uD to R,D to advance to 2000.0 */
	stars[Ntycho].average.uR  = atof (&line[41]);
	stars[Ntycho].average.uD  = atof (&line[49]);

	stars[Ntycho].average.duR = atof (&line[69]);
	stars[Ntycho].average.duD = atof (&line[75]);

	/* Tycho uses J2000 equinox and 1991.25 epoch for coordinates */
	/* the magnitudes have no temporal information */ 
	stars[Ntycho].measure.t   = TychoEpoch;

	/* one pass of addstar does either r or b */
	if (photcode == TYCHO_B) {
	  stars[Ntycho].measure.M     	 = atof (&line[110]);
	  stars[Ntycho].measure.dM    	 = atof (&line[117]);
	  stars[Ntycho].measure.photcode = TYCHO_B;
	} else {
	  stars[Ntycho].measure.M     	 = atof (&line[123]);
	  stars[Ntycho].measure.dM    	 = atof (&line[130]);
	  stars[Ntycho].measure.photcode = TYCHO_V;
	}	
      
	Ntycho ++;
	CHECK_REALLOCATE (stars, Stars, NTYCHO, Ntycho, 5000);
      }

      if (nitems != NITEM) break;
    }
  next_section:
    if (VERBOSE) fprintf (stderr, "%d stars\n", Ntycho);
  }
  fclose (f);

  SkyListFree (skylist);

  free (buffer);
  REALLOCATE (stars, Stars, MAX (1, Ntycho));

  *nstars = Ntycho;
  if (VERBOSE) fprintf (stderr, "%d stars from Tycho\n", Ntycho);
  return (stars);
}
