package PS::MOPS::JPLEPH;

use 5.008008;
use strict;
use warnings;
use Carp;

require Exporter;

our @ISA = qw(Exporter);
our @EXPORT = qw(
    jpleph_calcEphemeris
    jpleph_calcDetectionResiduals
    jpleph_calcEphemerides
    jpleph_propagateElements
);

our $VERSION = '0.01';

require XSLoader;
XSLoader::load('PS::MOPS::JPLEPH', $VERSION);

# Preloaded methods go here.
use PS::MOPS::Constants qw(:all);
use Params::Validate;
use Astro::SLA;


# Global enums.
our $EBADINPUT = 100;
our $EBADQ = 101;
our $EBADE = 102;
our $EBADI = 103;
our $ESUNGRAZER = 104;
our $EJPLEPH = 200;


our $init;
sub _check_init {
    if (!$init) {
        PS::MOPS::JPLEPH::ssd_init_c();     # initialize JPL integrator, heavy startup
        $init = 1;
    }
}


sub _check_elements {
    # Check various validity or performance criteria of the input orbital elements.  Return zero or undef
    # if the orbit is valid.
    my ($q_AU, $e, $i_deg, $node_deg, $argPeri_deg, $timePeri_MJD, $hV, $epoch_MJD, $t_eph_ut_MJD) = @_;

    return $EBADQ if ($q_AU < 0);
    return $EBADE if ($e < 0);
    return $EBADI if ($i_deg > 180 or $i_deg < -90);
    return $ESUNGRAZER if ($q_AU < 0.05 and abs($timePeri_MJD - $t_eph_ut_MJD) < 5);
    return 0;
}


sub jpleph_calcEphemeris {
    # Compute an ephemeris by contacting the local ephemeris server.
    my %args = validate(@_, {
        orbit => 1,                 # orbit of object
        obscode => 1,               # observatory code
        epoch_mjd => 0,             # epoch of prediction
        epoch => 0,                 # epoch of prediction (deprecated)
        trails => 0,                # include trailed detection info
        exposureTime_sec => 0,      # exposure time for trailed dets
        objectName => 0,            # override orbitId if specified
        withSRC => 0,               # use SRC to compute positional errors?
    });
    my $orb = $args{orbit};
    my $epoch_mjd = ($args{epoch_mjd} || $args{epoch});
    my $objectName = ($args{objectName} || $orb->orbitId);
    my $res;                        # result

    # Initialize the JPL integrator if necessary.
    _check_init();

    # Check our inptus.
    my ($q_AU, $e, $i_deg, $node_deg, $argPeri_deg, $timePeri_MJD, $hV, $epoch_MJD, $t_eph_ut_MJD) = 
        (@{$orb}{qw(q e i node argPeri timePeri hV epoch)}, $epoch_mjd);
    return undef if _check_elements($q_AU, $e, $i_deg, $node_deg, $argPeri_deg, $timePeri_MJD, $hV, $epoch_MJD, $t_eph_ut_MJD);


    # If trails is specified, get dual RA/Dec pairs, and compute the trail parameters.
    if ($args{trails}) {
        die "exposureTime_sec not specified" unless exists($args{exposureTime_sec});
        my ($ra_deg, $dec_deg, $mag, $ra2_deg, $dec2_deg, $errcod) = PS::MOPS::JPLEPH::ssd_compute_ephemeris_pair_c(
            @{$orb}{qw(epoch q e i node argPeri timePeri hV)}, 
            $args{obscode}, $epoch_mjd, $args{exposureTime_sec},
        );

        my ($length_deg, $orient_deg);      # trail info
        my ($ra_rad, $dec_rad, $ra2_rad, $dec2_rad) = 
            map { $_ / $DEG_PER_RAD } 
            ($ra_deg, $dec_deg, $ra2_deg, $dec2_deg);   # convert deg to rad

        $length_deg = slaDsep($ra_rad, $dec_rad, $ra2_rad, $dec2_rad) * $DEG_PER_RAD;
        $orient_deg = slaBear($ra_rad, $dec_rad, $ra2_rad, $dec2_rad) * $DEG_PER_RAD;
        if ($errcod) {
            $res = {
                objectName => $objectName,
                errcod => $errcod,
            };
        }
        else {
            $res = {
                objectName => $objectName,
                epoch => $epoch_mjd,
                ra => $ra_deg,
                dec => $dec_deg,
                mag => $mag,
                orient_deg => $orient_deg,
                length_deg => $length_deg,
            };
        }
    }
    else {
        my ($ra_deg, $dec_deg, $mag, $smaa_sec, $smia_sec, $pa_deg, $errcod);
        if(!$orb->{covariance}) {
            $args{withSRC} = 0;
        }
        if($args{withSRC}) {
            my ($cov_01, $cov_02, $cov_03, $cov_04, $cov_05, $cov_06, $cov_07, 
                $cov_08, $cov_09, $cov_10, $cov_11, $cov_12, $cov_13, $cov_14, 
                $cov_15, $cov_16, $cov_17, $cov_18, $cov_19, $cov_20, $cov_21) = @{$orb->{covariance}};
            ($ra_deg, 
             $dec_deg, 
             $mag, 
             $smaa_sec, 
             $smia_sec, 
             $pa_deg, 
             $errcod) = PS::MOPS::JPLEPH::ssd_compute_ephemeris_with_src_c(
                @{$orb}{qw(epoch q e i node argPeri timePeri hV)}, 
                $args{obscode}, $epoch_mjd,
                $cov_01, $cov_02, $cov_03, $cov_04, $cov_05, $cov_06, $cov_07, 
                $cov_08, $cov_09, $cov_10, $cov_11, $cov_12, $cov_13, $cov_14, 
                $cov_15, $cov_16, $cov_17, $cov_18, $cov_19, $cov_20, $cov_21,
            );
        } else {
            ($ra_deg, 
             $dec_deg, 
             $mag, 
             $errcod) = PS::MOPS::JPLEPH::ssd_compute_ephemeris_c(
                @{$orb}{qw(epoch q e i node argPeri timePeri hV)}, 
                $args{obscode}, $epoch_mjd,
            );
            $smaa_sec = -1;
            $smia_sec = -1;
            $pa_deg = -1;
        }
        
        if ($errcod) {
            $res = {
                objectName => $objectName,
                errcod => $errcod,
            };
        }
        else {
            $res = {
                objectName => $objectName,
                epoch => $epoch_mjd,
                ra => $ra_deg,
                dec => $dec_deg,
                mag => $mag,
                smaa_sec => $smaa_sec,
                smia_sec => $smia_sec,
                pa_deg => $pa_deg,
            };
        }
    }

    return $res;
}


sub jpleph_calcDetectionResiduals {
    # Compute per-detection residuals given an orbit and a list of detections.  Compute the
    # predicted position of the orbit using the epochs and obscodes of the supplied detections.
    my %args = validate(@_, {
        orbit => 1,                 # orbit of object
        detections_ref => 1,        # ARRAYREF of detections to analyze
    });

    # Connect to ephem server if not connected.
    my $det;
    my %resids;
    foreach $det (@{$args{detections_ref}})  {
        my $eph_href = PS::MOPS::JPLEPH::jpleph_calcEphemeris(
            orbit => $args{orbit},
            obscode => $det->obscode,
            epoch => $det->epoch,
        );
        $resids{$det->detId} = {
            epoch => $det->epoch,
            dra => $det->ra - $eph_href->{ra},
            ddec => $det->dec - $eph_href->{dec},
            dmag => $det->mag - $eph_href->{mag},
            dtot => slaDsep(
                $det->ra / $DEG_PER_RAD, 
                $det->dec / $DEG_PER_RAD, 
                $eph_href->{ra} / $DEG_PER_RAD, 
                $eph_href->{dec} / $DEG_PER_RAD
            ) * $DEG_PER_RAD,        # cvt to degrees
        };
    }

    return \%resids;
}


sub jpleph_calcEphemerides {
    # Most-basic and computationally expensive method of calculating
    # ephemerides.  Has to fetch something an orbit definition from 
    # the database, then crunch on it.
    use PS::MOPS::DC::SSM;
    use PS::MOPS::DC::Orbit;

    # Inputs
    #   objectNames : arrayref of objectNames
    #   epoch : epoch to calculate
    # Return
    #   hashref of results, keys objectName
    #           objectName => $objectName,
    #           epoch => $epoch,
    #           ra => $ra,
    #           dec => $dec,
    #           mag => $mag,
    my %args = validate(@_, {
        instance => 1,              # MOPS instance
        objectNames => 0,           # list of SSM object names
        derivedObjectNames => 0,    # list of derived object names
        epoch_mjd => 0,             # epoch for ephemeris calc
        epoch => 0,                 # alias for epoch_mjd
        obscode => 1,               # observatory code
        trails => 0,                # include trailed detection info
        exposureTime_sec => 0,      # exposure time for trailed dets
        host => 0,                  # ignored: optional hostname for remote job
        user => 0,                  # ignored: username@host
        withSRC => 0,               # shoall we use the SRC to compute errors?
    });
    my %results;    # calculated ephemerides
    croak "neither epoch_mjd nor epoch specified" unless exists($args{epoch_mjd}) or exists($args{epoch});
    $args{epoch_mjd} = $args{epoch} unless defined($args{epoch_mjd});

    my $orb;        # orbit object
    my $orbits;
    if ($args{objectNames}) {
        $orbits = modcs_retrieve($args{instance}, objectNames => $args{objectNames});
    }
    elsif ($args{derivedObjectNames}) {
        $orbits = modco_retrieve($args{instance}, derivedobjectNames => $args{derivedObjectNames});
    }
    else {
        return %results;        # pass buck, return empty something
    }

    # note: JPLEPH DIVA lib noise redirected to /dev/null
    foreach $orb (@{$orbits}) {
        $results{$orb->objectName} = jpleph_calcEphemeris(
            orbit => $orb,
            objectName => $orb->objectName,
            obscode => $args{obscode},
            epoch_mjd => $args{epoch_mjd},
            trails => $args{trails},
            exposureTime_sec => $args{exposureTime_sec},
            withSRC => $args{withSRC},
        )
    }

    return \%results;
}


sub jpleph_propagateElements {
    # Propagate the input orbit to the specified epoch, and return an orbit-like
    # object containing the propagated elements.
    my %args = validate(@_, {
        orbit => 1,                 # orbit of object
        epoch_mjd => 1,             # epoch of prediction
    });
    my $orb = $args{orbit};
    my $epoch_mjd = ($args{epoch_mjd} || $args{epoch});
    my $res;                        # result

    # Initialize the JPL integrator if necessary.
    _check_init();

    # Check our inptus.
    my ($q_AU, $e, $i_deg, $node_deg, $argPeri_deg, $timePeri_MJD, $hV, $epoch_MJD, $t_eph_ut_MJD) = 
        (@{$orb}{qw(q e i node argPeri timePeri hV epoch)}, $epoch_mjd);
    return undef if _check_elements($q_AU, $e, $i_deg, $node_deg, $argPeri_deg, $timePeri_MJD, $hV, $epoch_MJD, $t_eph_ut_MJD);


    my ($out_q, $out_e, $out_i, $out_node, $out_argPeri, $out_timePeri, $errcod) = PS::MOPS::JPLEPH::ssd_propagate_elements_c(
        @{$orb}{qw(epoch q e i node argPeri timePeri)}, 
        $epoch_mjd,
    );
    return undef if $errcod;

    $res = {
        epoch => $epoch_mjd,
        q => $out_q,
        e => $out_e,
        i => $out_i,
        node => $out_node,
        argPeri => $out_argPeri,
        timePeri => $out_timePeri,
        hV => $hV,      # decorative, not used for anything
        objectName => $orb->{objectName},
        orbitId => $orb->{orbitId},
    };

    return $res;
}



# Autoload methods go after =cut, and are processed by the autosplit program.

1;
__END__

=head1 NAME

PS::MOPS::JPLEPH - Perl extension for generating ephemerides using
JPL's Rapid Comet and Asteroid Orbit Determination Process software.

=head1 SYNOPSIS

  use PS::MOPS::JPLEPH;

  # Single ephemeris.
  $results = jpleph_calcEphemeris(
    orbit => $orbit,
    epoch_mjd => $mjd,
    obscode => 568,
  );

  # Batch/multiple.
  $results = jpleph_calcEphemerides(
    instance => $inst,
    derivedObjectNames => \@list,
    epoch => $mjd,
    obscode => 568,
  );

=head1 DESCRIPTION

=head1 SEE ALSO

=head1 AUTHOR

Larry Denneau, E<lt>denneau@ifa.hawaii.edu<gt>

=head1 COPYRIGHT AND LICENSE

Copyright 2007 Institute for Astronomy, University of Hawaii

=cut
