# include "astro.h"
# define TIME {gettimeofday (&now, (struct timeval *) NULL); \
fprintf (stderr, "start time: %f\n", now.tv_sec+now.tv_usec*1e-6);}
double rnd_gauss();
# define T970101 852076800
# define DRIFT_RATE 33.168 /* msec per row */ 
# define MAG_LIMIT 20.100

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

  int i, N1, N2, NX, NY, Ncat, Nregions, nstar;
  char **regions, **gregions();
  Stars *catalog, *stars1, *stars2, *stars3;
  struct timeval now;  
  Coords coords;
  Header header;
  FILE *f;
  double X, Y, M, Rsigma, Msigma, Moffset;

  if (argc < 7) {
    fprintf (stderr, "USAGE: simulate scale ra dec Moffset Msigma Rsigma template output\n");
    exit (0);
  }
  
  Moffset = atof (argv[4]);
  Msigma = atof (argv[5]);
  Rsigma = atof (argv[6]);

  /* load in template file, determine size of "image" */
  f = fopen (argv[7], "r");
  if (f == NULL) {
    fprintf (stderr, "could not find template file %s\n", argv[7]);
    exit(0);
  }
  ALLOCATE (header.buffer, char, 8640);
  fread (header.buffer, 1, 8640, f);
  header.size = 8640;
  fits_scan (&header, "NAXIS1", "%d", 1, &NX);
  fits_scan (&header, "NAXIS2", "%d", 1, &NY);
  gettimeofday (&now, (struct timeval *) NULL);
  fprintf (stderr, "now: %d.%06d\n", now.tv_sec, now.tv_usec);
  fits_modify (&header, "TZERO", "%d", 1, now.tv_sec - T970101);
  fits_modify (&header, "TRATE", "%lf", 1, DRIFT_RATE);
  fits_modify (&header, "FLIMIT", "%lf", 1, MAG_LIMIT);
  fclose (f);

  /* set up coordinates for image */
  strcpy (coords.ctype, "RA---TAN");
  coords.crpix1 = 0.5*NX;
  coords.crpix2 = 0.5*NY;
  coords.pc1_1 = 1; coords.pc1_2 = 0;
  coords.pc2_1 = 0; coords.pc2_2 = 1;
  coords.cdelt1 = coords.cdelt2 = atof (argv[1]) / 3600.0;
  coords.crval1 = atof (argv[2]);
  coords.crval2 = atof (argv[3]);
  fits_modify (&header, "RA", "%lf", 1, coords.crval1);
  fits_modify (&header, "DEC", "%lf", 1, coords.crval2);

  /* find appropriate catalogs, load files */
  regions = gregions (coords, &Nregions);
  for (N2 = i = 0; i < Nregions; i++) {
    gcatalog (regions[i], &catalog, &Ncat);
  }
  gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY);
  free (catalog);

  /* select catalog star within image region */ 
  ALLOCATE (stars3, Stars, N2);
  for (nstar = i = 0; i < N2; i++) {
    if ((stars2[i].X > 0) && (stars2[i].X < NX) && 
	(stars2[i].Y > 0) && (stars2[i].Y < NY)) {
      stars3[nstar] = stars2[i];
      nstar++;
    }
  }
  free (stars2);
  fits_modify (&header, "NSTARS", "%d", 1, nstar);

  /* write header */
  f = fopen (argv[8], "w");
  if (f == NULL) {
    fprintf (stderr, "could not open output file %s\n", argv[8]);
    exit(0);
  }
  fwrite (header.buffer, 1, 8640, f);
  
  /* write data */
  gauss_init ();
  for (i = 0; i < nstar; i++) {
    X = rnd_gauss (stars3[i].X, Rsigma);
    Y = rnd_gauss (stars3[i].Y, Rsigma);
    M = rnd_gauss (stars3[i].mag, Msigma) + Moffset;
    fprintf (f, "%6.1f %6.1f %6.3f %03.0f 1 3.2\n", 
	     X, Y, M, Msigma*1000.0);
  }
  fclose (f);

}

