/*! \brief Calculate MOPS coarse ephemerides.

  $Id: coarseEphem.c 1951 2007-10-03 00:01:43Z fpierfed $
  This program calculates MOPS coarse ephemerides for a supplied set of
  orbital parameters.  For each object, three positions (ephemerides)
  are calculated, separated by a specified amount of time.  The there
  positions are used by the calling script to estimate the proximity of
  a given object to a list of fields so that a reduced set of orbital
  parameters can be used to calculate accurate ephemerides.

  For each set of orbital parameters, three lines are generated, one for
  each ephemeris.  If the three ephemerides are all below the user-specified
  detection limit, no output is generated for the object.  If an error
  occurs during ephemeris calculation, no output is generated.

  The output is written to a file instead of STDOUT.  We can't use
  STDOUT because the JPL DIVA libraries spew stuff to STDOUT and
  it's difficult and annoying to filter.
*/

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>      // O_RDWR and friends

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <getopt.h>

#include <mopsephem.h>

/*! \brief Tokenize a line read in using fgets into orbital elements. */
/*!
  \param *line, the string to tokenize
  \param buffer for q, perihelion distance, in AU
  \param buffer for e, eccentricity, degrees
  \param buffer for i, inclination, degrees
  \param buffer for node, longitude of ascending node, degrees
  \param buffer for arg_peri, argument of perihelion, degrees
  \param buffer for time_peri, time of perihelion, MJD
  \param buffer for h_v, magnitude
  \param buffer for t_epoch, MJD
  \param buffer for objectName, this object's unique identifier

  Tokenize the input string, and store the tokenized values.
*/
int toke(char *line, double *q, double *e, double *i, double *node, double *arg_peri,
         double *time_peri, double *h_v, double *t_epoch, char *objectName) {

  char *tok;
  char *delim = " \t\n";
  tok = strtok(line, delim);
  if (NULL == tok)
    return 1;
  else
    *q = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *e = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *i = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *node = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *arg_peri = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *time_peri = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *h_v = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    *t_epoch = atof(tok);

  tok = strtok(NULL, delim);
  if (NULL == tok)
    return 1;
  else
    strcpy(objectName, tok);

  return 0;   // OK
}

int check_inputs(
                 double q, 
                 double e, 
                 double i, 
                 double node, 
                 double arg_peri, 
                 double time_peri, 
                 double h_v, 
                 double t_epoch) {

  /** Perform zero-order checking of inputs. */
  if (q < 0) 
    return 1;
  if (e < 0)      // non-sensical eccentricity 
    return 2;
  if (i > 180 || i < 0)
    return 3;
  //    if (time_peri < 0)
  //        return 4;

  return 0;
}


void usage(void) {
  /*! \brief Display usage message */
  char *msg = "\
usage: coarseEphem [options] ephem_time outfile < orbits\n\
\n\
REQUIRED\n\
  --obscode OBS : observatory code for calculations\n\
  --delta_mjd T : +/- delta days to use for ephemeris calculation\n\
  --limiting_mag M : upper magnitude threshold for inclusion\n\
  ephem_time : use ephem_time (in MJD) as ephemeris epoch\n\
  outfile : write output ephemerides to outfile\n\
  orbits : input orbits file\n\
OPTIONAL\n\
  --error_filename FILENAME : write errs to FILENAME\n\
  --diva_filename FILENAME : redirect JPL DIVA output to FILENAME,\n\
    or /dev/null\n\
";
  fprintf(stderr, msg);
}


static int get_opts(
                    int argc, 
                    char **argv, 
                    char **obscode, 
                    double *delta_mjd, 
                    double *limiting_mag, 
                    double *ephem_time, 
                    char **outfile,
                    char **errfile,
                    char **diva_file_ptr) 
{

  int c;
  int option_index = 0;
  while (1) {
    static struct option long_opts[] = {
      {"obscode", required_argument, NULL, 0},
      {"delta_mjd", required_argument, NULL, 0},
      {"limiting_mag", required_argument, NULL, 0},
      {"error_filename", required_argument, NULL, 0},
      {"diva_filename", required_argument, NULL, 0},
      {0, 0, 0, 0},
    };

    c = getopt_long(argc, argv, "", long_opts, &option_index);
    if (c == -1)
      break;

    switch(c) {
    case 0:
      if (strcmp(long_opts[option_index].name, "obscode") == 0) {
        *obscode = optarg;
      }
      else if (strcmp(long_opts[option_index].name, "delta_mjd") == 0) {
        *delta_mjd = atof(optarg);
      }
      else if (strcmp(long_opts[option_index].name, "limiting_mag") == 0) {
        *limiting_mag = atof(optarg);
      }
      else if (strcmp(long_opts[option_index].name, "error_filename") == 0) {
        *errfile = optarg;
      }
      else if (strcmp(long_opts[option_index].name, "diva_filename") == 0) {
        *diva_file_ptr = optarg;
      }
      break;

    default:
      usage();
      exit(100);
    }
  }

  if (optind < argc) {
    *ephem_time = atof(argv[optind++]);
  }
  if (optind < argc) {
    *outfile = argv[optind++];
  }

  return 0;
}


/* 
  The Lahey compiler requires the main to be called MAIN__ while the Intel 
  compiler wants it to be called main.
*/
#ifdef __INTEL_COMPILER
  int main(int argc, char **argv) {
#else
  int MAIN__(int argc, char **argv) {
#endif

  double elems[6];
  double ra, dec;
  char line[1000];
  int n;



  /* Handle options. */
  char *diva_filename = NULL;           
  char *outfile = NULL;                 /**< output filename */
  double ephem_time = -1;               /**< ephemeris time */
  char *obscode = "568";                /**< MPC observatory code */
  double delta_mjd = 1.0;               /**< delta time around ephem_time */
  double limiting_mag = 24.5;           /**< maximum allowed magnitude */
  char *errfile = NULL;                 /**< write errors to this file */
  get_opts(argc, argv, &obscode, &delta_mjd, &limiting_mag, &ephem_time, &outfile, &errfile, &diva_filename);


  /* Play some dirty tricks to hide DIVA noise. */
  int null_stdout_fd;
  if (diva_filename) {
    null_stdout_fd = open(diva_filename, 
        O_WRONLY | O_CREAT,
        S_IROTH | S_IWOTH | S_IRGRP | S_IWGRP | S_IRUSR | S_IWUSR
    );
  }
  else {
    null_stdout_fd = open("/dev/null", O_WRONLY);       // create new write FD for /dev/null
  }
  if (-1 == null_stdout_fd) {
    fprintf(stderr, "Could not open DIVA errors file: %s\n", diva_filename);
    exit(100);
  }
  close(1);                                             // close current STDOUT
  dup2(null_stdout_fd, 1);                              // dup /dev/null to STDOUT FD, confusing Fortran (ha ha!)

  /* Get lines describing objects from files. */
  double q, e, i, node, arg_peri, time_peri;
  double h_v, mag;
  double t_epoch;
  char objectName[100];
  double tol = 0.0001;
  FILE *out_fp = NULL;
  FILE *err_fp = NULL;

  if (!outfile || ephem_time < 0) {
    usage();
    exit(100);
  }

  /* Initialize JPL stuff integrator. */
  mopsephem_init();         /* initialize JPL integrator */


  if (errfile) {
    err_fp = fopen(errfile, "w");
    if (!err_fp) {
      fprintf(stderr, "could not open error file %s\n", errfile);
      exit(100);
    }
  }

  out_fp = fopen(outfile, "w");
  if (!out_fp) {
    fprintf(stderr, "could not open output file %s\n", outfile);
    exit(100);
  }

  /* PLAN:  For each orbit specified in the input, we will calculate three ephemerides:
     at ephem_time - delta_mjd, ephem_time + 0, and ephem_time + delta_mjd.  If the magnitude
     at any of these positions is less than our threshold, we will print out all three
     positions.
  */

  while (fgets(line, 999, stdin)) {
    if (line[0] == '#') {
      continue;   // comment line, ignore
    }
    if (toke(line, &q, &e, &i, &node, &arg_peri, &time_peri, &h_v, &t_epoch, objectName)) {
      fprintf(err_fp ? err_fp : stderr, "*** *** Line skipped: %s", line);
      continue;   // problem getting args, bail
    }

    /* Do some MOPS filtering of data.  Generally we want do not want to calculate ephemerides
       for anything close to the sun or otherwise unobservable.  Initially we have noticed that
       comets with small q and epoch near time_peri cause problems with JPL code.  So we're
       going to skip objects with q < 0.05 and fabs(time_peri - epoch) < 5. */
    int valid = check_inputs(q, e, i, node, arg_peri, time_peri, h_v, t_epoch);
    if (valid) {
      fprintf(err_fp ? err_fp : stderr, "EINVALIDINPUT %s %d\n", objectName, valid);
    }
    else if (q < 0.05 && fabs(time_peri - t_epoch) < 5) {
      fprintf(err_fp ? err_fp : stderr, "ESUNGRAZER JPLEPH-ERROR: SUN %s\n", objectName);
    }
    else {
      /* Okay.  Compute the three ephemerides. */
      EPHEM_BLK eph[3];    /**< placeholder for each ephemeris */
      int errcod;
      double ephem_times[3] = {
        ephem_time - delta_mjd,
        ephem_time,
        ephem_time + delta_mjd
      };

      ORBIT_BLK orb;
      orb.q_au = q;
      orb.e = e;
      orb.i_deg = i;
      orb.node_deg = node;
      orb.arg_peri_deg = arg_peri;
      orb.time_peri_mjd = time_peri;
      orb.h_v = h_v;
      orb.g = 0;
      orb.t_epoch_mjd = t_epoch;
      orb.covar = NULL;

      errcod = mopsephem_compute_ephemerides(&orb, eph, ephem_times, 3, obscode);

      if (errcod) {
        fprintf(err_fp ? err_fp : stderr, "JPLEPH-ERROR: %d %s\n", 
                errcod,
                objectName);   /* processing error occurred */
      }
      else if (eph[0].mag <= limiting_mag || eph[1].mag <= limiting_mag || eph[2].mag <= limiting_mag) {
        /* Object statisfies all criteria.  Print it. */
        fprintf(out_fp, "%s %.12f %.12f %.12f %.12f\n", objectName, ephem_times[0], eph[0].ra_deg, eph[0].dec_deg, eph[0].mag);
        fprintf(out_fp, "%s %.12f %.12f %.12f %.12f\n", objectName, ephem_times[1], eph[1].ra_deg, eph[1].dec_deg, eph[1].mag);
        fprintf(out_fp, "%s %.12f %.12f %.12f %.12f\n", objectName, ephem_times[2], eph[2].ra_deg, eph[2].dec_deg, eph[2].mag);
      }
    }
  }
    
  fclose(out_fp);
  if (err_fp) 
    fclose(err_fp);
  exit(0);
}
