# include "addstar.h"

// given a file with the pointer at the start of the table block and the 
// corresponding image header, load the stars from the table
Stars *ReadStarsFITS (FILE *f, Header *header, Header *in_theader, unsigned int *nstars) {

  int Nskip, Nstars;
  char type[80];
  Header theader;
  FTable table;
  Stars *stars;
  
  if (in_theader == NULL) {
    table.header = &theader;
    if (!gfits_fread_header (f, table.header)) Shutdown ("ERROR: can't read table header");
  } else {
    table.header = in_theader;
    Nskip = in_theader[0].size;
    fseek (f, Nskip, SEEK_CUR); 
  }

  /* load the table data */
  if (!gfits_fread_ftable_data (f, &table)) {
    fprintf (stderr, "ERROR: can't read table header\n");
    exit (1);
  }

  if (!gfits_scan (table.header, "EXTTYPE", "%s", 1, type)) {
    strcpy (type, "SMPDATA");
  }

  stars = NULL;
  if (!strcmp (type, "SMPDATA")) {
    stars = Convert_SMPDATA (&table, &Nstars);
  }
  if (!strcmp (type, "PS1_DEV_0")) {
    stars = Convert_PS1_DEV_0 (&table, &Nstars);
  }
  if (!strcmp (type, "PS1_DEV_1")) {
    stars = Convert_PS1_DEV_1 (&table, &Nstars);
  }
  if (stars == NULL) {
    fprintf (stderr, "ERROR: invalid table type %s\n", type);
    exit (1);
  }
  if (*nstars != Nstars) {
    fprintf (stderr, "WARNING: inconsistent number of stars? %d vs %d\n", *nstars, Nstars);
  }
  *nstars = Nstars;
  return stars;
}

Stars *Convert_SMPDATA (FTable *table, int *nstars) {

  int i, Nstars, swapped;
  double ZeroPt;
  Stars *stars = NULL;
  SMPData *smpdata = NULL;

  swapped = FALSE;
  smpdata = gfits_table_get_SMPData (table, &Nstars, &swapped);
  ZeroPt = GetZeroPoint();

  /* XXX we need to check at least the size of the loaded table */
  // XXX use memset to zero-out the elements before setting

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].X      = smpdata[i].X;
    stars[i].Y      = smpdata[i].Y;

    if ((smpdata[i].M >= ZeroPt) || isnan(smpdata[i].M)) {
      stars[i].M    = NAN;
      stars[i].Mgal = NAN;
      stars[i].Map  = NAN;
    } else {
      stars[i].M    = smpdata[i].M;
      stars[i].Mgal = smpdata[i].M;
      stars[i].Map  = smpdata[i].M;
    }

    stars[i].dM     = smpdata[i].dM*0.001;
    stars[i].dophot = smpdata[i].dophot;

    stars[i].fx     = smpdata[i].fx;
    stars[i].fy     = smpdata[i].fy;
    stars[i].df     = smpdata[i].df;
  }    
  *nstars = Nstars;
  return (stars);
}

Stars *Convert_PS1_DEV_0 (FTable *table, int *nstars) {

  int i, Nstars;
  double ZeroPt;
  Stars *stars;
  PS1_DEV_0 *ps1data;

  ps1data = gfits_table_get_PS1_DEV_0 (table, &Nstars, NULL);
  ZeroPt = GetZeroPoint();

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].X       = ps1data[i].X;
    stars[i].Y       = ps1data[i].Y;
    stars[i].dX      = ps1data[i].dX;
    stars[i].dY      = ps1data[i].dY;
    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
      stars[i].M     = NAN;
    } else {
      stars[i].M     = ps1data[i].M + ZeroPt;
    }
    stars[i].dM      = ps1data[i].dM;
    stars[i].Mpeak   = ps1data[i].Mpeak;

    stars[i].sky     = ps1data[i].sky;
    stars[i].dsky    = ps1data[i].dSky;

    stars[i].fx      = ps1data[i].fx;
    stars[i].fy      = ps1data[i].fy;
    stars[i].df      = ps1data[i].df;

    stars[i].psfProb = ps1data[i].psfProb;
    stars[i].psfQual = ps1data[i].psfQual;

    stars[i].detID   = ps1data[i].detID;

    /* these are set elsewhere */
    stars[i].R       = 0.0;
    stars[i].D       = 0.0;
    stars[i].dR      = 0.0;
    stars[i].dD      = 0.0;

    stars[i].uR      = 0.0;
    stars[i].uD      = 0.0;
    stars[i].duR     = 0.0;
    stars[i].duD     = 0.0;

    stars[i].P       = 0.0;
    stars[i].dP      = 0.0;

    stars[i].Mcal    = 0;
    stars[i].t       = 0;
    stars[i].dt      = 0;
    stars[i].airmass = 0;
    stars[i].code    = 0;
    stars[i].found   = 0;

    /* these are not used */
    stars[i].Map     = NAN;
    stars[i].Mgal    = NAN;
    stars[i].dophot  = 0;
  }    
  *nstars = Nstars;
  return (stars);
}

// XXX I need to make the IPP I/O functions and these functions
// consistent wrt ZERO_POINT....
Stars *Convert_PS1_DEV_1 (FTable *table, int *nstars) {

  int i, Nstars;
  double ZeroPt;
  Stars *stars;
  PS1_DEV_1 *ps1data;

  ps1data = gfits_table_get_PS1_DEV_1 (table, &Nstars, NULL);
  ZeroPt = GetZeroPoint();

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].X       = ps1data[i].X;
    stars[i].Y       = ps1data[i].Y;
    stars[i].dX      = ps1data[i].dX;
    stars[i].dY      = ps1data[i].dY;
    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
	stars[i].M   = NAN;
    } else {
	stars[i].M   = ps1data[i].M + ZeroPt;
    }
    stars[i].dM      = ps1data[i].dM;
    stars[i].Mpeak   = ps1data[i].Mpeak;

    stars[i].sky     = ps1data[i].sky;
    stars[i].dsky    = ps1data[i].dSky;

    stars[i].fx      = ps1data[i].fx;
    stars[i].fy      = ps1data[i].fy;
    stars[i].df      = ps1data[i].df;

    stars[i].psfProb   = ps1data[i].psfProb;
    stars[i].psfQual   = ps1data[i].psfQual;
    stars[i].crNsigma  = ps1data[i].crNsigma;
    stars[i].extNsigma = ps1data[i].extNsigma;

    stars[i].detID     = ps1data[i].detID;
    stars[i].flags     = ps1data[i].flags;

    // XXX not defined anyway
    // stars[i].stargal   = ps1data[i].stargal;

    /* these are set elsewhere */
    stars[i].R       = 0.0;
    stars[i].D       = 0.0;
    stars[i].dR      = 0.0;
    stars[i].dD      = 0.0;

    stars[i].uR      = 0.0;
    stars[i].uD      = 0.0;
    stars[i].duR     = 0.0;
    stars[i].duD     = 0.0;

    stars[i].P       = 0.0;
    stars[i].dP      = 0.0;

    stars[i].Mcal    = 0;
    stars[i].t       = 0;
    stars[i].dt      = 0;
    stars[i].airmass = 0;
    stars[i].az      = 0;
    stars[i].code    = 0;
    stars[i].found   = 0;

    /* these are not used */
    stars[i].Map     = NAN;
    stars[i].Mgal    = NAN;
    stars[i].dophot  = 0;
  }    
  *nstars = Nstars;
  return (stars);
}
