# Tests of streamgapdf implementation, impulse tests moved to
# test_streamgapdf_impulse.py
import numpy
from scipy import integrate
import pytest
sdf_sanders15= None #so we can set this up and then use in other tests
sdf_sanders15_unp= None #so we can set this up and then use in other tests
sdfl_sanders15= None #so we can set this up and then use in other tests
sdfl_sanders15_unp= None #so we can set this up and then use in other tests

# Put seed in first function, so the seed gets set even if other test files
# were run first
def test_setupimpact_error():
    numpy.random.seed(1)
    #Imports
    from galpy.df import streamgapdf
    from galpy.orbit import Orbit
    from galpy.potential import LogarithmicHaloPotential
    from galpy.actionAngle import actionAngleIsochroneApprox
    from galpy.util import conversion #for unit conversions
    lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
    aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
    prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
    V0, R0= 220., 8.
    sigv= 0.365*(10./2.)**(1./3.) # km/s
    with pytest.raises(IOError) as excinfo:
        dum= streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                         leading=False,nTrackChunks=26,
                         nTrackIterations=1,
                         sigMeanOffset=4.5,
                         tdisrupt=10.88\
                             /conversion.time_in_Gyr(V0,R0),
                         Vnorm=V0,Rnorm=R0,
                         impactb=0.,
                         subhalovel=numpy.array([6.82200571,132.7700529,
                                                 149.4174464])/V0,
                         timpact=0.88/conversion.time_in_Gyr(V0,R0),
                         impact_angle=-2.34)
    # Should be including these:
    #                 GM=10.**-2.\
    #                     /conversion.mass_in_1010msol(V0,R0),
    #                 rs=0.625/R0)
    return None

def test_leadingwtrailingimpact_error():
    #Imports
    from galpy.df import streamgapdf
    from galpy.orbit import Orbit
    from galpy.potential import LogarithmicHaloPotential
    from galpy.actionAngle import actionAngleIsochroneApprox
    from galpy.util import conversion #for unit conversions
    lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
    aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
    prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
    V0, R0= 220., 8.
    sigv= 0.365*(10./2.)**(1./3.) # km/s
    with pytest.raises(ValueError) as excinfo:
        dum= streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                         leading=True,nTrackChunks=26,
                         nTrackIterations=1,
                         sigMeanOffset=4.5,
                         tdisrupt=10.88\
                             /conversion.time_in_Gyr(V0,R0),
                         Vnorm=V0,Rnorm=R0,
                         impactb=0.,
                         subhalovel=numpy.array([6.82200571,132.7700529,
                                                 149.4174464])/V0,
                         timpact=0.88/conversion.time_in_Gyr(V0,R0),
                         impact_angle=-2.34,
                         GM=10.**-2.\
                             /conversion.mass_in_1010msol(V0,R0),
                         rs=0.625/R0)
    return None

def test_trailingwleadingimpact_error():
    #Imports
    from galpy.df import streamgapdf
    from galpy.orbit import Orbit
    from galpy.potential import LogarithmicHaloPotential
    from galpy.actionAngle import actionAngleIsochroneApprox
    from galpy.util import conversion #for unit conversions
    lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
    aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
    prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
    V0, R0= 220., 8.
    sigv= 0.365*(10./2.)**(1./3.) # km/s
    with pytest.raises(ValueError) as excinfo:
        dum= streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                         leading=False,nTrackChunks=26,
                         nTrackIterations=1,
                         sigMeanOffset=4.5,
                         tdisrupt=10.88\
                             /conversion.time_in_Gyr(V0,R0),
                         Vnorm=V0,Rnorm=R0,
                         impactb=0.,
                         subhalovel=numpy.array([6.82200571,132.7700529,
                                                 149.4174464])/V0,
                         timpact=0.88/conversion.time_in_Gyr(V0,R0),
                         impact_angle=2.34,
                         GM=10.**-2.\
                             /conversion.mass_in_1010msol(V0,R0),
                         rs=0.625/R0)
    return None

#Exact setup from Section 5 of Sanders, Bovy, and Erkal (2015); should reproduce those results (which have been checked against a simulation)
def test_sanders15_setup():
    #Imports
    from galpy.df import streamdf, streamgapdf
    from galpy.orbit import Orbit
    from galpy.potential import LogarithmicHaloPotential
    from galpy.actionAngle import actionAngleIsochroneApprox
    from galpy.util import conversion #for unit conversions
    lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
    aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
    prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
    global sdf_sanders15
    V0, R0= 220., 8.
    sigv= 0.365*(10./2.)**(1./3.) # km/s
    sdf_sanders15= streamgapdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                               leading=False,nTrackChunks=26,
                               nTrackIterations=1,
                               sigMeanOffset=4.5,
                               tdisrupt=10.88\
                                   /conversion.time_in_Gyr(V0,R0),
                               Vnorm=V0,Rnorm=R0,
                               impactb=0.,
                               subhalovel=numpy.array([6.82200571,132.7700529,
                                                       149.4174464])/V0,
                               timpact=0.88/conversion.time_in_Gyr(V0,R0),
                               impact_angle=-2.34,
                               GM=10.**-2.\
                                   /conversion.mass_in_1010msol(V0,R0),
                               rs=0.625/R0)
    assert not sdf_sanders15 is None, 'sanders15 streamgapdf setup did not work'
    # Also setup the unperturbed model
    global sdf_sanders15_unp
    sdf_sanders15_unp= streamdf(sigv/V0,progenitor=prog_unp_peri,pot=lp,aA=aAI,
                               leading=False,nTrackChunks=26,
                               nTrackIterations=1,
                               sigMeanOffset=4.5,
                               tdisrupt=10.88\
                                   /conversion.time_in_Gyr(V0,R0),
                               Vnorm=V0,Rnorm=R0)
    assert not sdf_sanders15_unp is None, \
        'sanders15 unperturbed streamdf setup did not work'
    return None

def test_sanders15_leading_setup():
    #Imports
    from galpy.df import streamdf, streamgapdf
    from galpy.orbit import Orbit
    from galpy.potential import LogarithmicHaloPotential, PlummerPotential
    from galpy.actionAngle import actionAngleIsochroneApprox
    from galpy.util import conversion #for unit conversions
    lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
    aAI= actionAngleIsochroneApprox(pot=lp,b=0.8)
    prog_unp_peri= Orbit([2.6556151742081835,
                          0.2183747276300308,
                          0.67876510797240575,
                          -2.0143395648974671,
                          -0.3273737682604374,
                          0.24218273922966019])
    global sdfl_sanders15
    V0, R0= 220., 8.
    sigv= 0.365*(10./2.)**(1./3.) # km/s
    # Use a Potential object for the impact
    pp= PlummerPotential(amp=10.**-2.\
                             /conversion.mass_in_1010msol(V0,R0),
                         b=0.625/R0)
    import warnings
    from galpy.util import galpyWarning
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always",galpyWarning)
        sdfl_sanders15= streamgapdf(sigv/V0,progenitor=prog_unp_peri,
                                    pot=lp,aA=aAI,
                                    leading=True,nTrackChunks=26,
                                    nTrackChunksImpact=29,
                                    nTrackIterations=1,
                                    sigMeanOffset=4.5,
                                    tdisrupt=10.88\
                                        /conversion.time_in_Gyr(V0,R0),
                                    Vnorm=V0,Rnorm=R0,
                                    impactb=0.,
                                    subhalovel=numpy.array([49.447319,
                                                            116.179436,
                                                            155.104156])/V0,
                                    timpact=0.88/conversion.time_in_Gyr(V0,R0),
                                    impact_angle=2.09,
                                    subhalopot=pp,
                                    nKickPoints=290,
                                    deltaAngleTrackImpact=4.5,
                                    multi=True) # test multi
        # Should raise warning bc of deltaAngleTrackImpact, might raise others
        raisedWarning= False
        for wa in w:
            raisedWarning= (str(wa.message) == "WARNING: deltaAngleTrackImpact angle range large compared to plausible value")
            if raisedWarning: break
        assert raisedWarning,  'deltaAngleTrackImpact warning not raised when it should have been'
    assert not sdfl_sanders15 is None, 'sanders15 trailing streamdf setup did not work'
    # Also setup the unperturbed model
    global sdfl_sanders15_unp
    sdfl_sanders15_unp= streamdf(sigv/V0,progenitor=prog_unp_peri,
                                 pot=lp,aA=aAI,
                                 leading=True,nTrackChunks=26,
                                 nTrackIterations=1,
                                 sigMeanOffset=4.5,
                                 tdisrupt=10.88\
                                     /conversion.time_in_Gyr(V0,R0),
                                 Vnorm=V0,Rnorm=R0)
    assert not sdfl_sanders15_unp is None, \
        'sanders15 unperturbed streamdf setup did not work'
    return None

# Some very basic tests
def test_nTrackIterations():
    assert sdf_sanders15.nTrackIterations == 1, 'nTrackIterations should have been 1'
    return None
def test_nTrackChunks():
    assert sdf_sanders15._nTrackChunks == 26, 'nTrackChunks should have been 26'
    return None
def test_deltaAngleTrackImpact():
    assert numpy.fabs(sdf_sanders15._deltaAngleTrackImpact-4.31) < 0.01, 'deltaAngleTrackImpact should have been ~4.31'
    return None

# Tests of the track near the impact
def test_trackNearImpact():
    # Sanity checks against numbers taken from plots of the simulation
    # Make sure we're near 14.5
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[14,0]*sdf_sanders15._ro
                      -14.5) < 0.2, '14th point along track near the impact is not near 14.5 kpc'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[14,1]*sdf_sanders15._vo
                      -80) < 3., 'Point along the track near impact near R=14.5 does not have the correct radial velocity'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[14,2]*sdf_sanders15._vo
                      -220.) < 3., 'Point along the track near impact near R=14.5 does not have the correct tangential velocity'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[14,3]*sdf_sanders15._ro
                      -0.) < 1., 'Point along the track near impact near R=14.5 does not have the correct vertical height'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[14,4]*sdf_sanders15._vo
                      -200.) < 5., 'Point along the track near impact near R=14.5 does not have the correct vertical velocity'
    # Another one!
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,0]*sdf_sanders15._ro
                      -16.25) < 0.2, '27th point along track near the impact is not near 16.25 kpc'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,1]*sdf_sanders15._vo
                      +130) < 3., 'Point along the track near impact near R=16.25 does not have the correct radial velocity'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,2]*sdf_sanders15._vo
                      -200.) < 3., 'Point along the track near impact near R=16.25 does not have the correct tangential velocity'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,3]*sdf_sanders15._ro
                      +12.) < 1., 'Point along the track near impact near R=16.25 does not have the correct vertical height'
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,4]*sdf_sanders15._vo
                      -25.) < 5., 'Point along the track near impact near R=16.25 does not have the correct vertical velocity'   
    assert numpy.fabs(sdf_sanders15._gap_ObsTrack[27,5]-1.2) < .2, 'Point along the track near impact near R=16.25 does not have the correct azimuth'   
    return None

def test_interpolatedTrackNearImpact():
    # Sanity checks against numbers taken from plots of the simulation
    # Make sure we're near X=-10.9
    theta= 2.7
    assert numpy.fabs(sdf_sanders15._kick_interpTrackX(theta)*sdf_sanders15._ro
                      +10.9) < 0.2, 'Point along track near the impact at theta=2.7 is not near X=-10.9 kpc'
    assert numpy.fabs(sdf_sanders15._kick_interpTrackY(theta)*sdf_sanders15._ro
                      -6.) < 0.5, 'Point along track near the impact at theta=2.7 is not near Y=6. kpc'
    assert numpy.fabs(sdf_sanders15._kick_interpTrackZ(theta)*sdf_sanders15._ro
                      +5.) < 0.5, 'Point along track near the impact at theta=2.7 is not near Z=5. kpc'
    assert numpy.fabs(sdf_sanders15._kick_interpTrackvX(theta)*sdf_sanders15._vo
                      +180.) < 5, 'Point along track near the impact at theta=2.7 is not near vX=-180 km/s'
    assert numpy.fabs(sdf_sanders15._kick_interpTrackvY(theta)*sdf_sanders15._vo
                      +190.) < 5., 'Point along track near the impact at theta=2.7 is not near vY=190 km/s'
    assert numpy.fabs(sdf_sanders15._kick_interpTrackvZ(theta)*sdf_sanders15._vo
                      -170.) < 5., 'Point along track near the impact at theta=2.7 is not near vZ=170 km/s'
    return None

# Test the calculation of the kicks in dv
def test_kickdv():
    # Closest one to the impact point, should be close to zero
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle))
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_deltav[tIndx]*sdf_sanders15._vo) < 0.3), 'Kick near the impact point not close to zero'
    # The peak, size and location
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_deltav[:,0]*sdf_sanders15._vo))-0.35) < 0.06, 'Peak dvx incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_deltav[:,0]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle < 0., 'Location of peak dvx incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_deltav[:,1]*sdf_sanders15._vo))-0.35) < 0.06, 'Peak dvy incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_deltav[:,1]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle > 0., 'Location of peak dvy incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_deltav[:,2]*sdf_sanders15._vo))-1.8) < 0.06, 'Peak dvz incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_deltav[:,2]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle > 0., 'Location of peak dvz incorrect'
    # Close to zero far from impact point
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle-1.5))
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_deltav[tIndx]*sdf_sanders15._vo) < 0.3), 'Kick far the impact point not close to zero'
    return None

# Test the calculation of the kicks in dO
def test_kickdO():
    from galpy.util import conversion
    # Closest one to the impact point, should be close to zero
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle))
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_dOap[tIndx,:3]*conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)) < 0.03), 'Kick near the impact point not close to zero'
    # The peak, size and location
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_dOap[:,0]*conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)))-0.085) < 0.01, 'Peak dOR incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_dOap[:,0])]-sdf_sanders15._impact_angle < 0., 'Location of peak dOR incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_dOap[:,1]*conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)))-0.07) < 0.01, 'Peak dOp incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_dOap[:,1])]-sdf_sanders15._impact_angle < 0., 'Location of peak dvy incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(sdf_sanders15._kick_dOap[:,2]*conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)))-0.075) < 0.01, 'Peak dOz incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(sdf_sanders15._kick_dOap[:,2])]-sdf_sanders15._impact_angle < 0., 'Location of peak dOz incorrect'
    # Close to zero far from impact point
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle-1.5))
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_dOap[tIndx,:3]*conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)) < 0.03), 'Kick far the impact point not close to zero'
    return None

def test_kickda():
    # All angle kicks should be small, just test that they are smaller than dO/O close to the impact
    nIndx= numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack-sdf_sanders15._impact_angle) < 0.75
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_dOap[nIndx,3:]) < 2.*(numpy.fabs(sdf_sanders15._kick_dOap[nIndx,:3]/sdf_sanders15._progenitor_Omega))), 'angle kicks not smaller than the frequency kicks'
    return None

# Test the interpolation of the kicks
def test_interpKickdO():
    from galpy.util import conversion
    freqConv= conversion.freq_in_Gyr(sdf_sanders15._vo,sdf_sanders15._ro)
    # Bunch of spot checks at some interesting angles
    # Impact angle
    theta= sdf_sanders15._impact_angle
    assert numpy.fabs(sdf_sanders15._kick_interpdOpar(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp0(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp1(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOr(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOp(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOz(theta)*freqConv) < 10.**-4., 'Frequency kick at the impact point is not zero'
    # random one
    theta= sdf_sanders15._impact_angle-0.25
    assert numpy.fabs(sdf_sanders15._kick_interpdOpar(theta)*freqConv+0.07) < 0.002, 'Frequency kick near the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp0(theta)*freqConv) < 0.002, 'Frequency kick near the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp1(theta)*freqConv) < 0.003, 'Frequency kick near the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOr(theta)*freqConv-0.05) < 0.01, 'Frequency kick near the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOp(theta)*freqConv-0.035) < 0.01, 'Frequency kick near the impact point is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOz(theta)*freqConv-0.04) < 0.01, 'Frequency kick near the impact point is not zero'
    # One beyond ._deltaAngleTrackImpact
    theta= sdf_sanders15._deltaAngleTrackImpact+0.1
    assert numpy.fabs(sdf_sanders15._kick_interpdOpar(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp0(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOperp1(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOr(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOp(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdOz(theta)*freqConv) < 10.**-16., 'Frequency kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdar(theta)) < 10.**-16., 'Angle kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdap(theta)) < 10.**-16., 'Angle kick beyond ._deltaAngleTrackImpact is not zero'
    assert numpy.fabs(sdf_sanders15._kick_interpdaz(theta)) < 10.**-16., 'Angle kick beyond ._deltaAngleTrackImpact is not zero'
    return None

def test_interpKickda():
    thetas= numpy.linspace(-0.75,0.75,10)+sdf_sanders15._impact_angle
    assert numpy.all(numpy.fabs(sdf_sanders15._kick_interpdar(thetas)) \
                         < 2.*numpy.fabs(sdf_sanders15._kick_interpdOr(thetas)/sdf_sanders15._progenitor_Omegar)), 'Interpolated angle kick not everywhere smaller than the frequency kick after one period'
    return None

# Test the sampling of present-day perturbed points based on the model
def test_sample():
    # Sample stars from the model and compare them to the stream
    xv_mock_per= sdf_sanders15.sample(n=100000,xy=True).T
    # Rough gap-density check
    ingap= numpy.sum((xv_mock_per[:,0]*sdf_sanders15._ro > 4.)\
                         *(xv_mock_per[:,0]*sdf_sanders15._ro < 5.))
    edgegap= numpy.sum((xv_mock_per[:,0]*sdf_sanders15._ro > 1.)\
                         *(xv_mock_per[:,0]*sdf_sanders15._ro < 2.))
    outgap= numpy.sum((xv_mock_per[:,0]*sdf_sanders15._ro > -2.5)\
                         *(xv_mock_per[:,0]*sdf_sanders15._ro < -1.5))
    assert numpy.fabs(ingap/float(edgegap)-0.015/0.05) < 0.05, 'gap density versus edge of the gap is incorect'
    assert numpy.fabs(ingap/float(outgap)-0.015/0.02) < 0.2, 'gap density versus outside of the gap is incorect'
    # Test track of the stream
    tIndx= (xv_mock_per[:,0]*sdf_sanders15._ro > 4.)\
        *(xv_mock_per[:,0]*sdf_sanders15._ro < 5.)\
        *(xv_mock_per[:,1]*sdf_sanders15._ro < 5.)
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,1])*sdf_sanders15._ro+12.25) < 0.1, 'Location of mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,2])*sdf_sanders15._ro-3.8) < 0.1, 'Location of mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,3])*sdf_sanders15._vo-255.) < 2., 'Location of mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,4])*sdf_sanders15._vo-20.) < 2., 'Location of mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,5])*sdf_sanders15._vo+185.) < 2., 'Location of mock track is incorrect near the gap'
    return None

# Test the sampling of present-day perturbed-unperturbed points
# (like in the paper)
def test_sample_offset():
    # Sample stars from the model and compare them to the stream
    numpy.random.seed(1)
    xv_mock_per= sdf_sanders15.sample(n=100000,xy=True).T
    numpy.random.seed(1) # should give same points
    xv_mock_unp= sdf_sanders15_unp.sample(n=100000,xy=True).T
    # Test perturbation as a function of unperturbed X
    tIndx= (xv_mock_unp[:,0]*sdf_sanders15._ro > 0.)\
        *(xv_mock_unp[:,0]*sdf_sanders15._ro < 1.)\
        *(xv_mock_unp[:,1]*sdf_sanders15._ro < 5.)
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,0]-xv_mock_unp[tIndx,0])*sdf_sanders15._ro+0.65) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,1]-xv_mock_unp[tIndx,1])*sdf_sanders15._ro-0.1) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,2]-xv_mock_unp[tIndx,2])*sdf_sanders15._ro-0.4) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,3]-xv_mock_unp[tIndx,3])*sdf_sanders15._vo) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,4]-xv_mock_unp[tIndx,4])*sdf_sanders15._vo+7.) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,5]-xv_mock_unp[tIndx,5])*sdf_sanders15._vo-4.) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    return None

# Test the sampling of present-day perturbed-unperturbed points
# (like in the paper, but for the leading stream impact)
def test_sample_offset_leading():
    # Sample stars from the model and compare them to the stream
    numpy.random.seed(1)
    xv_mock_per= sdfl_sanders15.sample(n=100000,xy=True).T
    numpy.random.seed(1) # should give same points
    xv_mock_unp= sdfl_sanders15_unp.sample(n=100000,xy=True).T
    # Test perturbation as a function of unperturbed X
    tIndx= (xv_mock_unp[:,0]*sdfl_sanders15._ro > 13.)\
        *(xv_mock_unp[:,0]*sdfl_sanders15._ro < 14.)\
        *(xv_mock_unp[:,1]*sdfl_sanders15._ro > 5.)
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,0]-xv_mock_unp[tIndx,0])*sdfl_sanders15._ro+0.5) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,1]-xv_mock_unp[tIndx,1])*sdfl_sanders15._ro-0.3) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,2]-xv_mock_unp[tIndx,2])*sdfl_sanders15._ro-0.45) < 0.1, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,3]-xv_mock_unp[tIndx,3])*sdfl_sanders15._vo+2.) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,4]-xv_mock_unp[tIndx,4])*sdfl_sanders15._vo+7.) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    assert numpy.fabs(numpy.median(xv_mock_per[tIndx,5]-xv_mock_unp[tIndx,5])*sdfl_sanders15._vo-6.) < 0.5, 'Location of perturbed mock track is incorrect near the gap'
    return None

# Tests of the density and meanOmega functions

def test_pOparapar():
    #Test that integrating pOparapar gives density_par
    dens_frompOpar_close=\
        integrate.quad(lambda x: sdf_sanders15.pOparapar(x,0.3),
                       sdf_sanders15._meandO\
                           -10.*numpy.sqrt(sdf_sanders15._sortedSigOEig[2]),
                       sdf_sanders15._meandO\
                           +10.*numpy.sqrt(sdf_sanders15._sortedSigOEig[2]))[0]
    # This is actually in the gap!
    dens_fromOpar_half=\
        integrate.quad(lambda x: sdf_sanders15.pOparapar(x,2.6),
                       sdf_sanders15._meandO\
                           -10.*numpy.sqrt(sdf_sanders15._sortedSigOEig[2]),
                       sdf_sanders15._meandO\
                           +10.*numpy.sqrt(sdf_sanders15._sortedSigOEig[2]))[0]
    assert numpy.fabs(dens_fromOpar_half/dens_frompOpar_close-sdf_sanders15.density_par(2.6)/sdf_sanders15.density_par(0.3)) < 10.**-4., 'density from integrating pOparapar not equal to that from density_par for Sanders15 stream'
    return None

def test_density_apar_approx():
    # Test that the approximate density agrees with the direct integration
    # Need to do this relatively to another density, because there is an
    # overall offset
    apar= 2.6
    assert numpy.fabs(sdf_sanders15.density_par(apar,approx=False)/sdf_sanders15.density_par(apar,approx=True)/sdf_sanders15.density_par(0.3,approx=False)*sdf_sanders15.density_par(0.3,approx=True)-1.) < 10.**-3., 'Approximate density does not agree with direct integration'
    apar= 2.3
    assert numpy.fabs(sdf_sanders15.density_par(apar,approx=False)/sdf_sanders15.density_par(apar,approx=True)/sdf_sanders15.density_par(0.3,approx=False)*sdf_sanders15.density_par(0.3,approx=True)-1.) < 10.**-3., 'Approximate density does not agree with direct integration'
    return None

def test_density_apar_approx_higherorder():
    # Test that the approximate density agrees with the direct integration
    # Need to do this relatively to another density, because there is an
    # overall offset
    apar= 2.6
    assert numpy.fabs(sdf_sanders15.density_par(apar,approx=False)/sdf_sanders15.density_par(apar,approx=True,higherorder=True)/sdf_sanders15.density_par(0.3,approx=False)*sdf_sanders15.density_par(0.3,approx=True,higherorder=True)-1.) < 10.**-3., 'Approximate density does not agree with direct integration'
    apar= 2.3
    assert numpy.fabs(sdf_sanders15.density_par(apar,approx=False)/sdf_sanders15.density_par(apar,approx=True,higherorder=True)/sdf_sanders15.density_par(0.3,approx=False)*sdf_sanders15.density_par(0.3,approx=True,higherorder=True)-1.) < 10.**-3., 'Approximate density does not agree with direct integration'
    return None

def test_minOpar():
    # Test that for Opar < minOpar, p(Opar,apar) is in fact zero!
    apar= 0.3
    dO= 10.**-4.
    assert numpy.fabs(sdf_sanders15.pOparapar(sdf_sanders15.minOpar(apar)-dO,
                                              apar)) < 10.**-16., 'Probability for Opar < minOpar is not zero'
    apar= 2.6
    dO= 10.**-4.
    assert numpy.fabs(sdf_sanders15.pOparapar(sdf_sanders15.minOpar(apar)-dO,
                                              apar)) < 10.**-16., 'Probability for Opar < minOpar is not zero'
    return None

def test_meanOmega_approx():
    # Test that the approximate meanOmega agrees with the direct integration
    # Need to do this relatively to another density, because there is an
    # overall offset
    apar= 2.6
    assert numpy.fabs(sdf_sanders15.meanOmega(apar,approx=False,oned=True)/sdf_sanders15.meanOmega(apar,approx=True,oned=True)-1.) < 10.**-3., 'Approximate meanOmega does not agree with direct integration'
    apar= 2.3
    assert numpy.fabs(sdf_sanders15.meanOmega(apar,approx=False,oned=True)/sdf_sanders15.meanOmega(apar,approx=True,oned=True)-1.) < 10.**-3., 'Approximate meanOmega does not agree with direct integration'
    return None

def test_meanOmega_approx_higherorder():
    # Test that the approximate meanOmega agrees with the direct integration
    # Need to do this relatively to another density, because there is an
    # overall offset
    apar= 2.6
    assert numpy.fabs(sdf_sanders15.meanOmega(apar,approx=False,oned=True)/sdf_sanders15.meanOmega(apar,approx=True,higherorder=True,oned=True)-1.) < 10.**-3., 'Approximate meanOmega does not agree with direct integration'
    apar= 2.3
    assert numpy.fabs(sdf_sanders15.meanOmega(apar,approx=False,oned=True)/sdf_sanders15.meanOmega(apar,approx=True,higherorder=True,oned=True)-1.) < 10.**-3., 'Approximate meanOmega does not agree with direct integration'
    return None

def test_hernquist():
    # Test that Hernquist kicks are similar to Plummer kicks, but are
    # different in understood ways (...)
    from galpy.util import conversion
    # Switch to Hernquist
    V0, R0= 220., 8.
    impactb=0.
    subhalovel=numpy.array([6.82200571,132.7700529,
                            149.4174464])/V0
    impact_angle=-2.34
    GM=10.**-2./conversion.mass_in_1010msol(V0,R0)
    rs=0.625/R0
    sdf_sanders15._determine_deltav_kick(impact_angle,impactb,subhalovel,
                                         GM,rs,None,
                                         3,True)
    hernquist_kicks= sdf_sanders15._kick_deltav
    # Back to Plummer
    sdf_sanders15._determine_deltav_kick(impact_angle,impactb,subhalovel,
                                         GM,rs,None,
                                         3,False)
    # Repeat some of the deltav tests from above
    # Closest one to the impact point, should be close to zero
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle))
    assert numpy.all(numpy.fabs(hernquist_kicks[tIndx]*sdf_sanders15._vo) < 0.4), 'Kick near the impact point not close to zero for Hernquist'
    # The peak, size and location
    # Peak should be slightly less (guessed these correct!)
    assert numpy.fabs(numpy.amax(numpy.fabs(hernquist_kicks[:,0]*sdf_sanders15._vo))-0.25) < 0.06, 'Peak dvx incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(hernquist_kicks[:,0]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle < 0., 'Location of peak dvx incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(hernquist_kicks[:,1]*sdf_sanders15._vo))-0.25) < 0.06, 'Peak dvy incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(hernquist_kicks[:,1]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle > 0., 'Location of peak dvy incorrect'
    assert numpy.fabs(numpy.amax(numpy.fabs(hernquist_kicks[:,2]*sdf_sanders15._vo))-1.3) < 0.06, 'Peak dvz incorrect'
    assert sdf_sanders15._kick_interpolatedThetasTrack[numpy.argmax(hernquist_kicks[:,2]*sdf_sanders15._vo)]-sdf_sanders15._impact_angle > 0., 'Location of peak dvz incorrect'
    # Close to zero far from impact point
    tIndx= numpy.argmin(numpy.fabs(sdf_sanders15._kick_interpolatedThetasTrack\
                                       -sdf_sanders15._impact_angle-1.5))
    assert numpy.all(numpy.fabs(hernquist_kicks[tIndx]*sdf_sanders15._vo) < 0.3), 'Kick far the impact point not close to zero'
    return None

def test_determine_deltav_valueerrort():
    # Test that modeling leading (trailing) impact for trailing (leading) arm
    # raises a ValueError when using _determine_deltav_kick
    from galpy.util import conversion
    # Switch to Hernquist
    V0, R0= 220., 8.
    impactb=0.
    subhalovel=numpy.array([6.82200571,132.7700529,
                            149.4174464])/V0
    impact_angle=-2.34
    GM=10.**-2./conversion.mass_in_1010msol(V0,R0)
    rs=0.625/R0
    # Can't do minus impact angle!
    with pytest.raises(ValueError) as excinfo:
        sdf_sanders15._determine_deltav_kick(-impact_angle,impactb,subhalovel,
                                              GM,rs,None,
                                              3,True)
    return None

# Test the routine that rotates vectors to an arbitrary vector
def test_rotate_to_arbitrary_vector():
    from galpy.df.streamgapdf import _rotate_to_arbitrary_vector
    tol= -10.
    v= numpy.array([[1.,0.,0.]])
    # Rotate to 90 deg off
    ma= _rotate_to_arbitrary_vector(v,[0,1.,0])
    assert numpy.fabs(ma[0,0,1]+1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,0]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,2]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    # Rotate to 90 deg off
    ma= _rotate_to_arbitrary_vector(v,[0,0,1.])
    assert numpy.fabs(ma[0,0,2]+1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,0]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,1]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    # Rotate to same should be unit matrix
    ma= _rotate_to_arbitrary_vector(v,v[0])
    assert numpy.all(numpy.fabs(numpy.diag(ma[0])-1.) < 10.**tol), \
        'Rotation matrix to same vector is not unity'
    assert numpy.fabs(numpy.sum(ma**2.)-3.)< 10.**tol, \
        'Rotation matrix to same vector is not unity'
    # Rotate to -same should be -unit matrix
    ma= _rotate_to_arbitrary_vector(v,-v[0])
    assert numpy.all(numpy.fabs(numpy.diag(ma[0])+1.) < 10.**tol), \
        'Rotation matrix to minus same vector is not minus unity'
    assert numpy.fabs(numpy.sum(ma**2.)-3.)< 10.**tol, \
        'Rotation matrix to minus same vector is not minus unity'
    return None

# Test that the rotation routine works for multiple vectors
def test_rotate_to_arbitrary_vector_multi():
    from galpy.df.streamgapdf import _rotate_to_arbitrary_vector
    tol= -10.
    v= numpy.array([[1.,0.,0.],[0.,1.,0.]])
    # Rotate to 90 deg off
    ma= _rotate_to_arbitrary_vector(v,[0,0,1.])
    assert numpy.fabs(ma[0,0,2]+1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,0]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,1]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    # 2nd
    assert numpy.fabs(ma[1,1,2]+1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,2,1]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,0,0]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,0,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,0,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,1,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,1,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,2,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[1,2,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    return None

# Test the inverse of the routine that rotates vectors to an arbitrary vector
def test_rotate_to_arbitrary_vector_inverse():
    from galpy.df.streamgapdf import _rotate_to_arbitrary_vector
    tol= -10.
    v= numpy.array([[1.,0.,0.]])
    # Rotate to random vector and back
    a= numpy.random.uniform(size=3)
    a/= numpy.sqrt(numpy.sum(a**2.))
    ma= _rotate_to_arbitrary_vector(v,a)
    ma_inv= _rotate_to_arbitrary_vector(v,a,inv=True)
    ma= numpy.dot(ma[0],ma_inv[0])
    assert numpy.all(numpy.fabs(ma-numpy.eye(3)) < 10.**tol), 'Inverse rotation matrix incorrect'
    return None

# Test that rotating to vy in particular works as expected
def test_rotation_vy():
    from galpy.df.streamgapdf import _rotation_vy
    tol= -10.
    v= numpy.array([[1.,0.,0.]])
    # Rotate to 90 deg off
    ma= _rotation_vy(v)
    assert numpy.fabs(ma[0,0,1]+1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,0]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,2]-1.) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,0,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,1,2]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,0]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'
    assert numpy.fabs(ma[0,2,1]) < 10.**tol, 'Rotation matrix to 90 deg off incorrect'

