/** JPL orbit determination library

$Id: mopsod.c 1806 2007-08-21 17:02:12Z fpierfed $
Calculate differentially-corrected orbital parameters from IOD file and 
observations file.
*/

#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include "mopsod.h"

int mopsod_init() {
    ssd_init();
    return 0;
}


/**
 * This routine is called from Fortran to retrieve observation data populated
 * by the C interface.  We pass an opaque pointer into Fortran and use this
 * routine to retrieve the data in the Fortran side.  Note the Fortran
 * calling convention -- may be compiler specific.
 */
void mopsod_get_obs_(char *obsdata, int *i, char *obs_code, double *time_utc, 
    double *ra, double *ra_sigma_arcsec, double *dec, double *dec_sigma_arcsec, 
    double *mag_real, double *mag_sigma, 
    int *errcod, int obsdata_len, int obs_code_len) 
{
    OBSERVATION_STRUCT *obs = &((OBSERVATION_STRUCT *) obsdata)[*i];
    strncpy(obs_code, obs->obscode, obs_code_len);
    *time_utc = obs->epoch_mjd;
    *ra = obs->ra_deg;
    *dec = obs->dec_deg;
    *mag_real = obs->mag;
    *ra_sigma_arcsec = obs->ra_sigma_arcsec;
    *dec_sigma_arcsec = obs->dec_sigma_arcsec;
    *mag_sigma = obs->mag_sigma;
    *errcod = 0;
}


int mopsod_diffcor(
    ORBIT_BLK *iod,             ///< input IOD
    OBSERVATION_BLK *obs,       ///< input observations
    ORBIT_BLK *orbit            ///< output differentially corrected orbit
) {

    if (NULL == iod || NULL == obs || NULL == orbit)
        return EMOPSODBADINPUT;

    strncpy(orbit->id, iod->id, ORBIT_ID_LEN);
    orbit->ccode = 99;
    orbit->residual_arcseconds = 0; //< initialize
    orbit->cov_valid = 0;
    orbit->arc_length_days = 0;

    /* All observations have been loaded.  Run the corrector. */
    int errcod;
    int i;
    int got_min, got_max;
    double min_epoch_mjd = 0;
    double max_epoch_mjd = 0;
    double epoch;

    ssd_run_diffcor(
        (char *) obs->obs,
        &obs->nobs,

        &iod->q_au, 
        &iod->e, 
        &iod->i_deg, 
        &iod->node_deg, 
        &iod->arg_peri_deg, 
        &iod->time_peri_mjd, 
        &iod->t_epoch_mjd,

        &orbit->q_au, 
        &orbit->e, 
        &orbit->i_deg, 
        &orbit->node_deg, 
        &orbit->arg_peri_deg, 
        &orbit->time_peri_mjd, 
        &orbit->h_v, 
        &orbit->t_epoch_mjd, 
        &orbit->residual_arcseconds, 
        &orbit->covariance,
        &orbit->ccode, 
        &errcod,
        sizeof(OBSERVATION_STRUCT) * obs->nobs      ///< length of obs "string"
    );

    /* Compute the arc length all by our selves. */
    got_min = 0;        // Bools
    got_max = 0;
    for (i = 0; i < obs->nobs; i++) {
        epoch = obs->obs[i].epoch_mjd;          ///< fetch epoch for ith obs
        if (!got_min) {
            got_min = 1;
            min_epoch_mjd = epoch;
        }
        else {
            if (epoch < min_epoch_mjd) {
                min_epoch_mjd = epoch;
            }
        }

        if (!got_max) {
            got_max = 1;
            max_epoch_mjd = epoch;
        }
        else {
            if (epoch > max_epoch_mjd) {
                max_epoch_mjd = epoch;
            }
        }
    }
    orbit->arc_length_days = max_epoch_mjd - min_epoch_mjd;

    return errcod;
}
