Various tests of the use of nested lists of potentials

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)

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,

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()

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

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