# include "mkcmf.h"

# define ZERO_POINT 25.0
# define SKY 100.0
# define DSKY 2.0
# define PSFCHI 1.3
# define CRN 3.0
# define EXTN 4.0
# define FX 1.5
# define FY 1.0
# define DF 10.0
# define PSFQUAL 0.98
# define FLAGS 0x1101

void gauss_init (int Nbin);
double rnd_gauss (double mean, double sigma);
void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_V4 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_V3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
void writeStars_PS1_SV3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
void writeStars_PS1_DEV_1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars);

int ADDNOISE = TRUE;
float BAD_PSFQF_FRAC = 0.0;

static float exptime = 1.0;
static Coords coords;

static char reserved[] =  "Reserved space.  This line can be used to add a new FITS card.";

int main (int argc, char **argv) {

  // generate a simple, fake cmf file
  // load a text table with X,Y,Mag (instrumental?)

  int N, Nstars, NSTARS, found;
  unsigned int *Flag;
  double *X, *Y, *M, Xmax, Ymax;

  FILE *f, *fits;
  Header header;
  Header theader;
  Matrix matrix;
  FTable ftable;

  int APPEND = FALSE;
  if ((N = get_argument (argc, argv, "-append"))) {
    remove_argument (N, &argc, argv);
    APPEND = TRUE;
  }

  static char *photcode = "SIMTEST.r.Chip";
  if ((N = get_argument (argc, argv, "-photcode"))) {
    remove_argument (N, &argc, argv);
    photcode = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  static double mjd;
  mjd = NAN;
  if ((N = get_argument (argc, argv, "-mjd"))) {
    remove_argument (N, &argc, argv);
    mjd = strtod (argv[N], NULL);
    remove_argument (N, &argc, argv);
  }

  static char *date = "2001-01-01";
  if ((N = get_argument (argc, argv, "-date"))) {
    remove_argument (N, &argc, argv);
    date = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  static char *time = "00:00:00";
  if ((N = get_argument (argc, argv, "-time"))) {
    remove_argument (N, &argc, argv);
    time = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  char extname[80], exthead[80];
  strcpy (extname, "Chip.psf");
  strcpy (exthead, "Chip.hdr");
  if ((N = get_argument (argc, argv, "-extname"))) {
    remove_argument (N, &argc, argv);
    snprintf (extname, 80, "%s.psf", argv[N]);
    snprintf (exthead, 80, "%s.hdr", argv[N]);
    remove_argument (N, &argc, argv);
  }

  double RA = 10.0;
  double DEC = 20.0;
  if ((N = get_argument (argc, argv, "-radec"))) {
    remove_argument (N, &argc, argv);
    RA = atof (argv[N]);
    remove_argument (N, &argc, argv);
    DEC = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  double CRPIX1 = 0.0;
  double CRPIX2 = 0.0;
  if ((N = get_argument (argc, argv, "-crpix"))) {
    remove_argument (N, &argc, argv);
    CRPIX1 = atof (argv[N]);
    remove_argument (N, &argc, argv);
    CRPIX2 = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  if ((N = get_argument (argc, argv, "-exptime"))) {
    remove_argument (N, &argc, argv);
    exptime = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }  

  // XXX note that the airmass and ra,dec,mjd can be inconsistent (for a given observatory location)
  float airmass = 1.0;
  if ((N = get_argument (argc, argv, "-airmass"))) {
    remove_argument (N, &argc, argv);
    airmass = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }  

  int imageID = -1;
  if ((N = get_argument (argc, argv, "-imageID"))) {
    remove_argument (N, &argc, argv);
    imageID = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }  
  int sourceID = -1;
  if ((N = get_argument (argc, argv, "-sourceID"))) {
    remove_argument (N, &argc, argv);
    sourceID = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }  

  int NX = 0;
  int NY = 0;
  if ((N = get_argument (argc, argv, "-size"))) {
    remove_argument (N, &argc, argv);
    NX = atof (argv[N]);
    remove_argument (N, &argc, argv);
    NY = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  // add support for all cmf types
  static char *type = "PS1_V2";
  if ((N = get_argument (argc, argv, "-type"))) {
    remove_argument (N, &argc, argv);
    type = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  // add support for all cmf types
  int FROM_COORDS = FALSE;
  if ((N = get_argument (argc, argv, "-coords"))) {
    remove_argument (N, &argc, argv);
    FROM_COORDS = TRUE;
  }

  // expect a field with the photom flags for each detection
  int READ_FLAGS = FALSE;
  if ((N = get_argument (argc, argv, "-flags"))) {
    remove_argument (N, &argc, argv);
    READ_FLAGS = TRUE;
  }

  // add support for all cmf types
  ADDNOISE = TRUE;
  if ((N = get_argument (argc, argv, "-no-noise"))) {
    remove_argument (N, &argc, argv);
    ADDNOISE = FALSE;
  }

  // random bad PSF_QF values
  if ((N = get_argument (argc, argv, "-bad-psfqf-frac"))) {
    remove_argument (N, &argc, argv);
    BAD_PSFQF_FRAC = atof(argv[N]);
    remove_argument (N, &argc, argv);
  }

  if (argc != 3) {
    fprintf (stderr, "USAGE mkcmf (input) (output) [-date date] [-time time] [-radec ra dec] [-cmftype type]\n");
    exit (2);
  }

  // init the random seed
  { 
    struct timeval now;
    gettimeofday (&now, NULL);
    long A = now.tv_sec + now.tv_usec * 1000000;
    srand48(A);
  }
    
  /* bore site center guess */
  InitCoords (&coords, "DEC--TAN");
  coords.crval1 = RA;
  coords.crval2 = DEC;
  coords.crpix1 = CRPIX1;
  coords.crpix2 = CRPIX2;
  
  coords.cdelt1 = 0.25/3600.0;
  coords.cdelt2 = 0.25/3600.0;

  // load stars and generate complete output fields
  f = fopen (argv[1], "r");
  if (f == NULL) {
    fprintf (stderr, "unable to open input file %s\n", argv[1]);
    exit (1);
  }
    
  gauss_init (2048);

  // load test stars from a file:
  Nstars = 0;
  NSTARS = 100;
  ALLOCATE (X, double, NSTARS);
  ALLOCATE (Y, double, NSTARS);
  ALLOCATE (M, double, NSTARS);
  ALLOCATE (Flag, unsigned int, NSTARS);

  Xmax = Ymax = 0;
  while (1) {
    int status;
    double ra, dec, mag, xobs, yobs, xraw, yraw, mraw;
    unsigned int flags;
    if (FROM_COORDS) {
      status = fscanf (f, "%lf %lf %lf %lf %lf %lf", &ra, &dec, &mag, &xraw, &yraw, &mraw);
      RD_to_XY (&xobs, &yobs, ra, dec, &coords);
    } else {
      status = fscanf (f, "%lf %lf %lf", &xobs, &yobs, &mraw);
    }
    if (status == EOF) break;
    if (READ_FLAGS) {
      status = fscanf (f, "%x", &flags);
      if (status == EOF) {
	  fprintf (stderr, "error: missing flag for last star?\n");
	  exit (1);
      }
      Flag[Nstars] = flags;
    } else {
      Flag[Nstars] = 0;
    }

    X[Nstars] = xobs;
    Y[Nstars] = yobs;
    M[Nstars] = mraw;

    Xmax  = MAX(Xmax, X[Nstars]);
    Ymax  = MAX(Ymax, Y[Nstars]);

    if (Nstars == NSTARS - 1) {
      NSTARS += 100;
      REALLOCATE (X, double, NSTARS);
      REALLOCATE (Y, double, NSTARS);
      REALLOCATE (M, double, NSTARS);
      REALLOCATE (Flag, unsigned int, NSTARS);
    }
    Nstars ++;
  }
  if (NX && NY) {
      Xmax = NX;
      Ymax = NY;
  }

  // create primary header
  gfits_init_header (&header);    
  header.extend = TRUE;
  gfits_create_header (&header);
  gfits_create_matrix (&header, &matrix);
  gfits_print (&header, "NEXTEND", "%d", 1, 1);

  // XXX add minimum needed header fields 
  gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax));
  gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax));

  gfits_modify (&header, "NSTARS",   "%d", 1, Nstars);
  gfits_modify (&header, "PHOTCODE", "%s", 1, photcode);
  if (isfinite(mjd)) {
    gfits_modify (&header, "MJD-OBS", "%lf", 1, mjd);
  } else {
    gfits_modify (&header, "DATE-OBS", "%s", 1, date);
    gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
  }
  gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT);
  gfits_modify (&header, "EXPTIME", "%lf", 1, exptime);
  gfits_modify (&header, "AIRMASS", "%lf", 1, airmass);
  gfits_modify (&header, "NASTRO",   "%d", 1, 100); 

  if (imageID == -1) imageID = 100000.0*drand48();
  if (sourceID == -1) sourceID = 1000.0*drand48();
  gfits_modify (&header, "IMAGEID",  "%d", 1, imageID);
  gfits_modify (&header, "SOURCEID", "%d", 1, sourceID);

  PutCoords (&coords, &header);
  gfits_modify (&header, "EXTNAME",   "%s", 1, exthead);

  int i;
  for (i = 1; i < 32; i++) {
    gfits_modify_alt (&header, "COMMENT", "%C", i, reserved);
  }

  ftable.header = &theader;

  // set up desired CMF type:
  found = FALSE;
  if (!strcmp(type, "PS1_DEV_0")) {
    writeStars_PS1_DEV_0 (&ftable, X, Y, M, Nstars);
    found = TRUE;
  }
  if (!strcmp(type, "PS1_DEV_1")) {
    writeStars_PS1_DEV_1 (&ftable, X, Y, M, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V1")) {
    writeStars_PS1_V1 (&ftable, X, Y, M, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V2")) {
    writeStars_PS1_V2 (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V3")) {
    writeStars_PS1_V3 (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V4")) {
    writeStars_PS1_V4 (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V5")) {
    writeStars_PS1_V5 (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_V5_Lensing")) {
    writeStars_PS1_V5_Lensing (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!strcmp(type, "PS1_SV3")) {
    writeStars_PS1_SV3 (&ftable, X, Y, M, Flag, Nstars); 
    found = TRUE;
  }
  if (!found) {
    fprintf (stderr, "ERROR: unknown CMF type %s\n", type);
    exit (1);
  }

  gfits_modify (ftable.header, "EXTHEAD",   "%s", 1, exthead);
  gfits_modify (ftable.header, "EXTNAME",   "%s", 1, extname);

  // open for append (if you want a new file, need to blow the old object)
  if (APPEND) {
    fits = fopen (argv[2], "a+");
  } else {
    fits = fopen (argv[2], "w");
  }
  if (fits == NULL) {
    fprintf (stderr, "ERROR: can't open output file %s\n", argv[2]);
    exit (1);
  }

  /* if we are appending, fix up header */
  if (APPEND) {
    static char simple[] = "XTENSION= 'IMAGE  '            / Image extension";
    int Ns, No;
    Ns = strlen (simple);
    No = 80 - Ns;
    strncpy (header.buffer, simple, Ns);
    memset (&header.buffer[Ns], ' ', No);
  }

  gfits_fwrite_header  (fits, &header);
  gfits_fwrite_matrix  (fits, &matrix);
  gfits_fwrite_Theader (fits, ftable.header);
  gfits_fwrite_table   (fits, &ftable);
  fclose (fits);

  exit (0);
}


static int Ngaussint = 0;
static double *gaussint;

extern double drand48();

double gaussian (double x, double mean, double sigma) {

  double f;

  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));

  return (f);

}

/* integrate a gaussian from -5 sigma to +5 sigma */
void gauss_init (int Nbin) {
 
  int i;
  double val, x, dx, dx1, dx2, dx3, df;
  double mean, sigma;
 
  /* no need to generate this if it already exists */
  if (Ngaussint == Nbin) return;

  // A = time(NULL);
  // // XXX this is expensive if called a lot (1 sec min)
  // // for (B = 0; A == time(NULL); B++);
  // B = A + 10000;
  // srand48(B);
 
  Ngaussint = Nbin;
  ALLOCATE (gaussint, double, Ngaussint + 1);

  val = 0;
  dx = 1.0 / Ngaussint;
  dx1 = dx / 3.0;
  dx2 = 2.0*dx/3.0;
  dx3 = dx;
  mean = 0.0;
  sigma = 1.0;
 
  for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
    df = (3.0*gaussian(x    , mean, sigma) + 
          9.0*gaussian(x+dx1, mean, sigma) +
          9.0*gaussian(x+dx2, mean, sigma) + 
          3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    val += df;
    if (val > (i + 0.5) / (double) Ngaussint) {
      gaussint[i] = x + dx / 2.0;
      i++;
    }
  }
}

double rnd_gauss (double mean, double sigma) {
 
  int i;
  double y;
 
  y = drand48();
  i = Ngaussint*y;
  y = gaussint[i]*sigma + mean;
 
  return (y);
 
}
 
double int_gauss (int i) {
  double y;
  y = gaussint[i];
  return (y);
}
 
void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {

  int i;
  PS1_DEV_0 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, PS1_DEV_0, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].psfQF   = PSFQUAL;
    stars[i].nFrames   = 1;
  }

  gfits_table_set_PS1_DEV_0 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_DEV_0");
}

void writeStars_PS1_DEV_1 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {

  int i;
  PS1_DEV_1 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, PS1_DEV_1, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].psfQF   = PSFQUAL;
    stars[i].nFrames   = 1;
    stars[i].flags     = FLAGS;
  }

  gfits_table_set_PS1_DEV_1 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_DEV_1");
}

void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {

  int i;
  CMF_PS1_V1 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V1, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].psfQF   = PSFQUAL;
    stars[i].nFrames   = 1;
    stars[i].flags     = FLAGS;
  }

  gfits_table_set_CMF_PS1_V1 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V1");
}

void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_V2 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V2, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF   = 0.25;
    } else {
      stars[i].psfQF   = PSFQUAL;
    }
    
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].nFrames   = 1;
    stars[i].flags     = Flag[i];
  }

  gfits_table_set_CMF_PS1_V2 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V2");
}

void writeStars_PS1_V3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_V3 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V3, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];
    stars[i].Map = M[i] - 0.05;

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF   = 0.25;
    } else {
      stars[i].psfQF   = PSFQUAL;
    }
    
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].nFrames   = 1;
    stars[i].flags     = Flag[i];

    stars[i].kronFlux  = flux * 1.25;
    stars[i].kronFluxErr = fSN * flux * 1.25;
  }

  gfits_table_set_CMF_PS1_V3 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V3");
}

void writeStars_PS1_V4 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_V4 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V4, Nstars);
  for (i = 0; i < Nstars; i++) {
    stars[i].detID = i;

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].M = M[i];
    stars[i].Map = M[i] - 0.05;
    stars[i].MapRaw = M[i] - 0.10;

    stars[i].Flux = flux;
    stars[i].dFlux = flux * fSN;

    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
    stars[i].apFluxErr = stars[i].apFlux * fSN;

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
    }

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF     = 0.25;
      stars[i].psfQFperf = 0.24;
    } else {
      stars[i].psfQF     = PSFQUAL;
      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
    }
    
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;
    stars[i].dM = fSN;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;
    stars[i].nFrames   = 1;
    stars[i].flags     = Flag[i];

    stars[i].kronFlux  = flux * 1.25;
    stars[i].kronFluxErr = fSN * flux * 1.25;
  }

  gfits_table_set_CMF_PS1_V4 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V4");
}

void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_V5 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V5, Nstars);
  for (i = 0; i < Nstars; i++) {

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].detID = i;
    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;

    stars[i].posangle = 10.0;
    stars[i].pltscale = 0.25;

    stars[i].M = M[i];
    stars[i].dM = fSN;

    stars[i].Flux = flux;
    stars[i].dFlux = flux * fSN;

    stars[i].Map = M[i] - 0.05;
    stars[i].MapRaw = M[i] - 0.10;

    stars[i].apRadius = 8.0;

    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
    stars[i].apFluxErr = stars[i].apFlux * fSN;

    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
    stars[i].dMcal = 0.05;

    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
    stars[i].apNpix = 3.14*8.0*8.0;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;

    stars[i].k         = 1.0;
    stars[i].fwhmMaj   = FX*2.8;
    stars[i].fwhmMin   = FY*2.8;

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF     = 0.25;
      stars[i].psfQFperf = 0.24;
    } else {
      stars[i].psfQF     = PSFQUAL;
      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
    }

    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;

    stars[i].Mxx       = FX;
    stars[i].Mxy       = 0.01;
    stars[i].Myy       = FX;
    stars[i].M3c       = FX;
    stars[i].M3s       = FX;
    stars[i].M4c       = FX;
    stars[i].M4s       = FX;
    stars[i].Mr1       = FX;
    stars[i].Mrh       = FX;

    stars[i].kronFlux  = flux * 1.25;
    stars[i].kronFluxErr = fSN * flux * 1.25;

    stars[i].kronInner = fSN * flux * 0.9;
    stars[i].kronOuter = fSN * flux * 1.5;

    stars[i].skyLimitRad = 1;
    stars[i].skyLimitFlux = 2;
    stars[i].skyLimitSlope = 0.2;

    stars[i].flags     = Flag[i];
    stars[i].flags2    = 0x80;

    stars[i].nFrames   = 1;

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      float Moff = fSN*rnd_gauss(0.0, 1.0);
      stars[i].M += Moff;
      stars[i].Map += Moff;
      stars[i].MapRaw += Moff;
      stars[i].Mcalib += Moff;
    }
  }

  gfits_table_set_CMF_PS1_V5 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V5");
}

void writeStars_PS1_SV3 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_SV3 *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_SV3, Nstars);
  for (i = 0; i < Nstars; i++) {

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].detID = i;
    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;

    stars[i].posangle = 10.0;
    stars[i].pltscale = 0.25;

    stars[i].M = M[i];
    stars[i].dM = fSN;

    stars[i].Flux = flux;
    stars[i].dFlux = flux * fSN;

    stars[i].Map = M[i] - 0.05;
    stars[i].MapRaw = M[i] - 0.10;

    stars[i].apRadius = 8.0;

    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
    stars[i].apFluxErr = stars[i].apFlux * fSN;

    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
    stars[i].dMcal = 0.05;

    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
    stars[i].apNpix = 3.14*8.0*8.0;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;

    stars[i].k         = 1.0;
    stars[i].fwhmMaj   = FX*2.8;
    stars[i].fwhmMin   = FY*2.8;

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF     = 0.25;
      stars[i].psfQFperf = 0.24;
    } else {
      stars[i].psfQF     = PSFQUAL;
      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
    }

    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;

    stars[i].Mxx       = FX;
    stars[i].Mxy       = 0.01;
    stars[i].Myy       = FX;
    stars[i].M3c       = FX;
    stars[i].M3s       = FX;
    stars[i].M4c       = FX;
    stars[i].M4s       = FX;
    stars[i].Mr1       = FX;
    stars[i].Mrh       = FX;

    stars[i].kronFlux  = flux * 1.25;
    stars[i].kronFluxErr = fSN * flux * 1.25;

    stars[i].kronInner = fSN * flux * 0.9;
    stars[i].kronOuter = fSN * flux * 1.5;

    // stars[i].skyLimitRad = 1;
    // stars[i].skyLimitFlux = 2;
    // stars[i].skyLimitSlope = 0.2;

    stars[i].flags     = Flag[i];
    stars[i].flags2    = 0x80;

    stars[i].nFrames   = 1;

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      float Moff = fSN*rnd_gauss(0.0, 1.0);
      stars[i].M += Moff;
      stars[i].Map += Moff;
      stars[i].MapRaw += Moff;
      stars[i].Mcalib += Moff;
    }
  }

  gfits_table_set_CMF_PS1_SV3 (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_SV3");
}

void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {

  int i;
  CMF_PS1_V5_R2_Lensing *stars;
  float flux, fSN;

  // XXX add gaussian-distributed noise based on counts
  // this needs to make different output 'stars' entries depending on the desired type
  ALLOCATE (stars, CMF_PS1_V5_R2_Lensing, Nstars);
  for (i = 0; i < Nstars; i++) {

    flux = pow (10.0, -0.4*M[i]);
    fSN = 1.0 / sqrt(flux);

    stars[i].detID = i;
    stars[i].X = X[i];
    stars[i].Y = Y[i];
    stars[i].dX = FX * fSN;
    stars[i].dY = FY * fSN;

    stars[i].posangle = 10.0;
    stars[i].pltscale = 0.25;

    stars[i].M = M[i];
    stars[i].dM = fSN;

    stars[i].Flux = flux;
    stars[i].dFlux = flux * fSN;

    stars[i].Map = M[i] - 0.05;
    stars[i].MapRaw = M[i] - 0.10;

    stars[i].apRadius = 8.0;

    stars[i].apFlux = pow(10.0, -0.4*stars[i].Map);
    stars[i].apFluxErr = stars[i].apFlux * fSN;

    stars[i].Mcalib = M[i] + ZERO_POINT + 2.5*log10(exptime);
    stars[i].dMcal = 0.05;

    XY_to_RD (&stars[i].RA, &stars[i].DEC, X[i], Y[i], &coords);
    stars[i].apNpix = 3.14*8.0*8.0;

    stars[i].Mpeak     = M[i] + 1.0;
    stars[i].sky       = SKY;
    stars[i].dSky      = DSKY;
    stars[i].psfChisq  = PSFCHI;
    stars[i].crNsigma  = CRN;
    stars[i].extNsigma = EXTN;
    stars[i].fx        = FX;
    stars[i].fy        = FY;
    stars[i].df        = DF;

    stars[i].k         = 1.0;
    stars[i].fwhmMaj   = FX*2.8;
    stars[i].fwhmMin   = FY*2.8;

    // randomly give poor PSFQF values
    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
      stars[i].psfQF     = 0.25;
      stars[i].psfQFperf = 0.24;
    } else {
      stars[i].psfQF     = PSFQUAL;
      stars[i].psfQFperf = MAX(PSFQUAL - 0.01, 0.0);
    }

    stars[i].psfNdof   = 1;
    stars[i].psfNpix   = 2;

    stars[i].Mxx       = FX;
    stars[i].Mxy       = 0.01;
    stars[i].Myy       = FX;
    stars[i].M3c       = FX;
    stars[i].M3s       = FX;
    stars[i].M4c       = FX;
    stars[i].M4s       = FX;
    stars[i].Mr1       = FX;
    stars[i].Mrh       = FX;

    stars[i].X11_sm_obj = FX;
    stars[i].X12_sm_obj = FX;
    stars[i].X22_sm_obj = FY;
    stars[i].E1_sm_obj  = FX;
    stars[i].E2_sm_obj  = FX;

    stars[i].X11_sh_obj = FX;
    stars[i].X12_sh_obj = 0.2;
    stars[i].X22_sh_obj = FY;
    stars[i].E1_sh_obj  = FX*0.9;
    stars[i].E2_sh_obj  = FX*0.8;

    stars[i].X11_sm_psf = FX;
    stars[i].X12_sm_psf = FX;
    stars[i].X22_sm_psf = FX;
    stars[i].E1_sm_psf  = FX;
    stars[i].E2_sm_psf  = FX;

    stars[i].X11_sh_psf = FX;
    stars[i].X12_sh_psf = FX;
    stars[i].X22_sh_psf = FX;
    stars[i].E1_sh_psf  = FX;
    stars[i].E2_sh_psf  = FX;

    stars[i].kronFlux  = flux * 1.25;
    stars[i].kronFluxErr = fSN * flux * 1.25;

    stars[i].kronInner = fSN * flux * 0.9;
    stars[i].kronOuter = fSN * flux * 1.5;

    stars[i].skyLimitRad = 1;
    stars[i].skyLimitFlux = 2;
    stars[i].skyLimitSlope = 0.2;

    stars[i].flags     = Flag[i];
    stars[i].flags2    = 0x80;

    stars[i].nFrames   = 1;

    if (ADDNOISE) {
      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
      float Moff = fSN*rnd_gauss(0.0, 1.0);
      stars[i].M += Moff;
      stars[i].Map += Moff;
      stars[i].MapRaw += Moff;
      stars[i].Mcalib += Moff;
    }
  }

  gfits_table_set_CMF_PS1_V5_R2_Lensing (ftable, stars, Nstars);
  gfits_modify (ftable->header, "EXTTYPE",   "%s", 1, "PS1_V5");
}

