# include "relphot.h"

// relphot_client reads from the local catalogs and generates a set of bright, subset catalogs
// we are going to save these as a single big FITS file, with only 3 extensions:

// measure, average, secfilt.  these only need to contain the fields relevant to the
// MeasureTiny, AverageTiny, and Secfilt tables, but with cat_id appended to each row so we can
// reassign correctly on read

# define GET_COLUMN(OUT,NAME,TYPE)					\
  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
  myAssert (!strcmp(type, #TYPE), "wrong column type");

// XXX double check the header extname values for each table
// XXX make sure we free things as we can
// XXX make sure we close files as we can
// XXX handle free and close on error return as well
BrightCatalog *BrightCatalogLoad(char *filename) {

  int i, Ncol;
  off_t Nrow;
  char type[16];
  Header header;
  Header theader;
  Matrix matrix;
  FTable ftable;

  header.buffer = NULL;
  matrix.buffer = NULL;
  ftable.buffer = NULL;
  theader.buffer = NULL;
  BrightCatalog *catalog = NULL;

  FILE *f = fopen (filename, "r");
  if (!f) {
    fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename);
    return NULL;
  }

  /* load in PHU segment (ignore) */
  if (!gfits_fread_header (f, &header)) {
    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
    goto escape;
  }
  if (!gfits_fread_matrix (f, &matrix, &header)) {
    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
    goto escape;
  }

  ALLOCATE (catalog, BrightCatalog, 1);
  catalog->Naverage = 0;
  catalog->Nmeasure = 0;

  ftable.header = &theader;

  // *** Measure Data *** 
  { 
    // load data for this header 
    if (!gfits_load_header (f, &theader)) goto escape;

    // read the fits table bytes
    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
 
    // need to create and assign to flat-field correction
    GET_COLUMN(dR,        "RA_OFF",   	float);
    GET_COLUMN(dD,        "DEC_OFF",  	float);
    GET_COLUMN(M,         "MAG_SYS",  	float);
    GET_COLUMN(Mcal,      "MAG_CAL",  	float);
    GET_COLUMN(dM,        "MAG_ERR",  	float);
    GET_COLUMN(airmass,   "AIRMASS",  	float);
    GET_COLUMN(Xccd,      "X_CCD",    	float);
    GET_COLUMN(Yccd,      "Y_CCD",    	float);
    GET_COLUMN(dt,        "EXPTIME",  	float);
    GET_COLUMN(t,         "TIME",     	int);
    GET_COLUMN(averef,    "AVE_REF",  	int); // XXX signed vs unsigned?
    GET_COLUMN(imageID,   "IMAGE_ID", 	int);
    GET_COLUMN(dbFlags,   "DB_FLAGS", 	int);
    GET_COLUMN(photFlags, "PHOT_FLAGS", int);
    GET_COLUMN(catID,     "CAT_ID",     int);
    GET_COLUMN(photcode,  "PHOTCODE",   short);
    gfits_free_header (&theader);
    gfits_free_table  (&ftable);

    MeasureTiny *measure = NULL;
    ALLOCATE (measure, MeasureTiny, Nrow);
    for (i = 0; i < Nrow; i++) {
      measure[i].dR        = dR[i];
      measure[i].dD        = dD[i];
      measure[i].M         = M[i];
      measure[i].Mcal      = Mcal[i];
      measure[i].dM        = dM[i];
      measure[i].airmass   = airmass[i];
      measure[i].Xccd      = Xccd[i];
      measure[i].Yccd      = Yccd[i];
      measure[i].dt        = dt[i];
      measure[i].t         = t[i];
      measure[i].averef    = averef[i];
      measure[i].imageID   = imageID[i];
      measure[i].dbFlags   = dbFlags[i];
      measure[i].photFlags = photFlags[i];
      measure[i].catID     = catID[i];
      measure[i].photcode  = photcode[i];
      measure[i].myDet     = FALSE;
    }
    fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow);

    free (dR      );
    free (dD      );
    free (M       );
    free (Mcal    );
    free (dM      );
    free (airmass );
    free (Xccd    );
    free (Yccd    );
    free (dt      );
    free (t       );
    free (averef  );
    free (imageID );
    free (dbFlags );
    free (photFlags );
    free (catID   );
    free (photcode);

    catalog->measure = measure;
    catalog->Nmeasure = Nrow;
  }

  /*** load Average values ***/
  { 

    // load data for this header 
    if (!gfits_load_header (f, &theader)) goto escape;

    // read the fits table bytes
    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
 
    // need to create and assign to flat-field correction
    GET_COLUMN(R,             "RA",   	   double);
    GET_COLUMN(D,             "DEC",  	   double);
    GET_COLUMN(Nmeasure,      "NMEAS",       int);
    GET_COLUMN(measureOffset, "MEASURE_OFF", int);
    GET_COLUMN(flags,         "FLAGS",       int);
    GET_COLUMN(catID,         "CAT_ID",      int);
    GET_COLUMN(objID,         "OBJ_ID",      int);
    gfits_free_header (&theader);
    gfits_free_table  (&ftable);

    AverageTiny *average = NULL;
    ALLOCATE (average, AverageTiny, Nrow);
    for (i = 0; i < Nrow; i++) {
      average[i].R              = R[i];
      average[i].D              = D[i];
      average[i].Nmeasure       = Nmeasure[i];
      average[i].measureOffset  = measureOffset[i];
      average[i].flags          = flags[i];
      average[i].catID          = catID[i];
      average[i].objID          = objID[i];
      average[i].nOwn           = 0;
    }
    fprintf (stderr, "loaded data for %lld average\n", (long long) Nrow);

    free (R             );
    free (D             );
    free (Nmeasure      );
    free (measureOffset );
    free (flags         );
    free (catID         );
    free (objID         );

    catalog->average = average;
    catalog->Naverage = Nrow;
  }

  /*** load Secfilt values ***/
  {

    // load data for this header 
    if (!gfits_load_header (f, &theader)) goto escape;

    // read the fits table bytes
    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
 
    // need to create and assign to flat-field correction
    GET_COLUMN(M,     "MAG",      float);
    GET_COLUMN(dM,    "MAG_ERR",  float);
    GET_COLUMN(Xm,    "MAG_CHI",  float);
    GET_COLUMN(flags, "FLAGS",    int);
    GET_COLUMN(Ncode, "NCODE",    short);
    GET_COLUMN(Nused, "NUSED",    short);
    GET_COLUMN(M_20,  "MAG_20",   short);
    GET_COLUMN(M_80,  "MAG_80",   short);
    gfits_free_header (&theader);
    gfits_free_table  (&ftable);

    SecFilt *secfilt = NULL;
    ALLOCATE (secfilt, SecFilt, Nrow);
    for (i = 0; i < Nrow; i++) {
      secfilt[i].M     = M[i];         
      secfilt[i].dM    = dM[i];
      secfilt[i].Xm    = Xm[i];
      secfilt[i].flags = flags[i];
      secfilt[i].Ncode = Ncode[i];
      secfilt[i].Nused = Nused[i];
      secfilt[i].M_20  = M_20[i];
      secfilt[i].M_80  = M_80[i];
    }
    fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow);

    free (M    );
    free (dM   );
    free (Xm   );
    free (flags);
    free (Ncode);
    free (Nused);
    free (M_20 );
    free (M_80 );
    catalog->secfilt = secfilt;
    // assert Nsecfilt * Naverage = Nrow?
  }

  gfits_free_header (&header);
  gfits_free_matrix (&matrix);
  fclose (f);

  return catalog;

escape:
  gfits_free_header (&header);
  gfits_free_matrix (&matrix);
  gfits_free_header (&theader);
  gfits_free_table  (&ftable);
  if (catalog) free (catalog);

  fclose (f);
  return NULL;
}

// we are passed a BrightCatalog structure, write it to a FITS table (3 ext)
int BrightCatalogSave(char *filename, BrightCatalog *catalog) {

  int i;
  Header header;
  Header theader;
  Matrix matrix;
  FTable ftable;

  gfits_init_header (&header);
  header.extend = TRUE;
  gfits_create_header (&header);
  gfits_create_matrix (&header, &matrix);

  FILE *f = fopen (filename, "w");
  if (!f) {
    fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename);
    return FALSE;
  }

  gfits_fwrite_header  (f, &header);
  gfits_fwrite_matrix  (f, &matrix);
  gfits_free_header (&header);
  gfits_free_matrix (&matrix);

  /*** MeasureTiny ***/
  {
    gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY");

    gfits_define_bintable_column (&theader, "E", "RA_OFF",   "ra offset",                  "arcsec", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "DEC_OFF",  "dec offset",                 "arcsec", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "magnitude (err)",             NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "AIRMASS",  "airmass",                	    NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "X_CCD",    "ccd x coord",            	   "pix",    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "Y_CCD",    "ccd y coord",            	   "pix",    1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "EXPTIME",  "-2.5 * log (exposure time)", "sec",    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "TIME",     "time of exp",                "sec",    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "J", "AVE_REF",  "pointer to average table",    NULL,    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image",                       NULL,    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags",                       NULL,    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags",                 NULL,    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "J", "CAT_ID",   "catalog",                     NULL,    1.0, 1.0*0x8000);
    gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode",                    NULL,    1.0, 1.0*0x80);

    // generate the output array that carries the data
    gfits_create_table (&theader, &ftable);

    // create intermediate storage arrays
    float *dR        ; ALLOCATE (dR       ,  float, catalog->Nmeasure);
    float *dD        ; ALLOCATE (dD       ,  float, catalog->Nmeasure);
    float *M         ; ALLOCATE (M        ,  float, catalog->Nmeasure);
    float *Mcal      ; ALLOCATE (Mcal     ,  float, catalog->Nmeasure);
    float *dM        ; ALLOCATE (dM       ,  float, catalog->Nmeasure);
    float *airmass   ; ALLOCATE (airmass  ,  float, catalog->Nmeasure);
    float *Xccd      ; ALLOCATE (Xccd     ,  float, catalog->Nmeasure);
    float *Yccd      ; ALLOCATE (Yccd     ,  float, catalog->Nmeasure);
    float *dt        ; ALLOCATE (dt       ,  float, catalog->Nmeasure);
    int   *t         ; ALLOCATE (t        ,  int  , catalog->Nmeasure);
    int   *averef    ; ALLOCATE (averef   ,  int  , catalog->Nmeasure);
    int   *imageID   ; ALLOCATE (imageID  ,  int  , catalog->Nmeasure);
    int   *dbFlags   ; ALLOCATE (dbFlags  ,  int  , catalog->Nmeasure);
    int   *photFlags ; ALLOCATE (photFlags,  int  , catalog->Nmeasure);
    int   *catID     ; ALLOCATE (catID    ,  int  , catalog->Nmeasure);
    short *photcode  ; ALLOCATE (photcode ,  short, catalog->Nmeasure);

    // assign the storage arrays
    MeasureTiny *measure = catalog->measure;
    for (i = 0; i < catalog->Nmeasure; i++) {
      dR[i]       = measure[i].dR       ;
      dD[i]       = measure[i].dD       ;
      M[i]  	  = measure[i].M        ;
      Mcal[i]     = measure[i].Mcal     ;
      dM[i]       = measure[i].dM       ;
      airmass[i]  = measure[i].airmass  ;
      Xccd[i]     = measure[i].Xccd     ;
      Yccd[i]     = measure[i].Yccd     ;
      dt[i]       = measure[i].dt       ;
      t[i] 	  = measure[i].t        ;
      averef[i]   = measure[i].averef   ;
      catID[i]    = measure[i].catID    ;
      imageID[i]  = measure[i].imageID  ;
      dbFlags[i]  = measure[i].dbFlags  ;
      photFlags[i]= measure[i].photFlags;
      photcode[i] = measure[i].photcode ;
    }

    // add the columns to the output array
    gfits_set_bintable_column (&theader, &ftable, "RA_OFF",   	dR,        catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "DEC_OFF",  	dD,        catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",  	M,         catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",  	Mcal,      catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  	dM,        catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "AIRMASS",  	airmass,   catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "X_CCD",    	Xccd,      catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "Y_CCD",    	Yccd,      catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "EXPTIME",  	dt,        catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "TIME",     	t,         catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "AVE_REF",  	averef,    catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", 	imageID,   catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", 	dbFlags,   catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "CAT_ID",     catID,     catalog->Nmeasure);
    gfits_set_bintable_column (&theader, &ftable, "PHOTCODE",   photcode,  catalog->Nmeasure);

    free (dR      );
    free (dD      );
    free (M       );
    free (Mcal    );
    free (dM      );
    free (airmass );
    free (Xccd    );
    free (Yccd    );
    free (dt      );
    free (t       );
    free (averef  );
    free (imageID );
    free (dbFlags );
    free (photFlags );
    free (catID   );
    free (photcode);

    gfits_fwrite_Theader (f, &theader);
    gfits_fwrite_table  (f, &ftable);
    gfits_free_header (&theader);
    gfits_free_table (&ftable);
  }

  /*** AverageTiny ***/
  {
    gfits_create_table_header (&theader, "BINTABLE", "AVERAGE_TINY");

    gfits_define_bintable_column (&theader, "D", "RA",   	"ra (J2000)", 	         "degree", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "D", "DEC",  	"dec (J2000)",	         "degree", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "NMEAS",       "number of measures",     NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "MEASURE_OFF", "index to measurements",  NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "FLAGS",       "flags",                  NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "CAT_ID",      "catalog ref",            NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "OBJ_ID",      "object ref",             NULL,    1.0, 0.0);

    // generate the output array that carries the data
    gfits_create_table (&theader, &ftable);

    // create intermediate storage arrays
    double *R             ; ALLOCATE (R,             double, catalog->Naverage);
    double *D             ; ALLOCATE (D,             double, catalog->Naverage);
    int   *Nmeasure       ; ALLOCATE (Nmeasure,      int,    catalog->Naverage);
    int   *measureOffset  ; ALLOCATE (measureOffset, int,    catalog->Naverage);
    int   *flags          ; ALLOCATE (flags,         int,    catalog->Naverage);
    int   *catID          ; ALLOCATE (catID,         int,    catalog->Naverage);
    int   *objID          ; ALLOCATE (objID,         int,    catalog->Naverage);

    // assign the storage arrays
    AverageTiny *average = catalog->average;
    for (i = 0; i < catalog->Naverage; i++) {
      R[i]           	= average[i].R       ;
      D[i]           	= average[i].D       ;
      Nmeasure[i]    	= average[i].Nmeasure;
      measureOffset[i]  = average[i].measureOffset;
      flags[i]          = average[i].flags;
      catID[i]          = average[i].catID;
      objID[i]          = average[i].objID;
    }

    // add the columns to the output array
    gfits_set_bintable_column (&theader, &ftable, "RA",          R,             catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "DEC",         D,             catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "NMEAS",       Nmeasure,      catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "MEASURE_OFF", measureOffset, catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "FLAGS",       flags,         catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "CAT_ID",      catID,         catalog->Naverage);
    gfits_set_bintable_column (&theader, &ftable, "OBJ_ID",      objID,         catalog->Naverage);

    free (R             );
    free (D             );
    free (Nmeasure      );
    free (measureOffset );
    free (flags         );
    free (catID         );
    free (objID         );

    gfits_fwrite_Theader (f, &theader);
    gfits_fwrite_table  (f, &ftable);
    gfits_free_header (&theader);
    gfits_free_table (&ftable);
  }

  /*** SecFilt ***/
  {
    gfits_create_table_header (&theader, "BINTABLE", "SECFILT");

    gfits_define_bintable_column (&theader, "E", "MAG",      "ra offset",                "arcsec", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "dec offset",               "arcsec", 1.0, 0.0);
    gfits_define_bintable_column (&theader, "E", "MAG_CHI",  "magnitude (sys)",           NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "J", "FLAGS",    "magnitude (cal)",           NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "I", "NCODE",    "magnitude (err)",           NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "I", "NUSED",    "airmass",                   NULL,    1.0, 0.0);
    gfits_define_bintable_column (&theader, "I", "MAG_20",   "ccd x coord",              "pix",    1.0, 0.0);
    gfits_define_bintable_column (&theader, "I", "MAG_80",   "ccd y coord",              "pix",    1.0, 0.0);

    // generate the output array that carries the data
    gfits_create_table (&theader, &ftable);

    // Nsecfilt is number of average filters, Nsec is total number of values
    int Nsecfilt = GetPhotcodeNsecfilt();
    int Nsec = Nsecfilt * catalog->Naverage;

    // create intermediate storage arrays
    float *M        ; ALLOCATE (M      ,  float, Nsec);
    float *dM       ; ALLOCATE (dM     ,  float, Nsec);
    float *Xm       ; ALLOCATE (Xm     ,  float, Nsec);
    int   *flags    ; ALLOCATE (flags  ,  int,   Nsec);
    short *Ncode    ; ALLOCATE (Ncode  ,  short, Nsec);
    short *Nused    ; ALLOCATE (Nused  ,  short, Nsec);
    short *M_20     ; ALLOCATE (M_20   ,  short, Nsec);
    short *M_80     ; ALLOCATE (M_80   ,  short, Nsec);

    // assign the storage arrays
    SecFilt *secfilt = catalog->secfilt;
    for (i = 0; i < Nsec; i++) {
      M     [i]       = secfilt[i]. M      ;
      dM    [i]       = secfilt[i]. dM     ;
      Xm    [i]       = secfilt[i]. Xm     ;
      flags [i]       = secfilt[i]. flags  ;
      Ncode [i]       = secfilt[i]. Ncode  ;
      Nused [i]       = secfilt[i]. Nused  ;
      M_20  [i]       = secfilt[i]. M_20   ;
      M_80  [i]       = secfilt[i]. M_80   ;
    }

    // add the columns to the output array
    gfits_set_bintable_column (&theader, &ftable, "MAG",      M    , Nsec);
    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM   , Nsec);
    gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Xm   , Nsec);
    gfits_set_bintable_column (&theader, &ftable, "FLAGS",    flags, Nsec);
    gfits_set_bintable_column (&theader, &ftable, "NCODE",    Ncode, Nsec);
    gfits_set_bintable_column (&theader, &ftable, "NUSED",    Nused, Nsec);
    gfits_set_bintable_column (&theader, &ftable, "MAG_20",   M_20 , Nsec);
    gfits_set_bintable_column (&theader, &ftable, "MAG_80",   M_80 , Nsec);

    free (M      );
    free (dM     );
    free (Xm     );
    free (flags  );
    free (Ncode  );
    free (Nused  );
    free (M_20   );
    free (M_80   );

    gfits_fwrite_Theader (f, &theader);
    gfits_fwrite_table  (f, &ftable);
    gfits_free_header (&theader);
    gfits_free_table (&ftable);
  }
  return TRUE;
}

// merge a list of catalogs into a single BrightCatalog array set
BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog) {

  off_t i, j, k;

  BrightCatalog *bcatalog = NULL;
  ALLOCATE (bcatalog, BrightCatalog, 1);
  
  int Nmeas = 0;
  int Naves = 0;
  for (i = 0; i < Ncatalog; i++) {
    Nmeas += catalog[i].Nmeasure;
    Naves += catalog[i].Naverage;
  }
    
  // XXX this prevents different catalogs from having different Nsecfilt values
  int Nsecfilt = GetPhotcodeNsecfilt();

  ALLOCATE (bcatalog[0].measure, MeasureTiny, Nmeas);
  ALLOCATE (bcatalog[0].average, AverageTiny, Naves);
  ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt);

  int Nm = 0;
  int Na = 0;
  for (i = 0; i < Ncatalog; i++) {
    if (!catalog[i].Naverage) continue;
    for (j = 0; j < catalog[i].Naverage; j++) {
      bcatalog[0].average[Na] = catalog[i].averageT[j];
      for (k = 0; k < Nsecfilt; k++) {
	bcatalog[0].secfilt[Nsecfilt*Na + k] = catalog[i].secfilt[Nsecfilt*j + k];
      }
      Na++;
      assert (Na <= Naves);
    }
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      bcatalog[0].measure[Nm] = catalog[i].measureT[j];
      Nm++;
      assert (Nm <= Nmeas);
    }
  }
  bcatalog->Naverage = Na;
  bcatalog->Nmeasure = Nm;
  return bcatalog;
}

// distribute a bright catalog across separate catalogs
CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt) {

  int i;

  CatalogSplitter *catalogs = NULL;

  ALLOCATE (catalogs, CatalogSplitter, 1);

  // as we see new BrightCatalogs, we will update maxID and extend this array
  catalogs->maxID = 0;
  ALLOCATE (catalogs->index, int, catalogs->maxID + 1);
  for (i = 0; i <= catalogs->maxID; i++) catalogs->index[i] = -1;

  catalogs->Nsecfilt = Nsecfilt;

  catalogs->Ncatalog =  0;
  catalogs->NCATALOG = 16;
  catalogs->catalog = NULL;
  ALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG);

  ALLOCATE (catalogs->catIDs,   int,   catalogs->NCATALOG);
  ALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG);
  ALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG);

  for (i = 0; i < catalogs->NCATALOG; i++) {
    dvo_catalog_init (&catalogs->catalog[i], TRUE);
    catalogs->catIDs[i] = 0;
    catalogs->NAVERAGE[i] = 256;
    catalogs->NMEASURE[i] = 256;
    catalogs->catalog[i].Naverage = 0;
    catalogs->catalog[i].Nmeasure = 0;
    ALLOCATE (catalogs->catalog[i].averageT, AverageTiny, catalogs->NAVERAGE[i]);
    ALLOCATE (catalogs->catalog[i].measureT, MeasureTiny, catalogs->NMEASURE[i]);
    ALLOCATE (catalogs->catalog[i].secfilt,  SecFilt,     catalogs->NAVERAGE[i]*Nsecfilt);
  }
  return catalogs;
}

// distribute a bright catalog across separate catalogs
int BrightCatalogSplitFree (CatalogSplitter *catalogs) {

  // XXX don't free the catalogs : 
  // for (i = 0; i < NCATALOG; i++) {
  //   free (catalogs->catalog[i].averageT);
  //   free (catalogs->catalog[i].measureT);
  //   free (catalogs->catalog[i].secfilt);
  // }
  // free (catalogs->catalog);

  free (catalogs->catIDs);
  free (catalogs->NAVERAGE);
  free (catalogs->NMEASURE);
  free (catalogs->index);
  free (catalogs);
  return TRUE;
}

// distribute a bright catalog across separate catalogs
int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog) {

  int i;

  int D_NCATALOG = 16;

  int Nsecfilt = catalogs->Nsecfilt;

  // find the max value of catID in this BrightCatalog
  int catIDmax = 0;
  for (i = 0; i < bcatalog->Naverage; i++) {
    catIDmax = MAX(catIDmax, bcatalog->average[i].catID);
  }
  // XXX validate the measure ID range here
    
  int maxIDold = catalogs->maxID;
  catalogs->maxID = MAX (maxIDold, catIDmax);
    
  // extend the index array and init
  REALLOCATE (catalogs->index, int, catalogs->maxID + 1);
  for (i = maxIDold + 1; i <= catalogs->maxID; i++) catalogs->index[i] = -1;

  // identify the new catID values
  for (i = 0; i < bcatalog->Naverage; i++) {
    int catID = bcatalog->average[i].catID;
    assert (catID > 0);
    assert (catID < catalogs->maxID + 1);
    int idx = catalogs->index[catID];
    if (idx != -1) continue; // already have seen this one

    // a new catID value
    int Ncat = catalogs->Ncatalog; // the next available slot
    catalogs->catIDs[Ncat] = catID;
    assert (Ncat >= 0);
    assert (Ncat < catalogs->NCATALOG);

    catalogs->catalog[Ncat].Nsecfilt = Nsecfilt;
    catalogs->catalog[Ncat].catID = catID;

    catalogs->index[catID] = Ncat;
    assert (catID > 0);
    assert (catID < catalogs->maxID + 1);

    catalogs->Ncatalog ++;

    if (catalogs->Ncatalog >= catalogs->NCATALOG) {
      catalogs->NCATALOG += D_NCATALOG;

      // fprintf (stderr, "realloc catalogs->catalog: old: %llx  ", (long long) catalogs->catalog);
      REALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG);
      // fprintf (stderr, "new: %llx  -  %llx\n", (long long) catalogs->catalog, (long long) (catalogs->catalog + sizeof(Catalog)*catalogs->NCATALOG));

      // fprintf (stderr, "realloc catalogs->NAVERAGE: old: %llx  ", (long long) catalogs->NAVERAGE);
      REALLOCATE (catalogs->NAVERAGE, off_t,  catalogs->NCATALOG);
      // fprintf (stderr, "new: %llx  -  %llx\n", (long long) catalogs->NAVERAGE, (long long) (catalogs->NAVERAGE + sizeof(off_t)*catalogs->NCATALOG));

      // fprintf (stderr, "realloc catalogs->NMEASURE: old: %llx  ", (long long) catalogs->NMEASURE);
      REALLOCATE (catalogs->NMEASURE, off_t,  catalogs->NCATALOG);
      // fprintf (stderr, "new: %llx  -  %llx\n", (long long) catalogs->NMEASURE, (long long) (catalogs->NMEASURE + sizeof(off_t)*catalogs->NCATALOG));

      REALLOCATE (catalogs->catIDs,   int,   catalogs->NCATALOG);

      int j;
      for (j = catalogs->NCATALOG - D_NCATALOG; j < catalogs->NCATALOG; j++) {
	dvo_catalog_init (&catalogs->catalog[j], TRUE);
	catalogs->catIDs[j] = 0;
	catalogs->NAVERAGE[j] = 256;
	catalogs->NMEASURE[j] = 256;
	catalogs->catalog[j].Naverage = 0;
	catalogs->catalog[j].Nmeasure = 0;
	ALLOCATE (catalogs->catalog[j].averageT, AverageTiny, catalogs->NAVERAGE[j]);
	ALLOCATE (catalogs->catalog[j].measureT, MeasureTiny, catalogs->NMEASURE[j]);
	ALLOCATE (catalogs->catalog[j].secfilt,  SecFilt,     catalogs->NAVERAGE[j]*Nsecfilt);
      }
      D_NCATALOG = MAX(2000, 2*D_NCATALOG);
    }
  }

  // each bcatalog (from each host) has a different set of catID ranges
  // they also (probably) have contiguous ranges.  But, it is a bit tricky to be sure I
  // can use that info, so perhaps ignore it

  // assign the averages to the corresponding catalog
  for (i = 0; i < bcatalog->Naverage; i++) {
    int ID = bcatalog->average[i].catID;
    int Nc = catalogs->index[ID];
    assert (Nc > -1);
    assert (Nc < catalogs->NCATALOG);

    int Na = catalogs->catalog[Nc].Naverage;
    catalogs->catalog[Nc].averageT[Na] = bcatalog->average[i];

    // secfilt entries are grouped in blocks for each average
    int k;
    for (k = 0; k < Nsecfilt; k++) {
      catalogs->catalog[Nc].secfilt[Na*Nsecfilt + k] = bcatalog->secfilt[i*Nsecfilt + k];
    }      

    catalogs->catalog[Nc].Naverage ++;
    if (catalogs->catalog[Nc].Naverage >= catalogs->NAVERAGE[Nc]) {
      catalogs->NAVERAGE[Nc] += MAX(catalogs->NAVERAGE[Nc], 1024);
      REALLOCATE (catalogs->catalog[Nc].averageT, AverageTiny, catalogs->NAVERAGE[Nc]);
      REALLOCATE (catalogs->catalog[Nc].secfilt,  SecFilt,     catalogs->NAVERAGE[Nc]*Nsecfilt);
    }	       
  }

  // assign the measures to the corresponding catalog
  // XXX what about averef and related links?  Do I need them?
  for (i = 0; i < bcatalog->Nmeasure; i++) {
    int ID = bcatalog->measure[i].catID;
    int Nc = catalogs->index[ID];
    assert (Nc > -1);
    assert (Nc < catalogs->NCATALOG);
    int Na = catalogs->catalog[Nc].Nmeasure;
    catalogs->catalog[Nc].measureT[Na] = bcatalog->measure[i];
    catalogs->catalog[Nc].Nmeasure ++;
    if (catalogs->catalog[Nc].Nmeasure >= catalogs->NMEASURE[Nc]) {
      catalogs->NMEASURE[Nc] += MAX(catalogs->NMEASURE[Nc], 4096);
      REALLOCATE (catalogs->catalog[Nc].measureT, MeasureTiny, catalogs->NMEASURE[Nc]);
    }	       
  }

  return TRUE;
}
