#!/usr/bin/env python
"""
Unit Tets for the ephem module

Given an orbit with its covariance matrix and a list of dates (in MJD),
ephem.ephemerides() computes the position of that orbit at each of the input
dates. Computed positions are given as a 6-element numpy array of the form
    [RA, Dec, mag, MJD, Err1, Err2, PA]
Where all angles are given in radians. Err1, Err2 and PA define the error
ellipses on the sky plane.

Input and ground truth data from these tests was obtained XXX
"""
import math
import time

import numpy
import ssd

import unittest



# Constants
# Ppercent results have to agree within (1 = 100%)
THRESH = 0.00001


class TestEphemerides(object):
    orbitType = None
    orbitElements = None
    covar = None
    
    def testEphemerides(self):
        t0 = time.time()
        res = ssd.ephemerides(orbitElements=self.orbitElements,
                              covariance=self.covar,
                              orbitEpoch=self.orbitEpoch,
                              absMag=self.absMag,
                              gMag=self.gMag,
                              obsCode=self.obscode,
                              times=self.times)
        
        print('Speed: %.02f ephemerides/sec' \
              %(float(len(self.times)) / float(time.time() - t0)))
        
        # Check the results (RA, Dec, mag, t, smaa, smia, PA).
        # Remember that internally we sort times in reverse chronological order.
	# Also, now we spit out 1-sigma errors, not 3-sigma like we used to. 
	# Multiply res[4:8] by 3. to match the 3-sigma truth.
        truth = self.truth.tolist()
        truth.sort(lambda x, y: cmp(x[3], y[3]))
        truth.reverse()
        
        for i in range(len(self.times)):
            # Convert the 1-sigma errors to 3-sigma.
	    if(res[i][4] != -1):
	        res[i][4:8] *= 3.
            for j in range(len(res[i])):
                if(j == 2):
                    try:
                        assert(self._testMags(truth[i][j], res[i][j]))
                    except:
                        raise(AssertionError('%f differs from %f more than %f%%' \
                                             %(res[i][j], truth[i][j],
                                               THRESH * 100.)))
                else:
                    if(truth[i][j]):
                        difference = abs(truth[i][j] - res[i][j]) / \
                                     abs(truth[i][j])
                    else:
                        difference = abs(truth[i][j] - res[i][j])
                    try:
                        self.assert_(difference < THRESH)
                    except:
                        raise(AssertionError('%f differs from %f more than %f%%' \
                                             %(res[i][j], truth[i][j],
                                               THRESH * 100.)))
            # <-- end for
        # <-- end for
        return

    def _testMags(self, trueMag, testMag):
        """
        Check whether testMag is close enough to trueMag. In reality, given
        floating point roundoff errors when dealing with different compilers and
        architectures, we test mags only up to the 1/100th of a mag (which is
        still too precose given the big uncertainties in mag prediction).
        """
        trueMag = round(trueMag, 2)
        testMag = round(testMag, 2)
        return(round(abs(trueMag - testMag), 2) <= 0.01)
        


class Test2004VD17(TestEphemerides, unittest.TestCase):
    """
    Unit test for ssd.ephemerides() and asteroid 2004 VD17

    Since there will be some rounding errors in the floats we get, we simply ask
    that numbers agree up to a fraction of a percent (see THRESH above).

    This uses cometary elements.
    """
    def setUp(self):
        self.obscode = 568              # MPC observatory code.
        # Elements in MOPS order (q, e, i, node, peri, tp).
        # Angles in degrees.
        self.orbitElements = numpy.array([.6201443674443176,
                                          .5887584958766938,
                                          4.222821075488327,
                                          224.2500776466688,
                                          90.6966378352936,
                                          53260.8239252879284])
        self.orbitEpoch = 53458.
        self.absMag = 18.848
        self.gMag = .150
        self.covar = numpy.array([8.918907774759833E-9, -1.142810787742955E-8,
                                  1.668718350363152E-8, 2.209344232000432E-8,
                                  -4.46902785491213E-8, 9.630403388871625E-6,
                                  -2.773490120559766E-8, 6.580845530891771E-8, 
                                  5.770190812395115E-6, 1.784156547157364E-7, 
                                  -1.53148919367182E-8, 3.494932893594178E-8,
                                  -9.96775894966775E-9, -3.115210512717402E-7, 
                                  3.788909051693304E-7, -1.976329090803309E-9, 
                                  5.428725864025983E-9, -4.674612067825687E-7,
                                  3.141170291580454E-7, -3.067995347802948E-7, 
                                  1.163429608362596E-7])
        self.norm = None
        self.times = numpy.array([53522.])
        self.truth = numpy.array([[107.436398, 19.334377, 23.913811, 53522.,
                                   4.748166e-02, 5.081293e-02,
                                   5.082780e-02, 4.746573e-02, 8.612144e+01]])
        return

class TestSungrazer(TestEphemerides, unittest.TestCase):
    """
    Unit test for ssd.ephemerides() and a synthetic Sungrazer

    Since there will be some rounding errors in the floats we get, we simply ask
    that numbers agree up to a fraction of a percent (see THRESH above).

    This uses cometary elements.
    """
    def setUp(self):
        self.obscode = 568              # MPC observatory code.
        # Elements in MOPS order (q, e, i, node, peri, tp).
        # Angles in degrees.
        self.orbitElements = numpy.array([0.039588737047,
                                          0.975067984003,
                                          9.588468981146,
                                          346.985560396840,
                                          252.214254881082,
                                          54416.518883420147])
        self.orbitEpoch = 54466.
        self.absMag = 18.187
        self.gMag = .150
        self.covar = None
        self.times = numpy.array([54066.])
        self.truth = numpy.array([[51.13008, 32.36019, 22.756190, 54066,
                                   -1, -1, -1, -1, -1]])
        return
    

class Test2007ML24(TestEphemerides, unittest.TestCase):
    """
    Unit test for ssd.ephemerides() and asteroid 2007 ML24

    Since there will be some rounding errors in the floats we get, we simply ask
    that numbers agree up to a fraction of a percent (see THRESH above).

    This uses cometary elements.
    """
    def setUp(self):
        self.obscode = 568              # MPC observatory code.
        # Elements in MOPS order (q, e, i, node, peri, tp).
        # Angles in degrees.
        self.orbitElements = numpy.array([.4861175044921673,
                                          .3589113515751249,
                                          33.43335660067169,
                                          281.894025165122,
                                          201.4705897673642,
                                          54374.3473513098434])
        self.orbitEpoch = 54287.
        self.absMag = 19.196
        self.gMag = .150
        self.covar = numpy.array([1.195015243069159E-6, -1.969853219494382E-7,
                                  1.189918645194846E-6, -6.085414811993752E-6,
                                  5.784915091616561E-6, .0006924064864479496,
                                  4.374037687849605E-6, -3.541425501712141E-6,
                                  -.0003471893039529744, 3.239378687216248E-6,
                                  9.893053318363044E-5, -.0001049576345900015,
                                  -.00630981769507034, 1.273886745657654E-5,
                                  .000117157125743492, -8.785751420913023E-5,
                                  .0001833699531905975, .006337409607040429,
                                  .0002944435635015719, -.0004222988703965268,
                                  .0001487121334005637])
        self.norm = None
        self.times = numpy.array([53522., 53642., 53762., 53882., 54337.])
        self.truth = numpy.array([[357.82575, 11.08036, 21.059745, 53522.,
                                   394.44772, 1419.21976,
                                   1472.9881, 8.9311176, 74.471141],
                                  [161.4937917, 9.9221111, 20.205281, 53642.,
                                   755.29840, 2797.54543,
                                   2897.5572, 29.974551, -74.901699],
                                  [301.5729167, -16.15650, 21.284532, 53762.,
                                   626.64772, 576.632116,
                                   851.54349, 8.1017622, 42.619378],
                                  [72.725750, 24.655806, 19.525814, 53882.,
                                   1704.5788, 856.357697,
                                   1907.5933, 5.0458806, -26.674046],
                                  [123.9570417, 43.721806, 20.875203, 54337.,
                                   27.945567, 47.2880455,
                                   51.519263, 19.049395, -64.713127]])
        return



if(__name__ == '__main__'):
    unittest.main()
