/*! \brief Calculate MOPS coarse ephemerides.

$Id: coarseEphem.c 1525 2007-04-13 03:02:07Z denneau $
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>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <getopt.h>

#include <mopsephem.h>


#define MAX_NUM 250

/*! \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: multiEphem [options] < orbits \n\n\
REQUIRED\n\
  --start_mjd T : starting time to compute ephemerides\n\
  --delta_mjd T : time step for each ephem\n\
  --num N : number of ephemerides to compute\n\
  --obscode OBS : observatory code for calculations, default 568 (Mauna Kea)\n\
  --limiting_mag M : upper magnitude threshold for inclusion\n\
  orbits : input orbits file\n\
OPTIONAL\n\
  --error_filename FILENAME : write errs to FILENAME\n\
  --diva_filename FILENAME : write errs to FILENAME,\n\
    or /dev/null\n\
";
    fprintf(stderr, msg);
}


static int get_opts(
        int argc, 
        char **argv, 
        char **obscode, 
        double *start_mjd,
        double *delta_mjd, 
        int *num,
        double *limiting_mag, 
        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},
            {"start_mjd", required_argument, NULL, 0},
            {"delta_mjd", required_argument, NULL, 0},
            {"num", 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, "start_mjd") == 0) {
                    *start_mjd = atof(optarg);
                }
                else if (strcmp(long_opts[option_index].name, "num") == 0) {
                    *num = atoi(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) {
        *outfile = argv[optind++];
    }

    return 0;
}


int MAIN__(int argc, char **argv) {
    double elems[6];
    double ra, dec;
    char line[1000];
    char out_buf[1000];
    int n;



    // Capture STDOUT so we can use it later under a different name.
    int stdout_fake_fid = dup(1);

    // Hide DIVA's stdout.
    int fd = open("/dev/null", O_RDWR);
    close(1);
    dup2(fd, 1);


    /* Handle options. */
    char *diva_filename = NULL;
    char *outfile = NULL;               /**< output filename */
    char *obscode = "";                 /**< MPC observatory code */
    double start_mjd = 0;               /**< start ephem time, MJD */
    double delta_mjd;                   /**< delta time for each ephem, days */
    int num;
    double limiting_mag;                /**< maximum allowed magnitude */
    char *errfile = NULL;               /**< write errors to this file */
    get_opts(argc, argv, &obscode, &start_mjd, &delta_mjd, &num, &limiting_mag, &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 || start_mjd <= 0) {
    if (start_mjd <= 0) {
        usage();
        exit(100);
    }

//    if (num > MAX_NUM) {
//        fprintf (stderr, "--num (%d) exceeds MAX_NUM (%d).\n", num, MAX_NUM);
//        exit(1);
//    }

    /* 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.
    */
    EPHEM_BLK *eph_blk;     /**< placeholder for each ephemeris */
    double *t_eph_ut;       /**< times */

    eph_blk = (EPHEM_BLK *) malloc(num * sizeof(EPHEM_BLK));
    t_eph_ut = (double *) malloc(num * sizeof(double));
    if (NULL == eph_blk || NULL == t_eph_ut) {
        fprintf(stderr, "Can't allocate eph_blk.\n");
        exit(3);
    }

    /* Copy eph_uts from command line. */
    int k;
    for (k = 0; k < num; k++) {
        t_eph_ut[k] = start_mjd + k * delta_mjd;
    }



    n = 0;
    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)) {
            // warn here
            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 our ephemerides. */
            int errcod;

            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.t_epoch_mjd = t_epoch;
            orb.g = 0;
            orb.covar = NULL;

            errcod = mopsephem_compute_ephemerides(&orb, eph_blk, t_eph_ut, num, obscode);
            if (errcod > 0) {
                fprintf(err_fp ? err_fp : stderr, "JPLEPH-ERROR: %d %s\n", errcod, objectName);   /* processing error occurred */
            }
            else {
                /* If any one of the ephemerides is visible, print the entire list. */
                int go = 0;
                int k;
                for (k = 0; k < num; k++) {
                    if (eph_blk[k].mag <= limiting_mag) {
                        go = 1;     // visible
                        break;
                    }
                }

                if (go) {
                    for (k = 0; k < num; k++) {
//                        fprintf(out_fp, "%s %.10f %.10f %.10f %.10f\n", 
                        snprintf(out_buf, 1000, "%s %.20f %.20f %.20f %.20f\n", 
                            objectName, 
                            t_eph_ut[k], 
                            eph_blk[k].ra_deg, eph_blk[k].dec_deg, eph_blk[k].mag);
                        write(stdout_fake_fid, out_buf, strlen(out_buf));
                    }
                }
            }
        }
    }

/*    fclose(out_fp);*/
    if (err_fp) 
        fclose(err_fp);
    exit(0);
}
