# 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) {

  off_t Nskip;
  unsigned int 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].datasize;
    fseeko (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 (!strcmp (type, "PS1_V1")) {
    stars = Convert_PS1_V1 (&table, &Nstars);
  }
  if (!strcmp (type, "PS1_V2")) {
    stars = Convert_PS1_V2 (&table, &Nstars);
  }
  if (stars == NULL) {
    fprintf (stderr, "ERROR: invalid table type %s\n", type);
    exit (1);
  }
  // Nstars is not necessarily == *nstars (The former is the number of detections, the
  // latter are the 'good' detections reported by the photometry system.
  *nstars = Nstars;

  return stars;
}

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

  off_t Nstars;
  unsigned int i;
  char swapped;
  double ZeroPt;
  Stars *stars = NULL;
  SMPData *smpdata = NULL;

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

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    InitStar (&stars[i]);

    stars[i].measure.Xccd      = smpdata[i].X;
    stars[i].measure.Yccd      = smpdata[i].Y;

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

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

    // the dophot type information gets pushed into the upper 2 bytes of photFlags
    stars[i].measure.photFlags = (smpdata[i].dophot << 16);

    stars[i].measure.FWx       = ShortPixels (smpdata[i].fx);
    stars[i].measure.FWy       = ShortPixels (smpdata[i].fy);
    stars[i].measure.theta     = ShortDegree (smpdata[i].df);
  }    
  *nstars = Nstars;
  return (stars);
}

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

  off_t Nstars;
  unsigned int i;
  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++) {
    InitStar (&stars[i]);
    stars[i].measure.Xccd     = ps1data[i].X;
    stars[i].measure.Yccd     = ps1data[i].Y;

    stars[i].measure.dXccd    = ShortPixels(ps1data[i].dX);
    stars[i].measure.dYccd    = ShortPixels(ps1data[i].dY);
   
    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
      stars[i].measure.M      = NAN;
    } else {
      stars[i].measure.M      = ps1data[i].M + ZeroPt;
    }
    stars[i].measure.Map      = NAN;
    stars[i].measure.dM       = ps1data[i].dM;
    stars[i].measure.Sky      = ps1data[i].sky;
    stars[i].measure.dSky     = ps1data[i].dSky;

    stars[i].measure.FWx      = ShortPixels(ps1data[i].fx);
    stars[i].measure.FWy      = ShortPixels(ps1data[i].fy);
    stars[i].measure.theta    = ShortDegree(ps1data[i].df);

    stars[i].measure.psfChisq = ps1data[i].psfChisq;
    stars[i].measure.psfQual  = ps1data[i].psfQual;

    stars[i].measure.detID    = ps1data[i].detID;
  }    
  *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, unsigned int *nstars) {

  off_t Nstars; 
  unsigned int i;
  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++) {
    InitStar (&stars[i]);
    stars[i].measure.Xccd       = ps1data[i].X;
    stars[i].measure.Yccd       = ps1data[i].Y;

    stars[i].measure.dXccd      = ShortPixels(ps1data[i].dX);
    stars[i].measure.dYccd      = ShortPixels(ps1data[i].dY);

    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
	stars[i].measure.M      = NAN;
    } else {
	stars[i].measure.M      = ps1data[i].M + ZeroPt;
    }
    stars[i].measure.Map        = NAN;
    stars[i].measure.dM         = ps1data[i].dM;
    stars[i].measure.Sky        = ps1data[i].sky;
    stars[i].measure.dSky       = ps1data[i].dSky;

    stars[i].measure.FWx        = ShortPixels(ps1data[i].fx);
    stars[i].measure.FWy        = ShortPixels(ps1data[i].fy);
    stars[i].measure.theta      = ShortDegree(ps1data[i].df);

    stars[i].measure.psfChisq  	= ps1data[i].psfChisq;
    stars[i].measure.psfQual   	= ps1data[i].psfQual;
    stars[i].measure.crNsigma  	= ps1data[i].crNsigma;
    stars[i].measure.extNsigma 	= ps1data[i].extNsigma;

    stars[i].measure.detID     	= ps1data[i].detID;
    stars[i].measure.photFlags  = ps1data[i].flags;
  }    
  *nstars = Nstars;
  return (stars);
}

Stars *Convert_PS1_V1 (FTable *table, unsigned int *nstars) {

  off_t Nstars; 
  unsigned int i;
  double ZeroPt;
  Stars *stars;
  CMF_PS1_V1 *ps1data;

  // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec.  this was a
  // mistake in two ways: a few format should have been defined (eg, CMF_PS1_V2), and the
  // layout used did not have clean byte-boundaries for the corresponding structure.  The
  // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the
  // autocode tools do not work to read in the new version, even if we recognize it.  Here we
  // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and
  // call a special conversion function if it is found.

  if (table[0].header[0].Naxis[0] == 136) {
      stars = Convert_PS1_V1_Alt (table, nstars);
      return (stars);
  }

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

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    InitStar (&stars[i]);
    stars[i].measure.Xccd       = ps1data[i].X;
    stars[i].measure.Yccd       = ps1data[i].Y;
    stars[i].measure.dXccd      = ShortPixels(ps1data[i].dX);
    stars[i].measure.dYccd      = ShortPixels(ps1data[i].dY);

    stars[i].measure.posangle   = ShortDegree(ps1data[i].posangle);
    stars[i].measure.pltscale   = ps1data[i].pltscale;

    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
	stars[i].measure.M      = NAN;
    } else {
	stars[i].measure.M      = ps1data[i].M + ZeroPt;
    }
    stars[i].measure.dM         = ps1data[i].dM;
    stars[i].measure.dMcal      = ps1data[i].dMcal;
    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
		        
    stars[i].measure.Sky        = ps1data[i].sky;
    stars[i].measure.dSky       = ps1data[i].dSky;
		        
    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
    stars[i].measure.psfQual    = ps1data[i].psfQual;
    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
    stars[i].measure.extNsigma  = ps1data[i].extNsigma;

    stars[i].measure.FWx        = ShortPixels(ps1data[i].fx);
    stars[i].measure.FWy        = ShortPixels(ps1data[i].fy);
    stars[i].measure.theta      = ShortDegree(ps1data[i].df);

    stars[i].measure.Mxx        = ShortPixels(ps1data[i].Mxx);
    stars[i].measure.Mxy        = ShortPixels(ps1data[i].Mxy);
    stars[i].measure.Myy        = ShortPixels(ps1data[i].Myy);
		        
    stars[i].measure.photFlags  = ps1data[i].flags;

    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
    stars[i].measure.detID      = ps1data[i].detID; 

    // the Average fields and the following Measure fields are set in FilterStars after
    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID, 
    // averef is set in find_matches, dbFlags is zero on ingest.

    // the following fields are currently not being set anywhere: t_msec
  }    
  *nstars = Nstars;
  return (stars);
}

Stars *Convert_PS1_V1_Alt (FTable *table, unsigned int *nstars) {

  off_t Nstars; 
  unsigned int i;
  double ZeroPt;
  Stars *stars;
  CMF_PS1_V2 *ps1data;

  // CMF_PS1_V1 was modified 2009.05.26 (r24251) to use doubles for ra & dec.  this was a
  // mistake in two ways: a few format should have been defined (eg, CMF_PS1_V2), and the
  // layout used did not have clean byte-boundaries for the corresponding structure.  The
  // former means we have two varieties of CMF_PS1_V1 out there; the latter means that the
  // autocode tools do not work to read in the new version, even if we recognize it.  Here we
  // test for the existence of the broken version (table[0].headers[0].Naxis[0] == 136), and
  // call a special conversion function if it is found.

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

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    InitStar (&stars[i]);
    stars[i].measure.Xccd       = ps1data[i].X;
    stars[i].measure.Yccd       = ps1data[i].Y;
    stars[i].measure.dXccd      = ShortPixels(ps1data[i].dX);
    stars[i].measure.dYccd      = ShortPixels(ps1data[i].dY);

    stars[i].measure.posangle   = ShortDegree(ps1data[i].posangle);
    stars[i].measure.pltscale   = ps1data[i].pltscale;

    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
	stars[i].measure.M      = NAN;
    } else {
	stars[i].measure.M      = ps1data[i].M + ZeroPt;
    }
    stars[i].measure.dM         = ps1data[i].dM;
    stars[i].measure.dMcal      = ps1data[i].dMcal;
    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
		        
    stars[i].measure.Sky        = ps1data[i].sky;
    stars[i].measure.dSky       = ps1data[i].dSky;
		        
    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
    stars[i].measure.psfQual    = ps1data[i].psfQual;
    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
    stars[i].measure.extNsigma  = ps1data[i].extNsigma;

    stars[i].measure.FWx        = ShortPixels(ps1data[i].fx);
    stars[i].measure.FWy        = ShortPixels(ps1data[i].fy);
    stars[i].measure.theta      = ShortDegree(ps1data[i].df);

    stars[i].measure.Mxx        = ShortPixels(ps1data[i].Mxx);
    stars[i].measure.Mxy        = ShortPixels(ps1data[i].Mxy);
    stars[i].measure.Myy        = ShortPixels(ps1data[i].Myy);
		        
    stars[i].measure.photFlags  = ps1data[i].flags;

    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
    stars[i].measure.detID      = ps1data[i].detID; 

    // the Average fields and the following Measure fields are set in FilterStars after
    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID, 
    // averef is set in find_matches, dbFlags is zero on ingest.

    // the following fields are currently not being set anywhere: t_msec
  }    
  *nstars = Nstars;
  return (stars);
}

Stars *Convert_PS1_V2 (FTable *table, unsigned int *nstars) {

  off_t Nstars; 
  unsigned int i;
  double ZeroPt;
  Stars *stars;
  CMF_PS1_V2 *ps1data;

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

  ALLOCATE (stars, Stars, Nstars);
  for (i = 0; i < Nstars; i++) {
    InitStar (&stars[i]);
    stars[i].measure.Xccd       = ps1data[i].X;
    stars[i].measure.Yccd       = ps1data[i].Y;
    stars[i].measure.dXccd      = ShortPixels(ps1data[i].dX);
    stars[i].measure.dYccd      = ShortPixels(ps1data[i].dY);

    stars[i].measure.posangle   = ShortDegree(ps1data[i].posangle);
    stars[i].measure.pltscale   = ps1data[i].pltscale;

    if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) {
	stars[i].measure.M      = NAN;
    } else {
	stars[i].measure.M      = ps1data[i].M + ZeroPt;
    }
    stars[i].measure.dM         = ps1data[i].dM;
    stars[i].measure.dMcal      = ps1data[i].dMcal;
    stars[i].measure.Map        = ps1data[i].Map + ZeroPt;
		        
    stars[i].measure.Sky        = ps1data[i].sky;
    stars[i].measure.dSky       = ps1data[i].dSky;
		        
    stars[i].measure.psfChisq   = ps1data[i].psfChisq;
    stars[i].measure.psfQual    = ps1data[i].psfQual;
    stars[i].measure.psfNdof    = ps1data[i].psfNdof;
    stars[i].measure.psfNpix    = ps1data[i].psfNpix;
    stars[i].measure.crNsigma   = ps1data[i].crNsigma;
    stars[i].measure.extNsigma  = ps1data[i].extNsigma;

    stars[i].measure.FWx        = ShortPixels(ps1data[i].fx);
    stars[i].measure.FWy        = ShortPixels(ps1data[i].fy);
    stars[i].measure.theta      = ShortDegree(ps1data[i].df);

    stars[i].measure.Mxx        = ShortPixels(ps1data[i].Mxx);
    stars[i].measure.Mxy        = ShortPixels(ps1data[i].Mxy);
    stars[i].measure.Myy        = ShortPixels(ps1data[i].Myy);
		        
    stars[i].measure.photFlags  = ps1data[i].flags;

    // this is may optionally be replaced by the internal sequence (see FilterStars.c)
    stars[i].measure.detID      = ps1data[i].detID; 

    // the Average fields and the following Measure fields are set in FilterStars after
    // the image metadata is in hand:  dR, dD, Mcal, dt, airmass, az, t, imageID, extID, 
    // averef is set in find_matches, dbFlags is zero on ingest.

    // the following fields are currently not being set anywhere: t_msec
  }    
  *nstars = Nstars;
  return (stars);
}
