From 1d2850f50e7b721b2db62184b881beff54154e0f Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 15 Mar 2018 13:22:53 -0400 Subject: [PATCH] Various tests of the use of nested lists of potentials --- tests/test_actionAngle.py | 13 +++++++++---- tests/test_orbit.py | 11 +++++++---- tests/test_potential.py | 17 ++++++++++++++--- tests/test_qdf.py | 4 +++- tests/test_streamdf.py | 3 ++- 5 files changed, 35 insertions(+), 13 deletions(-) diff --git a/tests/test_actionAngle.py b/tests/test_actionAngle.py index e350a652c..f9c5b8cc7 100644 --- a/tests/test_actionAngle.py +++ b/tests/test_actionAngle.py @@ -480,7 +480,7 @@ def test_actionAngleAdiabatic_basic_actions_gamma0(): from galpy.actionAngle import actionAngleAdiabatic from galpy.orbit import Orbit from galpy.potential import MWPotential - aAA= actionAngleAdiabatic(pot=MWPotential,gamma=0.) + aAA= actionAngleAdiabatic(pot=[MWPotential[0],MWPotential[1:]],gamma=0.) #circular orbit R,vR,vT,phi= 1.,0.,1.,2. js= aAA(Orbit([R,vR,vT,phi])) @@ -508,7 +508,8 @@ def test_actionAngleAdiabatic_basic_actions_c(): from galpy.actionAngle import actionAngleAdiabatic from galpy.orbit import Orbit from galpy.potential import MWPotential - aAA= actionAngleAdiabatic(pot=MWPotential,c=True) + # test nested list of potentials + aAA= actionAngleAdiabatic(pot=[MWPotential[0],MWPotential[1:]],c=True) #circular orbit R,vR,vT,z,vz= 1.,0.,1.,0.,0. js= aAA(R,vR,vT,z,vz) @@ -872,7 +873,9 @@ def test_actionAngleStaeckel_basic_actions_u0(): from galpy.actionAngle import actionAngleStaeckel from galpy.orbit import Orbit from galpy.potential import MWPotential - aAS= actionAngleStaeckel(pot=MWPotential,delta=0.71,c=False,useu0=True) + # test nested list of potentials + aAS= actionAngleStaeckel(pot=[MWPotential[0],MWPotential[1:]], + delta=0.71,c=False,useu0=True) #circular orbit R,vR,vT,z,vz= 1.,0.,1.,0.,0. js= aAS(R,vR,vT,z,vz) @@ -890,7 +893,9 @@ def test_actionAngleStaeckel_basic_actions_u0_c(): from galpy.actionAngle import actionAngleStaeckel from galpy.orbit import Orbit from galpy.potential import MWPotential - aAS= actionAngleStaeckel(pot=MWPotential,delta=0.71,c=True,useu0=True) + # test nested list of potentials + aAS= actionAngleStaeckel(pot=[MWPotential[0],MWPotential[1:]], + delta=0.71,c=True,useu0=True) #circular orbit R,vR,vT,z,vz= 1.,0.,1.,0.,0. js= aAS(R,vR,vT,z,vz) diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 754f2ad7b..8663d6aa4 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -53,7 +53,8 @@ mockFlatSolidBodyRotationSpiralArmsPotential, \ mockFlatSolidBodyRotationPlanarSpiralArmsPotential, \ triaxialLogarithmicHaloPotential, \ - testorbitHenonHeilesPotential + testorbitHenonHeilesPotential, \ + nestedListPotential _TRAVIS= bool(os.getenv('TRAVIS')) if not _TRAVIS: _QUICKTEST= True #Run a more limited set of tests @@ -123,6 +124,7 @@ def test_energy_jacobi_conservation(): pots.append('mockFlatSolidBodyRotationPlanarSpiralArmsPotential') pots.append('triaxialLogarithmicHaloPotential') pots.append('testorbitHenonHeilesPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -301,7 +303,7 @@ def test_energy_jacobi_conservation(): if isinstance(tp,testMWPotential) \ or isinstance(tp,testplanarMWPotential): o.integrate(ttimes, - [tmp.toPlanar() for tmp in tp._potlist], + potential.toPlanarPotential(tp._potlist), method=integrator) else: o.integrate(ttimes,ptp,method=integrator) @@ -347,7 +349,7 @@ def test_energy_jacobi_conservation(): if isinstance(tp,testMWPotential) \ or isinstance(tp,testplanarMWPotential): o.integrate(ttimes, - [tmp.toPlanar() for tmp in tp._potlist], + potential.toPlanarPotential(tp._potlist), method=integrator) else: o.integrate(ttimes,ptp,method=integrator) @@ -516,6 +518,7 @@ def test_liouville_planar(): pots.append('mockFlatSolidBodyRotationSpiralArmsPotential') pots.append('triaxialLogarithmicHaloPotential') pots.append('testorbitHenonHeilesPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -564,7 +567,7 @@ def test_liouville_planar(): #Calculate the Jacobian d x / d x if hasattr(tp,'_potlist'): if isinstance(tp,testMWPotential): - plist= [tmp.toPlanar() for tmp in tp._potlist] + plist= potential.toPlanarPotential(tp._potlist) else: plist= tp._potlist o.integrate_dxdv([1.,0.,0.,0.],ttimes,plist, diff --git a/tests/test_potential.py b/tests/test_potential.py index 2f9e9873c..478548e28 100644 --- a/tests/test_potential.py +++ b/tests/test_potential.py @@ -139,6 +139,7 @@ def test_forceAsDeriv_potential(): pots.append('mockDehnenSmoothBarPotentialTm5') pots.append('SolidBodyRotationSpiralArmsPotential') pots.append('triaxialLogarithmicHaloPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -299,6 +300,7 @@ def test_2ndDeriv_potential(): pots.append('mockDehnenSmoothBarPotentialTm5') pots.append('SolidBodyRotationSpiralArmsPotential') pots.append('triaxialLogarithmicHaloPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -515,6 +517,7 @@ def test_poisson_potential(): pots.append('DehnenSmoothDehnenBarPotential') pots.append('SolidBodyRotationSpiralArmsPotential') pots.append('triaxialLogarithmicHaloPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -539,6 +542,7 @@ def test_poisson_potential(): tol['rotatingSpiralArmsPotential']= -3 tol['specialSpiralArmsPotential']= -4 tol['SolidBodyRotationSpiralArmsPotential']= -2.9 #these are more difficult + tol['nestedListPotential']= -3 #these are more difficult #tol['RazorThinExponentialDiskPotential']= -6. for p in pots: #if not 'NFW' in p: continue #For testing the test @@ -629,6 +633,7 @@ def test_evaluateAndDerivs_potential(): pots.append('mockDehnenSmoothBarPotentialTm1') pots.append('mockDehnenSmoothBarPotentialTm5') pots.append('triaxialLogarithmicHaloPotential') + pots.append('nestedListPotential') rmpots= ['Potential','MWPotential','MWPotential2014', 'MovingObjectPotential', 'interpRZPotential', 'linearPotential', 'planarAxiPotential', @@ -2499,7 +2504,7 @@ def __init__(self): from galpy.potential import Potential, \ evaluatePotentials, evaluateRforces, evaluatezforces, evaluatephiforces, \ evaluateR2derivs, evaluatez2derivs, evaluateRzderivs, \ - evaluateDensities + evaluateDensities, _isNonAxi from galpy.potential import planarPotential, \ evaluateplanarPotentials, evaluateplanarRforces, evaluateplanarphiforces, \ evaluateplanarR2derivs @@ -2508,7 +2513,7 @@ class testMWPotential(Potential): def __init__(self,potlist=MWPotential): self._potlist= potlist Potential.__init__(self,amp=1.) - self.isNonAxi= True^numpy.prod([True^p.isNonAxi for p in self._potlist]) + self.isNonAxi= _isNonAxi(self._potlist) return None def _evaluate(self,R,z,phi=0,t=0,dR=0,dphi=0): return evaluatePotentials(self._potlist,R,z,phi=phi,t=t, @@ -2546,7 +2551,7 @@ def __init__(self,potlist=MWPotential): self._potlist= [p.toPlanar() for p in potlist if isinstance(p,Potential)] self._potlist.extend([p for p in potlist if isinstance(p,planarPotential)]) planarPotential.__init__(self,amp=1.) - self.isNonAxi= True^numpy.prod([True^p.isNonAxi for p in self._potlist]) + self.isNonAxi= _isNonAxi(self._potlist) return None def _evaluate(self,R,phi=0,t=0,dR=0,dphi=0): return evaluateplanarPotentials(self._potlist,R,phi=phi,t=t) @@ -2843,3 +2848,9 @@ def __init__(self): def OmegaP(self): # Non-axi, so need to set this to zero for Jacobi return 0. +class nestedListPotential(testMWPotential): + def __init__(self): + testMWPotential.__init__(self, + potlist=[potential.MWPotential2014,potential.SpiralArmsPotential()]) + def OmegaP(self): + return self._potlist[1].OmegaP() diff --git a/tests/test_qdf.py b/tests/test_qdf.py index 38e00b917..e259ecc02 100644 --- a/tests/test_qdf.py +++ b/tests/test_qdf.py @@ -20,8 +20,10 @@ def test_meanvR_adiabatic_gl(): def test_meanvR_adiabatic_mc(): numpy.random.seed(1) + # test nested list of potentials qdf= quasiisothermaldf(1./4.,0.2,0.1,1.,1., - pot=MWPotential,aA=aAA,cutcounter=True) + pot=[MWPotential[0],MWPotential[1:]], + aA=aAA,cutcounter=True) #In the mid-plane assert numpy.fabs(qdf.meanvR(0.9,0.,mc=True)) < 0.01, "qdf's meanvr is not equal to zero for adiabatic approx." #higher up diff --git a/tests/test_streamdf.py b/tests/test_streamdf.py index 7310945e5..db0b8e8fe 100644 --- a/tests/test_streamdf.py +++ b/tests/test_streamdf.py @@ -1261,8 +1261,9 @@ def test_fardalpot_trackaa(): from galpy.potential import IsochronePotential, FlattenedPowerPotential from galpy.actionAngle import actionAngleIsochroneApprox from galpy.util import bovy_conversion #for unit conversions + # test nested list of potentials pot= [IsochronePotential(b=0.8,normalize=0.8), - FlattenedPowerPotential(alpha=-0.7,q=0.6,normalize=0.2)] + [FlattenedPowerPotential(alpha=-0.7,q=0.6,normalize=0.2)]] aAI= actionAngleIsochroneApprox(pot=pot,b=0.9) obs= Orbit([1.10, 0.32, -1.15, 1.10, 0.31, 3.0]) sigv= 1.3 #km/s