# include "mkcmf.h"
void gauss_init (int Nbin);
double rnd_gauss (double mean, double sigma);

# 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

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

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

  int i, N, Nstars, NSTARS;
  double *X, *Y, *M;
  double Xmax, Ymax;
  float flux, fSN;

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

  PS1_DEV_1 *stars;

  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);
  }

  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);
  }

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

  f = fopen (argv[1], "r");
  if (f == NULL) {
    fprintf (stderr, "unable to open input file %s\n", argv[1]);
    exit (1);
  }
    
  Nstars = 0;
  NSTARS = 100;
  ALLOCATE (X, double, NSTARS);
  ALLOCATE (Y, double, NSTARS);
  ALLOCATE (M, double, NSTARS);

  Xmax = Ymax = 0;

  while (fscanf (f, "%lf %lf %lf", &X[Nstars], &Y[Nstars], &M[Nstars]) != EOF) {
    if (Nstars == NSTARS - 1) {
      NSTARS += 100;
      REALLOCATE (X, double, NSTARS);
      REALLOCATE (Y, double, NSTARS);
      REALLOCATE (M, double, NSTARS);
    }
    Nstars ++;
  }

  // XXX add gaussian-distributed noise based on counts
  ALLOCATE (stars, PS1_DEV_1, Nstars);
  gauss_init (2048);
  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] + FX * fSN * rnd_gauss(0.0, 1.0);
    stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0);
    stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0);

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

    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].psfQual   = PSFQUAL;
    stars[i].nFrames   = 1;
    stars[i].flags     = FLAGS;
  }

  // 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 + 50));
  gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax + 50));

  gfits_modify (&header, "NSTARS",   "%d", 1, Nstars);
  gfits_modify (&header, "PHOTCODE", "%s", 1, "SIMTEST.r.Chip");
  gfits_modify (&header, "DATE-OBS", "%s", 1, date);
  gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
  gfits_modify (&header, "ZERO_PT", "%lf", 1, 25.0);
  gfits_modify (&header, "EXPTIME", "%lf", 1, 1.0);
  gfits_modify (&header, "AIRMASS", "%lf", 1, 1.0);
  gfits_modify (&header, "NASTRO",   "%d", 1, 10); 

  gfits_modify (&header, "IMAGEID",  "%d", 1, (int)(1000*drand48())); 
  gfits_modify (&header, "SOURCEID", "%d", 1, (int)(100*drand48())); 

  /* bore site center guess */
  strcpy (coords.ctype, "DEC--TAN");
  coords.crval1 = RA;
  coords.crval2 = DEC;
  coords.crpix1 = 0.0;
  coords.crpix2 = 0.0;
  
  coords.cdelt1 = 0.25/3600.0;
  coords.cdelt2 = 0.25/3600.0;

  coords.pc1_1  = 1;
  coords.pc2_2  = 1;
  coords.pc1_2  = 0;
  coords.pc2_1  = 0;
  coords.Npolyterms = 1;

  PutCoords (&coords, &header);
  gfits_modify (&header, "EXTNAME",   "%s", 1, "Chip.hdr");

  ftable.header = &theader;
  gfits_table_set_PS1_DEV_1 (&ftable, stars, Nstars);
  gfits_modify (&theader, "EXTNAME",   "%s", 1, "Chip.psf");
  gfits_modify (&theader, "EXTTYPE",   "%s", 1, "PS1_DEV_1");
  gfits_modify (&theader, "EXTHEAD",   "%s", 1, "Chip.hdr");

  fits = fopen (argv[2], "w");
  if (fits == NULL) {
    fprintf (stderr, "ERROR: can't open output file %s\n", argv[2]);
    exit (1);
  }

  gfits_fwrite_header  (fits, &header);
  gfits_fwrite_matrix  (fits, &matrix);
  gfits_fwrite_Theader (fits, &theader);
  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;
  long A, B;
  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);
  for (B = 0; A == time(NULL); B++);
  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);
}
 
