#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "qfits.h"
#include <getopt.h>


struct image
{
  int nx, ny;
  char timestamp[64];
  unsigned char *c;
};
struct fourframe
{
  int nx, ny, nt, nlayer;
  char timestamp[64], observer[64];
  double mjd;
  float *z, *dt;
  int cospar, type;
};
int fgetline (FILE * file, char *s, int lim);
void write_fits (char *filename, struct fourframe ff);
struct image read_pgm (char *filename, int *status);
double nfd2mjd (char *date);
double date2mjd (int year, int month, double day);
void write_pgm (char *filename, struct fourframe ff);
void mjd2date (double mjd, char *date);

void
usage (void)
{
  printf ("pgm2fits p:w:h:s:n:Dd:x:y:c:o:gm:t:r:I\n\n");
  printf ("-p   image prefix\n");
  printf ("-w   image width in pixels\n");
  printf ("-h   image height in pixels\n");
  printf ("-s   number of first image to process\n");
  printf ("-n   number of images to process\n");
  printf ("-D   toggle for creating dark frame\n");
  printf ("-d   toggle for creating dark frame\n");
  printf ("-d   filename of dark frame to substract\n");
  printf ("-m   filename of mask frame to apply\n");
  printf ("-x   tracking rate in x (pix/s)\n");
  printf ("-y   tracking rate in y (pix/s)\n");
  printf ("-c   COSPAR [default 4553]\n");
  printf ("-o   observer [default \"Cees Bassa\"]\n");
  printf ("-g   toggle for guiding??\n");
  printf ("-t   time stamp of first image [YYYY-MM-DDTHH:MM:SS.SSS]\n");
  printf ("-r   frame rate (frames/s)\n");
  printf ("-I   integer output\n");
  exit (0);
}

int
main (int argc, char *argv[])
{
  int i, j, k, l, m, k0 = 1, status, nt, darkout = 0, darkin = 0, track =
    0, maskin = 0;
  int n, n0, di, dj, npix;
  struct image *img, drk, msk;
  struct fourframe ff;
  char filename[128], nfd[64];
  float s1, s2, z;
  float avg, std, max, cnt, *trk;
  int *wt;
  int arg = 0;
  char *path, *darkfile, *maskfile;
  double mjd, mjd0 = 0.0;
  float dxdn = 0.0, dydn = 0.0, dx, dy;
  int guide = 0, timereset = 0;
  float framerate = 25.0;
  char *env;

  // Set defaults
  env = getenv ("ST_COSPAR");
  ff.cospar = atoi (env);
  strcpy (ff.observer, "Cees Bassa");
  ff.type = 0;

  // Decode options
  if (argc > 1)
    {
      while ((arg =
	      getopt (argc, argv, "p:w:h:s:n:Dd:x:y:c:o:gm:t:r:I")) != -1)
	{
	  switch (arg)
	    {
	    case 'p':
	      path = optarg;
	      break;

	    case 'g':
	      guide = 1;
	      break;

	    case 'w':
	      ff.nx = atoi (optarg);
	      break;

	    case 'h':
	      ff.ny = atoi (optarg);
	      break;

	    case 'I':
	      ff.type = 1;
	      break;

	    case 'c':
	      ff.cospar = atoi (optarg);
	      break;

	    case 'o':
	      strcpy (ff.observer, optarg);
	      break;

	    case 's':
	      k0 = atoi (optarg);
	      break;

	    case 'n':
	      nt = atoi (optarg);
	      break;

	    case 'D':
	      darkout = 1;
	      break;

	    case 'd':
	      darkin = 1;
	      darkfile = optarg;
	      break;

	    case 'm':
	      maskin = 1;
	      maskfile = optarg;
	      break;

	    case 'x':
	      dxdn = atof (optarg);
	      track = 1;
	      break;

	    case 'y':
	      dydn = atof (optarg);
	      track = 1;
	      break;

	    case 't':
	      strcpy (nfd, optarg);
	      mjd0 = nfd2mjd (nfd);
	      timereset = 1;
	      break;

	    case 'r':
	      framerate = atof (optarg);
	      timereset = 1;
	      break;

	    default:
	      usage ();
	    }
	}
    }
  else
    {
      usage ();
    }

  // Add layer
  if (track == 1)
    ff.nlayer = 5;
  else
    ff.nlayer = 4;

  // Allocate
  ff.z = (float *) malloc (ff.nlayer * sizeof (float) * ff.nx * ff.ny);
  ff.dt = (float *) malloc (sizeof (float) * nt);
  trk = (float *) malloc (sizeof (float) * ff.nx * ff.ny);
  wt = (int *) malloc (sizeof (float) * ff.nx * ff.ny);
  img = (struct image *) malloc (sizeof (struct image) * nt);

  // Read dark file
  if (darkin == 1)
    drk = read_pgm (darkfile, &status);

  // Read mask file
  if (maskin == 1)
    msk = read_pgm (maskfile, &status);

  // Loop over files
  for (k = 0, l = 0; k < nt; k++)
    {
      sprintf (filename, "%s%06d.pgm", path, k + k0);
      img[l] = read_pgm (filename, &status);

      // Reset time
      if (timereset == 1)
	{
	  mjd = mjd0 + (double) (k + k0) / (86400.0 * framerate);
	  mjd2date (mjd, img[l].timestamp);
	}

      ff.dt[l] =
	86400.0 * (nfd2mjd (img[l].timestamp) - nfd2mjd (img[0].timestamp));
      if (status == 0)
	{
	  printf ("Read %s\n", filename);
	  l++;
	}
      else
	{
	  break;
	}
    }
  ff.nt = l;
  strcpy (ff.timestamp, img[0].timestamp);
  ff.mjd = nfd2mjd (img[0].timestamp);

  printf ("Accumulating image statistics\n");
  // Loop over pixels
  for (i = 0; i < ff.nx; i++)
    {
      for (j = 0; j < ff.ny; j++)
	{
	  n = i + ff.nx * j;

	  s1 = 0.0;
	  s2 = 0.0;
	  max = 0.0;
	  cnt = 0.0;

	  // Loop over images
	  for (k = 0; k < ff.nt; k++)
	    {
	      if (darkin == 0)
		z = (float) img[k].c[n];
	      else if (darkin == 1)
		z = (float) (img[k].c[n] - drk.c[n]);
	      s1 += z;
	      s2 += z * z;
	      if (z > max)
		{
		  max = z;
		  cnt = (float) k;
		}
	    }
	  s1 -= max;
	  s2 -= max * max;
	  avg = s1 / (float) (ff.nt - 1);
	  std = sqrt ((s2 - s1 * avg) / (float) (ff.nt - 2));

	  // Reset masked pixels
	  if (maskin == 1 && msk.c[n] == 0.0)
	    {
	      avg = 0.0;
	      std = 0.0;
	      max = 0.0;
	      cnt = 128.0;
	    }

	  for (m = 0; m < ff.nlayer; m++)
	    {
	      l = i + (ff.ny - j - 1) * ff.nx + m * ff.nx * ff.ny;
	      if (m == 0)
		ff.z[l] = avg;
	      if (m == 1)
		ff.z[l] = std;
	      if (m == 2)
		ff.z[l] = max;
	      if (m == 3)
		ff.z[l] = cnt;
	      if (m == 4)
		ff.z[l] = avg;
	    }
	}
    }

  // Create tracked layer
  if (track == 1)
    {
      printf ("Creating tracked layer\n");

      // Set weights
      for (i = 0; i < ff.nx * ff.ny; i++)
	{
	  wt[i] = 0;
	  trk[i] = 0.0;
	}

      // Loop over frames
      for (l = 0; l < ff.nt; l++)
	{
	  // Offset
	  dx = dxdn * (l - ff.nt / 2);
	  dy = dydn * (l - ff.nt / 2);

	  // Integer offset
	  di = (int) floor (dx + 0.5);
	  dj = (int) floor (dy + 0.5);

	  // Set
	  for (i = 0; i < ff.nx; i++)
	    {
	      for (j = 0; j < ff.ny; j++)
		{
		  k = i + ff.nx * j;
		  k0 = i + di + ff.nx * (j + dj);
		  if (i + di > 0 && i + di < ff.nx && j + dj > 0
		      && j + dj < ff.ny)
		    {
		      wt[k] += 1;
		      trk[k] += (float) img[l].c[k0];
		    }
		}
	    }
	}
      // Save layer
      for (i = 0; i < ff.nx; i++)
	{
	  for (j = 0; j < ff.ny; j++)
	    {
	      k = i + ff.nx * j;
	      if (guide == 0)
		l = i + (ff.ny - j - 1) * ff.nx;
	      else
		l = i + (ff.ny - j - 1) * ff.nx + 4 * ff.nx * ff.ny;
	      if (wt[k] > 0)
		ff.z[l] = trk[k] / (float) wt[k];
	      else
		ff.z[l] = trk[k];
	    }
	}
    }

  // Write fits
  if (ff.timestamp != NULL)
    sprintf (filename, "%s.fits", ff.timestamp);
  else
    strcpy (filename, "test.fits");
  write_fits (filename, ff);

  // Write dark frame
  if (darkout == 1)
    write_pgm ("dark.pgm", ff);

  // Free
  for (k = 0; k < ff.nt; k++)
    free (img[k].c);
  free (trk);
  free (img);
  free (wt);
  free (ff.dt);
  free (ff.z);

  return 0;
}

// Read pgm file
struct image
read_pgm (char *filename, int *status)
{
  int i;
  struct image img;
  FILE *file;
  char hbuf[64];

  // Open file
  file = fopen (filename, "rb");
  if (file == NULL)
    {
      *status = 1;
      return img;
    }

  // Read PGM format
  fgetline (file, hbuf, 64);
  if (strcmp (hbuf, "P5") != 0)
    {
      printf ("Not a valid PGM file!\n");
      exit (0);
    }

  // Read timestamp/image size
  fgetline (file, hbuf, 64);
  if (strstr (hbuf, "#") != NULL)
    {
      strcpy (img.timestamp, hbuf + 2);
      fgetline (file, hbuf, 64);
      sscanf (hbuf, "%d %d", &img.nx, &img.ny);
    }
  else
    {
      strcpy (img.timestamp, "2012-01-01T00:00:00");
      sscanf (hbuf, "%d %d", &img.nx, &img.ny);
    }
  fgetline (file, hbuf, 64);

  // Allocate
  img.c = (unsigned char *) malloc (sizeof (unsigned char) * img.nx * img.ny);

  // Read
  fread (img.c, 1, img.nx * img.ny, file);

  // Close file
  fclose (file);

  *status = 0;
  return img;
}

// Get line
int
fgetline (FILE * file, char *s, int lim)
{
  int c, i = 0;

  while (--lim > 0 && (c = fgetc (file)) != EOF && c != '\n')
    s[i++] = c;
  //  if (c=='\n')
  //    s[i++]=c;
  s[i] = '\0';

  return i;
}

// Write fits file
void
write_fits (char *filename, struct fourframe ff)
{
  int i, j, k, l;
  float *fbuf;
  int *ibuf;
  qfitsdumper qd;
  qfits_header *qh;
  char key[FITS_LINESZ + 1];
  char val[FITS_LINESZ + 1];
  char com[FITS_LINESZ + 1];
  char lin[FITS_LINESZ + 1];
  FILE *file;

  // Create FITS header
  qh = qfits_header_default ();

  // Add stuff
  if (ff.type == 1)
    qfits_header_add (qh, "BITPIX", "8", " ", NULL);
  else
    qfits_header_add (qh, "BITPIX", "-32", " ", NULL);
  qfits_header_add (qh, "NAXIS", "3", " ", NULL);
  sprintf (val, "%i", ff.nx);
  qfits_header_add (qh, "NAXIS1", val, " ", NULL);
  sprintf (val, "%i", ff.ny);
  qfits_header_add (qh, "NAXIS2", val, " ", NULL);
  sprintf (val, "%i", ff.nlayer);
  qfits_header_add (qh, "NAXIS3", val, " ", NULL);
  qfits_header_add (qh, "BSCALE", "1.0", " ", NULL);
  qfits_header_add (qh, "BZERO", "0.0", " ", NULL);
  qfits_header_add (qh, "DATAMAX", "255.0", " ", NULL);
  qfits_header_add (qh, "DATAMIN", "0.0", " ", NULL);
  sprintf (val, "%s", ff.timestamp);
  qfits_header_add (qh, "DATE-OBS", val, " ", NULL);

  // MJD-OBS
  sprintf (val, "%lf", ff.mjd);
  qfits_header_add (qh, "MJD-OBS", val, " ", NULL);
  sprintf (val, "%f", ff.dt[ff.nt - 1] - ff.dt[0]);
  qfits_header_add (qh, "EXPTIME", val, " ", NULL);
  sprintf (val, "%d", ff.nt);
  qfits_header_add (qh, "NFRAMES", val, " ", NULL);

  // Astrometry keywors
  sprintf (val, "%f", ff.nx / 2.0);
  qfits_header_add (qh, "CRPIX1", val, " ", NULL);
  sprintf (val, "%f", ff.ny / 2.0);
  qfits_header_add (qh, "CRPIX2", val, " ", NULL);
  qfits_header_add (qh, "CRVAL1", "0.0", " ", NULL);
  qfits_header_add (qh, "CRVAL2", "0.0", " ", NULL);
  qfits_header_add (qh, "CD1_1", "0.0", " ", NULL);
  qfits_header_add (qh, "CD1_2", "0.0", " ", NULL);
  qfits_header_add (qh, "CD2_1", "0.0", " ", NULL);
  qfits_header_add (qh, "CD2_2", "0.0", " ", NULL);
  qfits_header_add (qh, "CTYPE1", "'RA---TAN'", " ", NULL);
  qfits_header_add (qh, "CTYPE2", "'DEC--TAN'", " ", NULL);
  qfits_header_add (qh, "CUNIT1", "'deg'", " ", NULL);
  qfits_header_add (qh, "CUNIT2", "'deg'", " ", NULL);
  qfits_header_add (qh, "CRRES1", "0.0", " ", NULL);
  qfits_header_add (qh, "CRRES2", "0.0", " ", NULL);
  qfits_header_add (qh, "EQUINOX", "2000.0", " ", NULL);
  qfits_header_add (qh, "RADECSYS", "ICRS", " ", NULL);
  sprintf (val, "%d", ff.cospar);
  qfits_header_add (qh, "COSPAR", val, " ", NULL);
  sprintf (val, "'%s'", ff.observer);
  qfits_header_add (qh, "OBSERVER", val, " ", NULL);

  // Add timestamps
  for (k = 0; k < ff.nt; k++)
    {
      sprintf (key, "DT%04d", k);
      sprintf (val, "%f", ff.dt[k]);
      qfits_header_add (qh, key, val, " ", NULL);
    }

  // Dummy keywords
  for (k = 0; k < 10; k++)
    {
      sprintf (key, "DUMY%03d", k);
      qfits_header_add (qh, key, "0.0", " ", NULL);
    }

  // Dump fitsheader
  //  qfits_header_dump(qh,stdout);

  // Dump to file
  file = fopen (filename, "w");
  qfits_header_dump (qh, file);
  fclose (file);


  // Fill buffer
  fbuf = malloc (ff.nlayer * ff.nx * ff.ny * sizeof (float));
  ibuf = malloc (ff.nlayer * ff.nx * ff.ny * sizeof (int));
  for (i = 0, l = 0; i < ff.nx; i++)
    {
      for (j = ff.ny - 1; j >= 0; j--)
	{
	  for (k = 0; k < ff.nlayer; k++)
	    {
	      fbuf[l] = ff.z[l];
	      ibuf[l] = (int) ff.z[l];

	      l++;
	    }
	}
    }

  // Set parameters
  qd.filename = filename;
  qd.npix = ff.nlayer * ff.nx * ff.ny;
  if (ff.type == 1)
    {
      qd.ptype = PTYPE_INT;
      qd.ibuf = ibuf;
      qd.out_ptype = BPP_8_UNSIGNED;
    }
  else
    {
      qd.ptype = PTYPE_FLOAT;
      qd.fbuf = fbuf;
      qd.out_ptype = -32;
    }

  // Dump
  qfits_pixdump (&qd);

  free (fbuf);
  free (ibuf);

  return;
}

// Compute Julian Day from Date
double
date2mjd (int year, int month, double day)
{
  int a, b;
  double jd;

  if (month < 3)
    {
      year--;
      month += 12;
    }

  a = floor (year / 100.);
  b = 2. - a + floor (a / 4.);

  if (year < 1582)
    b = 0;
  if (year == 1582 && month < 10)
    b = 0;
  if (year == 1582 && month == 10 && day <= 4)
    b = 0;

  jd =
    floor (365.25 * (year + 4716)) + floor (30.6001 * (month + 1)) + day + b -
    1524.5;

  return jd - 2400000.5;
}

// nfd2mjd
double
nfd2mjd (char *date)
{
  int year, month, day, hour, min;
  double mjd, dday;
  float sec;

  sscanf (date, "%04d-%02d-%02dT%02d:%02d:%f", &year, &month, &day, &hour,
	  &min, &sec);
  dday = day + hour / 24.0 + min / 1440.0 + sec / 86400.0;
  mjd = date2mjd (year, month, dday);

  return mjd;
}

// Write pgm file
void
write_pgm (char *filename, struct fourframe ff)
{
  int i, j, k;
  FILE *file;

  file = fopen (filename, "w");
  fprintf (file, "P5\n# 2013-01-01T00:00:00\n%d %d\n255\n", ff.nx, ff.ny);
  for (j = 0; j < ff.ny; j++)
    {
      for (i = 0; i < ff.nx; i++)
	{
	  k = i + (ff.ny - j - 1) * ff.nx;
	  fprintf (file, "%c", (char) ff.z[k]);
	}
    }
  fclose (file);

  return;
}

// Compute Date from Julian Day
void
mjd2date (double mjd, char *date)
{
  double f, jd, dday;
  int z, alpha, a, b, c, d, e;
  int year, month, day, hour, min;
  double sec, x;

  jd = mjd + 2400000.5;
  jd += 0.5;

  z = floor (jd);
  f = fmod (jd, 1.);

  if (z < 2299161)
    a = z;
  else
    {
      alpha = floor ((z - 1867216.25) / 36524.25);
      a = z + 1 + alpha - floor (alpha / 4.);
    }
  b = a + 1524;
  c = floor ((b - 122.1) / 365.25);
  d = floor (365.25 * c);
  e = floor ((b - d) / 30.6001);

  dday = b - d - floor (30.6001 * e) + f;
  if (e < 14)
    month = e - 1;
  else
    month = e - 13;

  if (month > 2)
    year = c - 4716;
  else
    year = c - 4715;

  day = (int) floor (dday);
  x = 24.0 * (dday - day);
  x = 3600. * fabs (x);
  sec = fmod (x, 60.);
  x = (x - sec) / 60.;
  min = fmod (x, 60.);
  x = (x - min) / 60.;
  hour = x;

  sprintf (date, "%04d-%02d-%02dT%02d:%02d:%06.3f", year, month, day, hour,
	   min, sec);

  return;
}