# include "addstar.h"
int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2);

# define NFILTER 5
# define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
# define dSIN(A)   ((double) sin ((double)RAD_DEG*A))

# define GET_COLUMN_5(NAME,TYPE) \
  TYPE *NAME; \
  NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \
  assert (NAME); assert (!strcmp (type, #TYPE)); assert (Nrow == Nstars); assert (Ncol == NFILTER);

# define GET_COLUMN_1(NAME,TYPE) \
  TYPE *NAME; \
  NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \
  assert (NAME); assert (!strcmp (type, #TYPE)); assert (Nrow == Nstars); assert (Ncol == 1);

/* grab named photcode */
# define NAMED_PHOTCODE_AND_ZP(CODE,ZP,NAME) { \
  PhotCode *code; \
  code = GetPhotcodebyName (NAME); \
  if (code == NULL) { \
    fprintf (stderr, "ERROR:  photcode %s not found in photcode table\n", NAME); \
    exit (0); } \
  CODE = code[0].code; \
  ZP = 0.001*code[0].C; }

// XXX NOTE : as of 2008.02.27, the zero point is still carried internally in millimags

// 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 *ReadStarsSDSS (FILE *f, char *name, Header *header, Header *in_theader, Image *images, int *nimages, unsigned int *nstars) {

  int i, j, N, Nskip, Nstars, camcol;
  char type[80];
  Header theader;
  FTable table;
  Stars *stars;
  double clockRate, mjd[5], jd, sidtime, alt, az;
  float seeing[5], photErr[5], zeropt[5], ZeroPt;
  time_t tzero[5];
  char filtname[16][5];
  int photcode[5];
  int Nrow, Ncol; // used in the GET_COLUMN_1,5 macros above
  
  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);
  }

  strcpy (filtname[0], "u");
  strcpy (filtname[1], "g");
  strcpy (filtname[2], "r");
  strcpy (filtname[3], "i");
  strcpy (filtname[4], "z");

  NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS");
  NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS");
  NAMED_PHOTCODE_AND_ZP (photcode[2], zeropt[2], "R_SDSS");
  NAMED_PHOTCODE_AND_ZP (photcode[3], zeropt[3], "I_SDSS");
  NAMED_PHOTCODE_AND_ZP (photcode[4], zeropt[4], "Z_SDSS");

  // various header values needed to calculate per-star data below
  gfits_scan (header, "C_OBS", "%lf", 1, &clockRate); // value in header is usec / unbinned row
  clockRate *= 1e-6; // convert to seconds / unbinned row

  gfits_scan (table.header, "MJD_U", "%lf", 1, &mjd[0]);
  gfits_scan (table.header, "MJD_G", "%lf", 1, &mjd[1]);
  gfits_scan (table.header, "MJD_R", "%lf", 1, &mjd[2]);
  gfits_scan (table.header, "MJD_I", "%lf", 1, &mjd[3]);
  gfits_scan (table.header, "MJD_Z", "%lf", 1, &mjd[4]);
  tzero[0] = ohana_mjd_to_sec (mjd[0]);
  tzero[1] = ohana_mjd_to_sec (mjd[1]);
  tzero[2] = ohana_mjd_to_sec (mjd[2]);
  tzero[3] = ohana_mjd_to_sec (mjd[3]);
  tzero[4] = ohana_mjd_to_sec (mjd[4]);

  gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]);
  gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]);
  gfits_scan (table.header, "SEEING_R", "%f", 1, &seeing[2]);
  gfits_scan (table.header, "SEEING_I", "%f", 1, &seeing[3]);
  gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]);

  gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]);
  gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]);
  gfits_scan (table.header, "PSFERR_R", "%f", 1, &photErr[2]);
  gfits_scan (table.header, "PSFERR_I", "%f", 1, &photErr[3]);
  gfits_scan (table.header, "PSFERR_Z", "%f", 1, &photErr[4]);

  gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row

  ZeroPt = GetZeroPoint();

  // create a Star entry for each filter and detection
  Nstars = table.header[0].Naxis[1];
  ALLOCATE (stars, Stars, NFILTER*Nstars);

  GET_COLUMN_5 (rowc, float);
  GET_COLUMN_5 (colc, float);
  GET_COLUMN_5 (sky, float);
  GET_COLUMN_5 (psfCounts, float);
  GET_COLUMN_5 (fiberCounts, float);
  GET_COLUMN_5 (offsetRa, float);
  GET_COLUMN_5 (offsetDec, float);
  GET_COLUMN_5 (flags, int);
  GET_COLUMN_5 (flags2, int);

  GET_COLUMN_5 (prob_psf, float);

  GET_COLUMN_1 (ra, double);
  GET_COLUMN_1 (dec, double);

  GET_COLUMN_5 (rowcErr, float);
  GET_COLUMN_5 (colcErr, float);
  GET_COLUMN_5 (skyErr, float);
  GET_COLUMN_5 (psfCountsErr, float);

  // the value of stars[].M is supposed to be the instrumental magnitude offset by the
  // default zero point 25.0 (-2.5*log_10(counts/sec) + ZeroPt).  The magnitude reported
  // by SDSS is the calibrated mag: -2.5*log_10(counts/sec) + C_0.  Adjust magnitudes to
  // compensate for the difference.

  for (i = 0; i < Nstars; i++) {
    for (j = 0; j < NFILTER; j++) {
      N = NFILTER*i + j;
      stars[N].X         = colc[N];
      stars[N].Y         = rowc[N];
      stars[N].dX        = colcErr[N];
      stars[N].dY        = rowcErr[N];
      stars[N].M         = psfCounts[N] + ZeroPt - zeropt[j];
      stars[N].dM        = psfCountsErr[N];
      stars[N].Map       = fiberCounts[N] + ZeroPt - zeropt[j];
      stars[N].Mpeak     = NAN;
      stars[N].sky       = sky[N]; // adjust this to counts?
      stars[N].dsky      = skyErr[N];
      stars[N].fx        = seeing[j]; // reported in arcsec?
      stars[N].fy        = seeing[j];
      stars[N].df        = 0.0;
      stars[N].psfChisq  = prob_psf[N]; // XXX not really the correct value...
      stars[N].psfQual   = 0.0;
      stars[N].crNsigma  = 0.0;
      stars[N].extNsigma = 0.0;
      stars[N].detID     = N;
      stars[N].R         = ra[i] + dCOS(dec[i]) * offsetRa[N] / 3600.0;
      stars[N].D         = dec[i] + offsetDec[N] / 3600.0;
      stars[N].dR        = NAN;
      stars[N].dD        = NAN;
      stars[N].uR        = 0.0; // I don't trust the SDSS proper motions
      stars[N].uD        = 0.0; // I don't trust the SDSS proper motions
      stars[N].duR       = 0.0;
      stars[N].duD       = 0.0;
      stars[N].P         = 0.0;
      stars[N].dP        = 0.0;
      stars[N].Mcal      = 0.0;
      stars[N].t         = tzero[j] + clockRate*rowc[N]; // time since row 0
      stars[N].dt        = 53.907456; // is this 2048*clockRate ?

      SetSDSSFlags (&stars[N], flags[N], flags2[N]);

      // longitude and latitude for SDSS:
      // Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W
      // longitude = 105.820419312 deg = 7.05469417
      // latitude = 32.7803611755 deg

      double Longitude = 7.05469417;   // hours (+ = W)
      double Latitude = 32.7803611755; // degrees

      jd = ohana_sec_to_jd (stars[N].t);
      sidtime  = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees
      altaz (&alt, &az, sidtime - stars[N].R, stars[N].D, Latitude);

      stars[N].airmass   = 1.0 / dCOS(90.0 - alt);
      stars[N].az        = az;

      stars[N].code      = photcode[j];
      stars[N].found     = -1;
      stars[N].dophot    = 0;
    }
  }    

  for (i = 0; i < NFILTER; i++) {

    N = i + *nimages;
    
    // XXX for now, we define a totally fake coordinate system centered on the first listed star
    strcpy (images[N].coords.ctype, "RA---TAN");
    
    images[N].coords.crval1 = stars[0].R;
    images[N].coords.crval2 = stars[0].D;
    
    images[N].coords.crpix1 = stars[0].X;
    images[N].coords.crpix2 = stars[0].Y;
    images[N].coords.cdelt1 = images[N].coords.cdelt2 = 0.4 / 3600.0;

    images[N].coords.pc1_1 = 0.0;
    images[N].coords.pc1_2 = 1.0;
    images[N].coords.pc2_1 = 1.0;
    images[N].coords.pc2_2 = 0.0;

    images[N].coords.Npolyterms = 0;
    memset (images[N].coords.polyterms, 0, 2*7*sizeof(float));

    images[N].NX = 2048;
    images[N].NY = 1490;

    images[N].tzero = tzero[i];
    images[N].cerror = 0.0;
 
    // set photcodes for the 5 images (SDSS_U,G,R,I,Z)
    images[N].photcode = photcode[i];

    // calculate this from : C_OBS, TRACKING, and NY
    images[N].exptime = 2048*clockRate;
  
    images[N].apmifit = 0.0;
    images[N].dapmifit = 0.0;
    images[N].detection_limit = 0.0; 
    images[N].saturation_limit = 0.0;
    images[N].fwhm_x = seeing[i];
    images[N].fwhm_y = seeing[i];

    // XXX longitude and latitude are known for SDSS
    // SDSS is at : Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W
    // longitude = 105.820419312 deg = 7.05469417
    // latitude = 32.7803611755 deg

    double Longitude = 7.05469417;   // hours (+ = W)
    double Latitude = 32.7803611755; // degrees

    jd = ohana_sec_to_jd (images[N].tzero);
    images[N].sidtime  = ohana_lst (jd, Longitude);
    images[N].latitude = Latitude;

    altaz (&alt, &az, 15.0*images[N].sidtime - images[N].coords.crval1, images[N].coords.crval2, Latitude);

    images[N].trate = clockRate * 1e-4;
    images[N].secz = stars[0].airmass;
    images[N].ccdnum = camcol;

    // secz is in units milli-airmass
    images[N].Mcal = 0.0;
    images[N].Xm   = NAN_S_SHORT;
    images[N].code = 0;

    images[N].nstar = Nstars;
  
    images[N].imageID = N;
    images[N].externID = 0;
    images[N].sourceID = 0;

    // save the filename
    snprintf (images[N].name, 64, "%s[%s]", name, filtname[i]);
  }

  *nimages += NFILTER;
  *nstars = Nstars*NFILTER;
  return (stars);
}

// ha/dec -> alt/az
int altaz (double *alt, double *az, double ha, double dec, double latitude) {

  double sind, sinh, cosh;

  sind = dSIN (dec) * dSIN (latitude) + dCOS (dec) * dCOS (ha) * dCOS (latitude);
  *alt  = DEG_RAD * asin (sind);

  sinh = - dCOS (dec) * dSIN (ha);
  cosh =   dSIN (dec) * dCOS (latitude) - dCOS (dec) * dCOS (ha) * dSIN (latitude);

  *az = DEG_RAD * atan2 (sinh, cosh);

  // I may need to use the parallactic angle to get the right plate rotation...
  // sinh = -dCOS(az) * dSIN(alt) * dSIN(ha) * dSIN(latitude) + dSIN(az) * dSIN(alt) * dCOS(ha) - dSIN(ha) * dCOS(alt) * dCOS(latitude);
  // cosh = -dSIN(az) * dSIN(ha) * dSIN(latitude) - dCOS(az) * dCOS(ha);
  // rot = -DEG_RAD * atan2 (sinh, cosh);
  return TRUE;
}

int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) {

  // XXX this is wrong, need to roll left to set the correct bit 
  if (flags1 & 0x00000002) star[0].flags |= 0x0001; // BRIGHT            - 1  1
  if (flags1 & 0x00000004) star[0].flags |= 0x0002; // EDGE              - 1  2
  if (flags1 & 0x00000008) star[0].flags |= 0x0004; // BLENDED           - 1  3
  if (flags1 & 0x00000010) star[0].flags |= 0x0008; // CHILD             - 1  4
  if (flags1 & 0x00000020) star[0].flags |= 0x0010; // PEAKCENTER        - 1  5
  if (flags1 & 0x00000040) star[0].flags |= 0x0020; // NODEBLEND         - 1  6
  if (flags1 & 0x00040000) star[0].flags |= 0x0040; // SATUR             - 1 18
  if (flags1 & 0x00080000) star[0].flags |= 0x0080; // NOTCHECKED        - 1 19
  if (flags1 & 0x10000000) star[0].flags |= 0x0100; // BINNED1           - 1 28
  if (flags1 & 0x20000000) star[0].flags |= 0x0200; // BINNED2           - 1 29
  if (flags1 & 0x40000000) star[0].flags |= 0x0400; // BINNED4           - 1 30
  if (flags2 & 0x00000040) star[0].flags |= 0x0800; // LOCAL_EDGE        - 2  7
  if (flags2 & 0x00000800) star[0].flags |= 0x1000; // INTERP_CENTER     - 2 12
  if (flags2 & 0x00002000) star[0].flags |= 0x2000; // DEBLEND_NOPEAK    - 2 14
  if (flags2 & 0x02000000) star[0].flags |= 0x4000; // NOTCHECKED_CENTER - 2 26
  return (TRUE);

}
