diff --git a/galpy/actionAngle/actionAngleAdiabaticGrid.py b/galpy/actionAngle/actionAngleAdiabaticGrid.py index 8d5d70277..776146b45 100644 --- a/galpy/actionAngle/actionAngleAdiabaticGrid.py +++ b/galpy/actionAngle/actionAngleAdiabaticGrid.py @@ -37,7 +37,7 @@ def __init__( nEr=31, nLz=31, numcores=1, - **kwargs + **kwargs, ): """ Initialize an actionAngleAdiabaticGrid object @@ -103,7 +103,7 @@ def __init__( numpy.ones(len(thisRs)), # these two r dummies numpy.zeros(len(thisRs)), numpy.sqrt(2.0 * this * thisEzZmaxs), - **kwargs + **kwargs, )[2] jz = numpy.reshape(jz, (nR, nEz)) jzEzzmax[0:nR] = jz[:, nEz - 1] @@ -118,7 +118,7 @@ def __init__( 0.0, numpy.sqrt(2.0 * this[x] * thisEzZmaxs[x]), _justjz=True, - **kwargs + **kwargs, )[2] ), range(nR * nEz), @@ -137,7 +137,7 @@ def __init__( 0.0, numpy.sqrt(2.0 * y[jj] * self._EzZmaxs[ii]), _justjz=True, - **kwargs + **kwargs, )[2][0] if jj == nEz - 1: jzEzzmax[ii] = jz[ii, jj] @@ -204,7 +204,7 @@ def __init__( thisLzs / thisRL, numpy.zeros(len(thisRL)), numpy.zeros(len(thisRL)), - **kwargs + **kwargs, )[0] jr[:, 0:-1] = numpy.reshape(mjr, (nLz, nEr - 1)) jrERRa[0:nLz] = jr[:, 0] @@ -227,7 +227,7 @@ def __init__( 0.0, 0.0, _justjr=True, - **kwargs + **kwargs, )[0] ), range((nEr - 1) * nLz), @@ -256,7 +256,7 @@ def __init__( 0.0, 0.0, _justjr=True, - **kwargs + **kwargs, )[0][0] except UnboundError: # pragma: no cover raise @@ -338,7 +338,7 @@ def _evaluate(self, *args, **kwargs): numpy.zeros(numpy.sum(indx)), numpy.sqrt(2.0 * Ez[indx]), _justjz=True, - **kwargs + **kwargs, )[2] else: if ( @@ -360,7 +360,7 @@ def _evaluate(self, *args, **kwargs): 0.0, numpy.sqrt(2.0 * Ez), _justjz=True, - **kwargs + **kwargs, )[2] else: jz = ( @@ -405,7 +405,7 @@ def _evaluate(self, *args, **kwargs): numpy.zeros(len(thisRL)), numpy.zeros(len(thisRL)), _justjr=True, - **kwargs + **kwargs, )[0] else: if (ER - thisERRa) / (thisERRL - thisERRa) > 1.0 and ( @@ -445,7 +445,7 @@ def _evaluate(self, *args, **kwargs): 0.0, 0.0, _justjr=True, - **kwargs + **kwargs, )[0] else: jr = ( @@ -498,7 +498,7 @@ def Jz(self, *args, **kwargs): 0.0, numpy.sqrt(2.0 * Ez), _justjz=True, - **kwargs + **kwargs, )[2] else: jz = ( diff --git a/galpy/actionAngle/actionAngleIsochroneApprox.py b/galpy/actionAngle/actionAngleIsochroneApprox.py index 64003e349..9d60702ff 100644 --- a/galpy/actionAngle/actionAngleIsochroneApprox.py +++ b/galpy/actionAngle/actionAngleIsochroneApprox.py @@ -547,7 +547,7 @@ def plot(self, *args, **kwargs): vmax=2.0 * numpy.pi, crange=[0.0, 2.0 * numpy.pi], colorbar=True, - **kwargs + **kwargs, ) elif type == "lz": if downsample: @@ -572,7 +572,7 @@ def plot(self, *args, **kwargs): vmax=2.0 * numpy.pi, crange=[0.0, 2.0 * numpy.pi], colorbar=True, - **kwargs + **kwargs, ) elif type == "jz": if downsample: @@ -597,7 +597,7 @@ def plot(self, *args, **kwargs): vmax=2.0 * numpy.pi, crange=[0.0, 2.0 * numpy.pi], colorbar=True, - **kwargs + **kwargs, ) else: if deperiod: @@ -658,7 +658,7 @@ def plot(self, *args, **kwargs): vmax=vmax, crange=crange, colorbar=True, - **kwargs + **kwargs, ) elif type == "araphi": if downsample: @@ -685,7 +685,7 @@ def plot(self, *args, **kwargs): vmax=vmax, crange=crange, colorbar=True, - **kwargs + **kwargs, ) elif type == "azaphi": if downsample: @@ -712,7 +712,7 @@ def plot(self, *args, **kwargs): vmax=vmax, crange=crange, colorbar=True, - **kwargs + **kwargs, ) return None diff --git a/galpy/actionAngle/actionAngleSpherical.py b/galpy/actionAngle/actionAngleSpherical.py index 08aeb0c63..1ceeb7cb6 100644 --- a/galpy/actionAngle/actionAngleSpherical.py +++ b/galpy/actionAngle/actionAngleSpherical.py @@ -342,7 +342,7 @@ def _actionsFreqsAngles(self, *args, **kwargs): E[ii], L[ii], fixed_quad, - **kwargs + **kwargs, ) ) # Angles @@ -357,7 +357,7 @@ def _actionsFreqsAngles(self, *args, **kwargs): L[ii], vr[ii], fixed_quad, - **kwargs + **kwargs, ) ) az.append( @@ -377,7 +377,7 @@ def _actionsFreqsAngles(self, *args, **kwargs): vtheta[ii], phi[ii], fixed_quad, - **kwargs + **kwargs, ) ) Op = numpy.array(Op) @@ -535,7 +535,7 @@ def _calc_jr(self, rperi, rap, E, L, fixed_quad, **kwargs): rap, args=(E, L, self._2dpot), n=10, - **kwargs + **kwargs, )[0] / numpy.pi ) @@ -547,7 +547,7 @@ def _calc_jr(self, rperi, rap, E, L, fixed_quad, **kwargs): rperi, rap, args=(E, L, self._2dpot), - **kwargs + **kwargs, ) ) )[0] / numpy.pi @@ -561,7 +561,7 @@ def _calc_or(self, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): 0.0, numpy.sqrt(Rmean - rperi), args=(E, L, self._2dpot, rperi), - **kwargs + **kwargs, ) )[0] elif Rmean > rperi and fixed_quad: @@ -571,7 +571,7 @@ def _calc_or(self, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): numpy.sqrt(Rmean - rperi), args=(E, L, self._2dpot, rperi), n=10, - **kwargs + **kwargs, )[0] if Rmean < rap and not fixed_quad: Tr += numpy.array( @@ -580,7 +580,7 @@ def _calc_or(self, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): 0.0, numpy.sqrt(rap - Rmean), args=(E, L, self._2dpot, rap), - **kwargs + **kwargs, ) )[0] elif Rmean < rap and fixed_quad: @@ -590,7 +590,7 @@ def _calc_or(self, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): numpy.sqrt(rap - Rmean), args=(E, L, self._2dpot, rap), n=10, - **kwargs + **kwargs, )[0] Tr = 2.0 * Tr return 2.0 * numpy.pi / Tr @@ -605,7 +605,7 @@ def _calc_op(self, Or, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): 0.0, numpy.sqrt(Rmean - rperi), args=(E, L, self._2dpot, rperi), - **kwargs + **kwargs, ) )[0] elif Rmean > rperi and fixed_quad: @@ -615,7 +615,7 @@ def _calc_op(self, Or, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): numpy.sqrt(Rmean - rperi), args=(E, L, self._2dpot, rperi), n=10, - **kwargs + **kwargs, )[0] if Rmean < rap and not fixed_quad: I += numpy.array( @@ -624,7 +624,7 @@ def _calc_op(self, Or, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): 0.0, numpy.sqrt(rap - Rmean), args=(E, L, self._2dpot, rap), - **kwargs + **kwargs, ) )[0] elif Rmean < rap and fixed_quad: @@ -634,7 +634,7 @@ def _calc_op(self, Or, Rmean, rperi, rap, E, L, fixed_quad, **kwargs): numpy.sqrt(rap - Rmean), args=(E, L, self._2dpot, rap), n=10, - **kwargs + **kwargs, )[0] I *= 2 * L return I * Or / 2.0 / numpy.pi @@ -663,7 +663,7 @@ def _calc_angler(self, Or, r, Rmean, rperi, rap, E, L, vr, fixed_quad, **kwargs) 0.0, numpy.sqrt(r - rperi), args=(E, L, self._2dpot, rperi), - **kwargs + **kwargs, )[0] ) elif r > rperi and fixed_quad: @@ -675,7 +675,7 @@ def _calc_angler(self, Or, r, Rmean, rperi, rap, E, L, vr, fixed_quad, **kwargs) numpy.sqrt(r - rperi), args=(E, L, self._2dpot, rperi), n=10, - **kwargs + **kwargs, )[0] ) else: @@ -691,7 +691,7 @@ def _calc_angler(self, Or, r, Rmean, rperi, rap, E, L, vr, fixed_quad, **kwargs) 0.0, numpy.sqrt(rap - r), args=(E, L, self._2dpot, rap), - **kwargs + **kwargs, )[0] ) elif r < rap and fixed_quad: @@ -703,7 +703,7 @@ def _calc_angler(self, Or, r, Rmean, rperi, rap, E, L, vr, fixed_quad, **kwargs) numpy.sqrt(rap - r), args=(E, L, self._2dpot, rap), n=10, - **kwargs + **kwargs, )[0] ) else: @@ -731,7 +731,7 @@ def _calc_anglez( vtheta, phi, fixed_quad, - **kwargs + **kwargs, ): # First calculate psi i = numpy.arccos(Lz / L) @@ -755,7 +755,7 @@ def _calc_anglez( 0.0, numpy.sqrt(r - rperi), args=(E, L, self._2dpot, rperi), - **kwargs + **kwargs, )[0] ) elif fixed_quad: @@ -767,7 +767,7 @@ def _calc_anglez( numpy.sqrt(r - rperi), args=(E, L, self._2dpot, rperi), n=10, - **kwargs + **kwargs, )[0] ) if vr < 0.0: @@ -781,7 +781,7 @@ def _calc_anglez( 0.0, numpy.sqrt(rap - r), args=(E, L, self._2dpot, rap), - **kwargs + **kwargs, )[0] ) elif fixed_quad: @@ -793,7 +793,7 @@ def _calc_anglez( numpy.sqrt(rap - r), args=(E, L, self._2dpot, rap), n=10, - **kwargs + **kwargs, )[0] ) if vr < 0.0: diff --git a/galpy/actionAngle/actionAngleStaeckel.py b/galpy/actionAngle/actionAngleStaeckel.py index a91d1689f..b0ad34bfe 100644 --- a/galpy/actionAngle/actionAngleStaeckel.py +++ b/galpy/actionAngle/actionAngleStaeckel.py @@ -768,7 +768,7 @@ def JR(self, **kwargs): self._pot, ), n=order, - **kwargs + **kwargs, )[0] ) else: @@ -793,7 +793,7 @@ def JR(self, **kwargs): self._potu0v0, self._pot, ), - **kwargs + **kwargs, )[0] ) return self._JR @@ -847,7 +847,7 @@ def Jz(self, **kwargs): self._pot, ), n=order, - **kwargs + **kwargs, )[0] ) else: @@ -872,7 +872,7 @@ def Jz(self, **kwargs): self._potupi2, self._pot, ), - **kwargs + **kwargs, )[0] ) return self._JZ diff --git a/galpy/actionAngle/actionAngleStaeckelGrid.py b/galpy/actionAngle/actionAngleStaeckelGrid.py index 8a36515c5..3c36d70d0 100644 --- a/galpy/actionAngle/actionAngleStaeckelGrid.py +++ b/galpy/actionAngle/actionAngleStaeckelGrid.py @@ -37,7 +37,7 @@ def __init__( nLz=30, numcores=1, interpecc=False, - **kwargs + **kwargs, ): """ Initialize an actionAngleStaeckelGrid object @@ -507,7 +507,7 @@ def _evaluate(self, *args, **kwargs): numpy.array([vT]), numpy.array([z]), numpy.array([vz]), - **kwargs + **kwargs, ) return (jr[0], Lz[0], jz[0]) jr[jr < 0.0] = 0.0 @@ -746,7 +746,7 @@ def _EccZmaxRperiRap(self, *args, **kwargs): numpy.array([vT]), numpy.array([z]), numpy.array([vz]), - **kwargs + **kwargs, ) return (ecc[0], zmax[0], rperi[0], rap[0]) ecc[ecc < 0.0] = 0.0 diff --git a/galpy/df/constantbetadf.py b/galpy/df/constantbetadf.py index 5d724ebbd..2c72016ba 100644 --- a/galpy/df/constantbetadf.py +++ b/galpy/df/constantbetadf.py @@ -123,9 +123,9 @@ def _p_v_at_r(self, v, r): _evaluatePotentials(self._pot, r, 0) + 0.5 * v**2.0 ) * v ** (2.0 - 2.0 * self._beta) else: - return self.fE( - _evaluatePotentials(self._pot, r, 0) + 0.5 * v**2.0 - ) * v ** (2.0 - 2.0 * self._beta) + return self.fE(_evaluatePotentials(self._pot, r, 0) + 0.5 * v**2.0) * v ** ( + 2.0 - 2.0 * self._beta + ) def _vmomentdensity(self, r, n, m): if m % 2 == 1 or n % 2 == 1: diff --git a/galpy/df/diskdf.py b/galpy/df/diskdf.py index e4ecd151e..bcccf7bf4 100644 --- a/galpy/df/diskdf.py +++ b/galpy/df/diskdf.py @@ -24,9 +24,7 @@ import numpy import scipy -numpylog = ( - numpy.lib.scimath.log -) # somehow, this code produces log(negative), which scipy (now numpy.lib.scimath.log) implements as log(|negative|) + i pi while numpy gives NaN and we want the scipy behavior; not sure where the log(negative) comes from though! I think it's for sigma=0 DFs (this test fails with numpy.log) where the DF eval has a log(~zero) that can be slightly negative because of numerical precision issues +numpylog = numpy.lib.scimath.log # somehow, this code produces log(negative), which scipy (now numpy.lib.scimath.log) implements as log(|negative|) + i pi while numpy gives NaN and we want the scipy behavior; not sure where the log(negative) comes from though! I think it's for sigma=0 DFs (this test fails with numpy.log) where the DF eval has a log(~zero) that can be slightly negative because of numerical precision issues from scipy import integrate, interpolate, optimize, stats from ..actionAngle import actionAngleAdiabatic @@ -67,7 +65,7 @@ def __init__( beta=0.0, ro=None, vo=None, - **kwargs + **kwargs, ): """ Initialize a DF @@ -129,7 +127,7 @@ def __init__( dftype=self.__class__, surfaceSigmaProfile=self._surfaceSigmaProfile, beta=beta, - **kwargs + **kwargs, ) else: self._correct = False @@ -282,7 +280,7 @@ def _call_marginalizevperp(self, o, **kwargs): vcirc, sigmaR1 / self._gamma, ), - **kwargs + **kwargs, )[0] / numpy.fabs(cosalphalos) * sigmaR1 @@ -305,7 +303,7 @@ def _call_marginalizevperp(self, o, **kwargs): vcirc, sigmaR1, ), - **kwargs + **kwargs, )[0] / numpy.fabs(sinalphalos) * sigmaR1 @@ -363,7 +361,7 @@ def _call_marginalizevlos(self, o, **kwargs): vcirc, sigmaR1 / self._gamma, ), - **kwargs + **kwargs, )[0] / numpy.fabs(cosalphaperp) * sigmaR1 @@ -387,7 +385,7 @@ def _call_marginalizevlos(self, o, **kwargs): vcirc, sigmaR1, ), - **kwargs + **kwargs, )[0] / numpy.fabs(sinalphaperp) * sigmaR1 @@ -1535,9 +1533,9 @@ def kurtosisvT(self, R, romberg=False, nsigma=None, phi=0.0): self._vmomentsurfacemass(R, 0, 4, romberg=romberg, nsigma=nsigma) / surfmass ) s2 = vt2 - vt**2.0 - return ( - vt4 - 4.0 * vt * vt3 + 6.0 * vt**2.0 * vt2 - 3.0 * vt**4.0 - ) * s2 ** (-2.0) - 3.0 + return (vt4 - 4.0 * vt * vt3 + 6.0 * vt**2.0 * vt2 - 3.0 * vt**4.0) * s2 ** ( + -2.0 + ) - 3.0 @potential_physical_input def kurtosisvR(self, R, romberg=False, nsigma=None, phi=0.0): @@ -1580,9 +1578,9 @@ def kurtosisvR(self, R, romberg=False, nsigma=None, phi=0.0): self._vmomentsurfacemass(R, 4, 0, romberg=romberg, nsigma=nsigma) / surfmass ) s2 = vr2 - vr**2.0 - return ( - vr4 - 4.0 * vr * vr3 + 6.0 * vr**2.0 * vr2 - 3.0 * vr**4.0 - ) * s2 ** (-2.0) - 3.0 + return (vr4 - 4.0 * vr * vr3 + 6.0 * vr**2.0 * vr2 - 3.0 * vr**4.0) * s2 ** ( + -2.0 + ) - 3.0 def _ELtowRRapRperi(self, E, L): """ @@ -1794,7 +1792,7 @@ def __init__( profileParams=(1.0 / 3.0, 1.0, 0.2), correct=False, beta=0.0, - **kwargs + **kwargs, ): """ Initialize a Dehnen 'new' DF. @@ -1832,7 +1830,7 @@ def __init__( correct=correct, dftype="dehnen", beta=beta, - **kwargs + **kwargs, ) def eval(self, E, L, logSigmaR=0.0, logsigmaR2=0.0): @@ -1938,7 +1936,7 @@ def sample( targetSurfmass=True, targetSigma2=True, maxd=None, - **kwargs + **kwargs, ): """ Sample n*nphi points from this DF. @@ -2245,7 +2243,7 @@ def __init__( profileParams=(1.0 / 3.0, 1.0, 0.2), correct=False, beta=0.0, - **kwargs + **kwargs, ): """ Initialize a Shu DF. @@ -2283,7 +2281,7 @@ def __init__( correct=correct, dftype="shu", beta=beta, - **kwargs + **kwargs, ) def eval(self, E, L, logSigmaR=0.0, logsigmaR2=0.0): @@ -2355,7 +2353,7 @@ def sample( maxd=None, targetSurfmass=True, targetSigma2=True, - **kwargs + **kwargs, ): """ Sample n*nphi points from this DF. @@ -2593,7 +2591,7 @@ def __init__( profileParams=(1.0 / 3.0, 1.0, 0.2), correct=False, beta=0.0, - **kwargs + **kwargs, ): """ Initialize a Schwarzschild DF. @@ -2634,7 +2632,7 @@ def __init__( correct=correct, dftype="schwarzschild", beta=beta, - **kwargs + **kwargs, ) diff --git a/galpy/df/evolveddiskdf.py b/galpy/df/evolveddiskdf.py index 49214c0ef..a4eaeb904 100644 --- a/galpy/df/evolveddiskdf.py +++ b/galpy/df/evolveddiskdf.py @@ -3003,7 +3003,7 @@ def _call_marginalizevperp(self, o, integrate_method="dopr54_c", **kwargs): sigmaR1, phi, ), - **kwargs + **kwargs, )[0] / numpy.fabs(cosalphalos) * sigmaR1 @@ -3027,7 +3027,7 @@ def _call_marginalizevperp(self, o, integrate_method="dopr54_c", **kwargs): sigmaR1, phi, ), - **kwargs + **kwargs, )[0] / numpy.fabs(sinalphalos) * sigmaR1 @@ -3077,7 +3077,7 @@ def _call_marginalizevlos(self, o, integrate_method="dopr54_c", **kwargs): sigmaR1, phi, ), - **kwargs + **kwargs, )[0] / numpy.fabs(cosalphaperp) * sigmaR1 @@ -3102,7 +3102,7 @@ def _call_marginalizevlos(self, o, integrate_method="dopr54_c", **kwargs): sigmaR1, phi, ), - **kwargs + **kwargs, )[0] / numpy.fabs(sinalphaperp) * sigmaR1 diff --git a/galpy/df/quasiisothermaldf.py b/galpy/df/quasiisothermaldf.py index 7b8451917..8afb22002 100644 --- a/galpy/df/quasiisothermaldf.py +++ b/galpy/df/quasiisothermaldf.py @@ -604,7 +604,7 @@ def _vmomentdensity( _Omega=None, _sigmaR1=None, _sigmaz1=None, - **kwargs + **kwargs, ): """Non-physical version of vmomentdensity, otherwise the same""" if isinstance(R, numpy.ndarray): @@ -621,7 +621,7 @@ def _vmomentdensity( nmc=nmc, gl=gl, ngl=ngl, - **kwargs + **kwargs, ) for r, zz in zip(R, z) ] @@ -904,7 +904,7 @@ def _vmomentdensity( lambda x, y: 0.0, lambda x, y: nsigma, (R, z, self, sigmaR1, gamma, sigmaz1, n, m, o), - **kwargs + **kwargs, )[0] * sigmaR1 ** (2.0 + n + m) * gamma ** (1.0 + m) @@ -983,7 +983,7 @@ def _jmomentdensity( _vrs=None, _vts=None, _vzs=None, - **kwargs + **kwargs, ): """Non-physical version of jmomentdensity, otherwise the same""" if nsigma == None: @@ -1058,7 +1058,7 @@ def _jmomentdensity( lambda x, y: 0.0, lambda x, y: nsigma, (R, z, self, sigmaR1, gamma, sigmaz1, n, m, o), - **kwargs + **kwargs, )[0] * sigmaR1**2.0 * gamma @@ -1154,7 +1154,7 @@ def sigmaR2( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1170,7 +1170,7 @@ def sigmaR2( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1239,7 +1239,7 @@ def sigmaRz( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1255,7 +1255,7 @@ def sigmaRz( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1322,7 +1322,7 @@ def tilt( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) tsigmar2 = ( self._vmomentdensity( @@ -1338,7 +1338,7 @@ def tilt( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1356,7 +1356,7 @@ def tilt( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1374,7 +1374,7 @@ def tilt( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1453,7 +1453,7 @@ def sigmaz2( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1469,7 +1469,7 @@ def sigmaz2( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1538,7 +1538,7 @@ def meanvT( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1554,7 +1554,7 @@ def meanvT( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1623,7 +1623,7 @@ def meanvR( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1639,7 +1639,7 @@ def meanvR( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1707,7 +1707,7 @@ def meanvz( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._vmomentdensity( @@ -1723,7 +1723,7 @@ def meanvz( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1792,7 +1792,7 @@ def sigmaT2( mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) mvt = ( self._vmomentdensity( @@ -1808,7 +1808,7 @@ def sigmaT2( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1826,7 +1826,7 @@ def sigmaT2( _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass - mvt**2.0 @@ -1909,7 +1909,7 @@ def meanjr(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._jmomentdensity( @@ -1925,7 +1925,7 @@ def meanjr(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -1979,7 +1979,7 @@ def meanlz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._jmomentdensity( @@ -1995,7 +1995,7 @@ def meanlz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -2048,7 +2048,7 @@ def meanjz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): mc=mc, nmc=nmc, _returnmc=True, - **kwargs + **kwargs, ) return ( self._jmomentdensity( @@ -2064,7 +2064,7 @@ def meanjz(self, R, z, nsigma=None, mc=True, nmc=10000, **kwargs): _vrs=vrs, _vts=vts, _vzs=vzs, - **kwargs + **kwargs, ) / surfmass ) @@ -2169,7 +2169,7 @@ def sampleV_interpolate( R_min=None, R_max=None, z_max=None, - **kwargs + **kwargs, ): """ Sample radial, azimuthal, and vertical velocity at R,z using interpolation. diff --git a/galpy/df/streamdf.py b/galpy/df/streamdf.py index 704f741ff..0564f31bb 100644 --- a/galpy/df/streamdf.py +++ b/galpy/df/streamdf.py @@ -638,7 +638,7 @@ def plotTrack( *args, xlabel=_labelDict[d1.lower()], ylabel=_labelDict[d2.lower()], - **kwargs + **kwargs, ) if spread: addx, addy = self._parse_track_spread( @@ -718,7 +718,7 @@ def plotProgenitor(self, d1="x", d2="z", *args, **kwargs): *args, xlabel=_labelDict[d1.lower()], ylabel=_labelDict[d2.lower()], - **kwargs + **kwargs, ) return None @@ -1042,7 +1042,7 @@ def plotCompareTrackAAModel(self, **kwargs): xlabel=xlabel, ylabel=ylabel, yrange=yrange, - **kwargs + **kwargs, ) plot.plot(track_adiff, track_operp, "ko", overplot=True, **kwargs) return None @@ -3823,9 +3823,7 @@ def _determine_stream_spread_single( ) # angles sigangle2 = ( - sigAngle(thetasTrack) ** 2.0 - if hasattr(sigAngle, "__call__") - else sigAngle**2.0 + sigAngle(thetasTrack) ** 2.0 if hasattr(sigAngle, "__call__") else sigAngle**2.0 ) tsigadiag = numpy.ones(3) * sigangle2 tsigadiag[numpy.argmax(tsigOdiag)] = 1.0 diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 8eb1f91bc..d463abd40 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -1073,9 +1073,7 @@ def _plot(*args, **kwargs): Gravitational field to use. Default is the gravitational field used to integrate the orbit. normed : bool, optional if set, plot {quant}(t)/{quant}(0) rather than {quant}(t) - """.format( - quant=name.split("plot")[1] - ) + """.format(quant=name.split("plot")[1]) else: Estring = "" _plot.__doc__ = """Plot {quant}(t) along the orbit. @@ -1104,9 +1102,7 @@ def _plot(*args, **kwargs): ----- - 2019-04-13 - Written - Bovy (UofT) - """.format( - quant=name.split("plot")[1], Estring=Estring - ) + """.format(quant=name.split("plot")[1], Estring=Estring) return _plot else: raise AttributeError( @@ -6235,14 +6231,18 @@ def animate(self, **kwargs): # pragma: no cover margin: {{t: 20}}, hovermode: 'closest', showlegend: false, -""".format( - xlabel=xlabels[0], ylabel=ylabels[0], xmin=xmin[0], xmax=xmax[0] - ) +""".format(xlabel=xlabels[0], ylabel=ylabels[0], xmin=xmin[0], xmax=xmax[0]) hovertemplate = ( - lambda name, xlabel, ylabel, tlabel: f"""'{name}' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" + lambda name, + xlabel, + ylabel, + tlabel: f"""'{name}' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" ) hovertemplate_current = ( - lambda name, xlabel, ylabel, tlabel: f"""'{name} (Current location)' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" + lambda name, + xlabel, + ylabel, + tlabel: f"""'{name} (Current location)' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" ) for ii in range(1, nplots): layout += """ xaxis{idx}: {{ @@ -6458,9 +6458,7 @@ def animate(self, **kwargs): # pragma: no cover else: # else for "if there is a 2nd panel" setup_trace2 = """ let traces= [{traces_cumul}]; -""".format( - traces_cumul=traces_cumul - ) +""".format(traces_cumul=traces_cumul) if len(d1s) > 2: setup_trace3 = """ let trace{trace_num_1}= {{ @@ -6558,15 +6556,11 @@ def animate(self, **kwargs): # pragma: no cover traces_cumul += f""",trace{str(4*self.size+2*ii+1)},trace{str(4*self.size+2*ii+2)}""" setup_trace3 += """ let traces= [{traces_cumul}]; - """.format( - traces_cumul=traces_cumul - ) + """.format(traces_cumul=traces_cumul) elif len(d1s) > 1: # elif for "if there is a 3rd panel setup_trace3 = """ let traces= [{traces_cumul}]; -""".format( - traces_cumul=traces_cumul - ) +""".format(traces_cumul=traces_cumul) else: # else for "if there is a 3rd or 2nd panel" (don't think we can get here!) setup_trace3 = "" return HTML( @@ -7007,10 +7001,18 @@ def animate3d(self, mw_plane_bg=False, **kwargs): # pragma: no cover ] ) hovertemplate = ( - lambda name, xlabel, ylabel, zlabel, tlabel: f"""'{name}' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{zlabel}: %{{z:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" + lambda name, + xlabel, + ylabel, + zlabel, + tlabel: f"""'{name}' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{zlabel}: %{{z:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" ) hovertemplate_current = ( - lambda name, xlabel, ylabel, zlabel, tlabel: f"""'{name} (Current location)' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{zlabel}: %{{z:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" + lambda name, + xlabel, + ylabel, + zlabel, + tlabel: f"""'{name} (Current location)' + '
{xlabel}: %{{x:.2f}}' + '
{ylabel}: %{{y:.2f}}' + '
{zlabel}: %{{z:.2f}}' + '
{tlabel}: %{{customdata:.2f}}'""" ) layout = """{{ scene:{{ diff --git a/galpy/orbit/integratePlanarOrbit.py b/galpy/orbit/integratePlanarOrbit.py index 0108192d3..c9b79e4a1 100644 --- a/galpy/orbit/integratePlanarOrbit.py +++ b/galpy/orbit/integratePlanarOrbit.py @@ -560,9 +560,7 @@ def _parse_pot(pot): or isinstance(p, planarPotentialFromRZPotential) ) and isinstance(p._Pot, potential.RotateAndTiltWrapperPotential) - ) or isinstance( - p, potential.RotateAndTiltWrapperPotential - ): # pragma: no cover + ) or isinstance(p, potential.RotateAndTiltWrapperPotential): # pragma: no cover raise NotImplementedError( "Planar orbit integration in C for RotateAndTiltWrapperPotential not implemented; please integrate an orbit with (z,vz) = (0,0) instead" ) @@ -656,7 +654,8 @@ def _prep_tfuncs(pot_tfuncs): pot_tfuncs = None # NULL else: func_ctype = ctypes.CFUNCTYPE( - ctypes.c_double, ctypes.c_double # Return type + ctypes.c_double, + ctypes.c_double, # Return type ) # time try: # using numba if not _NUMBA_LOADED: @@ -1436,7 +1435,7 @@ def integrate_for_map(vxvv, psi, t0): this_vxvv, t=psi + init_psi, args=(pot,), - **extra_kwargs + **extra_kwargs, ) # go back to the cylindrical frame x = intOut[:, 2] * numpy.sin(psi + init_psi) @@ -1461,7 +1460,7 @@ def integrate_for_map(vxvv, psi, t0): this_vxvv, t=psi + init_psi, args=(pot,), - **extra_kwargs + **extra_kwargs, ) # go back to the cylindrical frame x = intOut[:, 0] diff --git a/galpy/potential/BurkertPotential.py b/galpy/potential/BurkertPotential.py index 87ad5627f..da375a39f 100644 --- a/galpy/potential/BurkertPotential.py +++ b/galpy/potential/BurkertPotential.py @@ -60,7 +60,7 @@ def _revaluate(self, r, t=0.0): """Potential as a function of r and time""" x = r / self.a return ( - -self.a**2.0 + -(self.a**2.0) * numpy.pi * ( -numpy.pi / x diff --git a/galpy/potential/CorotatingRotationWrapperPotential.py b/galpy/potential/CorotatingRotationWrapperPotential.py index 7176de7bc..83ecd4751 100644 --- a/galpy/potential/CorotatingRotationWrapperPotential.py +++ b/galpy/potential/CorotatingRotationWrapperPotential.py @@ -134,12 +134,6 @@ def _Rphideriv(self, *args, **kwargs): self._pot, *args, **kwargs ) - self._wrap_pot_func("_phi2deriv")( self._pot, *args, **kwargs - ) * self._vpo * ( - self._beta - 1.0 - ) * args[ - 0 - ] ** ( - self._beta - 2.0 - ) * ( + ) * self._vpo * (self._beta - 1.0) * args[0] ** (self._beta - 2.0) * ( kwargs.get("t", 0.0) - self._to ) diff --git a/galpy/potential/FerrersPotential.py b/galpy/potential/FerrersPotential.py index c323c1cbc..03dc0989e 100644 --- a/galpy/potential/FerrersPotential.py +++ b/galpy/potential/FerrersPotential.py @@ -373,9 +373,9 @@ def _potInt(x, y, z, a2, b2, c2, n): def integrand(tau): return _FracInt(x, y, z, a2, b2, c2, tau, n + 1) - return integrate.quad( - integrand, lowerlim(x**2, y**2, z**2, a2, b2, c2), numpy.inf - )[0] + return integrate.quad(integrand, lowerlim(x**2, y**2, z**2, a2, b2, c2), numpy.inf)[ + 0 + ] def _forceInt(x, y, z, a2, b2, c2, n, i): @@ -427,9 +427,9 @@ def integrand(tau): tau + coef2 ) ** 2 + _FracInt(x, y, z, a2, b2, c2, tau, n) * (-2.0 / (tau + coef2)) - return integrate.quad( - integrand, lowerlim(x**2, y**2, z**2, a2, b2, c2), numpy.inf - )[0] + return integrate.quad(integrand, lowerlim(x**2, y**2, z**2, a2, b2, c2), numpy.inf)[ + 0 + ] def _FracInt(x, y, z, a, b, c, tau, n): diff --git a/galpy/potential/IsochronePotential.py b/galpy/potential/IsochronePotential.py index 371a2da73..031317bee 100644 --- a/galpy/potential/IsochronePotential.py +++ b/galpy/potential/IsochronePotential.py @@ -77,7 +77,7 @@ def _R2deriv(self, R, z, phi=0.0, t=0.0): rb = numpy.sqrt(r2 + self.b2) return ( -( - -self.b**3.0 + -(self.b**3.0) - self.b * z**2.0 + (2.0 * R**2.0 - z**2.0 - self.b**2.0) * rb ) @@ -90,7 +90,7 @@ def _z2deriv(self, R, z, phi=0.0, t=0.0): rb = numpy.sqrt(r2 + self.b2) return ( -( - -self.b**3.0 + -(self.b**3.0) - self.b * R**2.0 - (R**2.0 - 2.0 * z**2.0 + self.b**2.0) * rb ) diff --git a/galpy/potential/KuzminDiskPotential.py b/galpy/potential/KuzminDiskPotential.py index 877a2f8d8..08ee4e004 100644 --- a/galpy/potential/KuzminDiskPotential.py +++ b/galpy/potential/KuzminDiskPotential.py @@ -55,10 +55,10 @@ def __init__(self, amp=1.0, a=1.0, normalize=False, ro=None, vo=None): return None def _evaluate(self, R, z, phi=0.0, t=0.0): - return -self._denom(R, z) ** -0.5 + return -(self._denom(R, z) ** -0.5) def _Rforce(self, R, z, phi=0.0, t=0.0): - return -self._denom(R, z) ** -1.5 * R + return -(self._denom(R, z) ** -1.5) * R def _zforce(self, R, z, phi=0.0, t=0.0): return -numpy.sign(z) * self._denom(R, z) ** -1.5 * (self._a + numpy.fabs(z)) diff --git a/galpy/potential/KuzminLikeWrapperPotential.py b/galpy/potential/KuzminLikeWrapperPotential.py index e08558730..8b6a60312 100644 --- a/galpy/potential/KuzminLikeWrapperPotential.py +++ b/galpy/potential/KuzminLikeWrapperPotential.py @@ -133,24 +133,18 @@ def _R2deriv(self, R, z, phi=0.0, t=0.0): self._pot, self._xi(R, z), 0.0, phi=phi, t=t ) * self._dxidR(R, z) ** 2.0 - _evaluateRforces( self._pot, self._xi(R, z), 0.0, phi=phi, t=t - ) * self._d2xidR2( - R, z - ) + ) * self._d2xidR2(R, z) def _z2deriv(self, R, z, phi=0.0, t=0.0): return evaluateR2derivs( self._pot, self._xi(R, z), 0.0, phi=phi, t=t ) * self._dxidz(R, z) ** 2.0 - _evaluateRforces( self._pot, self._xi(R, z), 0.0, phi=phi, t=t - ) * self._d2xidz2( - R, z - ) + ) * self._d2xidz2(R, z) def _Rzderiv(self, R, z, phi=0.0, t=0.0): return evaluateR2derivs( self._pot, self._xi(R, z), 0.0, phi=phi, t=t ) * self._dxidR(R, z) * self._dxidz(R, z) - _evaluateRforces( self._pot, self._xi(R, z), 0.0, phi=phi, t=t - ) * self._d2xidRdz( - R, z - ) + ) * self._d2xidRdz(R, z) diff --git a/galpy/potential/MiyamotoNagaiPotential.py b/galpy/potential/MiyamotoNagaiPotential.py index 871358299..db8aa7b7f 100644 --- a/galpy/potential/MiyamotoNagaiPotential.py +++ b/galpy/potential/MiyamotoNagaiPotential.py @@ -78,9 +78,9 @@ def _zforce(self, R, z, phi=0.0, t=0.0): sqrtbz = numpy.sqrt(self._b2 + z**2.0) asqrtbz = self._a + sqrtbz if isinstance(R, float) and sqrtbz == asqrtbz: - return -z / ( - R**2.0 + (self._a + numpy.sqrt(z**2.0 + self._b2)) ** 2.0 - ) ** (3.0 / 2.0) + return -z / (R**2.0 + (self._a + numpy.sqrt(z**2.0 + self._b2)) ** 2.0) ** ( + 3.0 / 2.0 + ) else: return ( -z diff --git a/galpy/potential/PowerSphericalPotential.py b/galpy/potential/PowerSphericalPotential.py index 8435aaeac..cc90d805b 100644 --- a/galpy/potential/PowerSphericalPotential.py +++ b/galpy/potential/PowerSphericalPotential.py @@ -191,9 +191,9 @@ def _R2deriv(self, R, z, phi=0.0, t=0.0): ----- - 2011-10-09 - Written - Bovy (NYU) """ - return 1.0 / (R**2.0 + z**2.0) ** ( - self.alpha / 2.0 - ) - self.alpha * R**2.0 / (R**2.0 + z**2.0) ** (self.alpha / 2.0 + 1.0) + return 1.0 / (R**2.0 + z**2.0) ** (self.alpha / 2.0) - self.alpha * R**2.0 / ( + R**2.0 + z**2.0 + ) ** (self.alpha / 2.0 + 1.0) def _z2deriv(self, R, z, phi=0.0, t=0.0): """ diff --git a/galpy/potential/RingPotential.py b/galpy/potential/RingPotential.py index b16d1bafd..121dde8a9 100644 --- a/galpy/potential/RingPotential.py +++ b/galpy/potential/RingPotential.py @@ -95,9 +95,9 @@ def _R2deriv(self, R, z, phi=0.0, t=0.0): Raz = numpy.sqrt(Raz2) m = 4.0 * R * self.a / Raz2 R2ma2mz2o4aR1m = (R**2 - self.a2 - z**2) / 4.0 / self.a / R / (1.0 - m) - return ( - 2 * R**2 + self.a2 + 3 * R * self.a + z**2 - ) / R / Raz2 * self._Rforce(R, z) + 2.0 * self.a / R / Raz * ( + return (2 * R**2 + self.a2 + 3 * R * self.a + z**2) / R / Raz2 * self._Rforce( + R, z + ) + 2.0 * self.a / R / Raz * ( m * (R**2 + self.a2 + z**2) / 4.0 diff --git a/galpy/potential/SpiralArmsPotential.py b/galpy/potential/SpiralArmsPotential.py index 85efb8c4e..efb3fb5ab 100644 --- a/galpy/potential/SpiralArmsPotential.py +++ b/galpy/potential/SpiralArmsPotential.py @@ -95,12 +95,8 @@ def __init__( Rs = conversion.parse_length(Rs, ro=self._ro) H = conversion.parse_length(H, ro=self._ro) omega = conversion.parse_frequency(omega, ro=self._ro, vo=self._vo) - self._N = ( - -N - ) # trick to flip to left handed coordinate system; flips sign for phi and phi_ref, but also alpha. - self._alpha = ( - -alpha - ) # we don't want sign for alpha to change, so flip alpha. (see eqn. 3 in the paper) + self._N = -N # trick to flip to left handed coordinate system; flips sign for phi and phi_ref, but also alpha. + self._alpha = -alpha # we don't want sign for alpha to change, so flip alpha. (see eqn. 3 in the paper) self._sin_alpha = numpy.sin(-alpha) self._tan_alpha = numpy.tan(-alpha) self._r_ref = r_ref diff --git a/galpy/potential/interpSphericalPotential.py b/galpy/potential/interpSphericalPotential.py index 1cb0579d4..b21cf82b4 100644 --- a/galpy/potential/interpSphericalPotential.py +++ b/galpy/potential/interpSphericalPotential.py @@ -91,7 +91,7 @@ def __init__( # Extrapolate as mass within rgrid[-1] self._rmin = rgrid[0] self._rmax = rgrid[-1] - self._total_mass = -self._rmax**2.0 * self._force_spline(self._rmax) + self._total_mass = -(self._rmax**2.0) * self._force_spline(self._rmax) self._Phimax = ( -self._pot_spline(self._rmax) + self._Phi0 + self._total_mass / self._rmax ) diff --git a/galpy/util/coords.py b/galpy/util/coords.py index 6c0690c1a..47c9113e3 100644 --- a/galpy/util/coords.py +++ b/galpy/util/coords.py @@ -2440,7 +2440,10 @@ def custom_to_radec(phi1, phi2, T=None, degree=False): if T is None: raise ValueError("Must set T= for custom_to_radec") return radec_to_custom( - phi1, phi2, T=numpy.transpose(T), degree=degree # T.T = inv(T) + phi1, + phi2, + T=numpy.transpose(T), + degree=degree, # T.T = inv(T) ) @@ -2475,7 +2478,12 @@ def custom_to_pmrapmdec(pmphi1, pmphi2, phi1, phi2, T=None, degree=False): if T is None: raise ValueError("Must set T= for custom_to_pmrapmdec") return pmrapmdec_to_custom( - pmphi1, pmphi2, phi1, phi2, T=numpy.transpose(T), degree=degree # T.T = inv(T) + pmphi1, + pmphi2, + phi1, + phi2, + T=numpy.transpose(T), + degree=degree, # T.T = inv(T) ) diff --git a/galpy/util/plot.py b/galpy/util/plot.py index 2895a5997..dcd0e3856 100644 --- a/galpy/util/plot.py +++ b/galpy/util/plot.py @@ -837,7 +837,7 @@ def text(*args, **kwargs): xycoords="axes fraction", horizontalalignment="center", verticalalignment="top", - **kwargs + **kwargs, ) elif kwargs.pop("bottom_left", False): pyplot.annotate(args[0], (0.05, 0.05), xycoords="axes fraction", **kwargs) @@ -847,7 +847,7 @@ def text(*args, **kwargs): (0.95, 0.05), xycoords="axes fraction", horizontalalignment="right", - **kwargs + **kwargs, ) elif kwargs.pop("top_right", False): pyplot.annotate( @@ -856,7 +856,7 @@ def text(*args, **kwargs): xycoords="axes fraction", horizontalalignment="right", verticalalignment="top", - **kwargs + **kwargs, ) elif kwargs.pop("top_left", False): pyplot.annotate( @@ -864,7 +864,7 @@ def text(*args, **kwargs): (0.05, 0.95), xycoords="axes fraction", verticalalignment="top", - **kwargs + **kwargs, ) else: pyplot.text(*args, **kwargs) @@ -1108,7 +1108,7 @@ def scatterplot(x, y, *args, **kwargs): overplot=True, color="%.2f" % (1.0 - w8[ii]), *args, - **kwargs + **kwargs, ) else: plot(plotx, ploty, overplot=True, zorder=1, *args, **kwargs) diff --git a/tests/conftest.py b/tests/conftest.py index adcb5d9fc..19c3f12b5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -103,18 +103,18 @@ def pytest_generate_tests(metafunc): jactol["mockFlatDehnenBarPotential"] = -8.0 # these are more difficult jactol["mockFlatDehnenSmoothBarPotential"] = -8.0 # these are more difficult jactol["mockMovingObjectLongIntPotential"] = -8.0 # these are more difficult - jactol["mockSlowFlatEllipticalDiskPotential"] = ( - -6.0 - ) # these are more difficult (and also not quite conserved) - jactol["mockSlowFlatSteadyLogSpiralPotential"] = ( - -8.0 - ) # these are more difficult (and also not quite conserved) - jactol["mockSlowFlatDehnenSmoothBarPotential"] = ( - -8.0 - ) # these are more difficult (and also not quite conserved) - jactol["mockSlowFlatDecayingDehnenSmoothBarPotential"] = ( - -8.0 - ) # these are more difficult (and also not quite conserved) + jactol[ + "mockSlowFlatEllipticalDiskPotential" + ] = -6.0 # these are more difficult (and also not quite conserved) + jactol[ + "mockSlowFlatSteadyLogSpiralPotential" + ] = -8.0 # these are more difficult (and also not quite conserved) + jactol[ + "mockSlowFlatDehnenSmoothBarPotential" + ] = -8.0 # these are more difficult (and also not quite conserved) + jactol[ + "mockSlowFlatDecayingDehnenSmoothBarPotential" + ] = -8.0 # these are more difficult (and also not quite conserved) # Now generate all inputs and run tests tols = [tol[p] if p in tol else tol["default"] for p in pots] jactols = [jactol[p] if p in jactol else tol["default"] for p in pots] diff --git a/tests/test_actionAngle.py b/tests/test_actionAngle.py index 057d6dbc9..2eff953f7 100644 --- a/tests/test_actionAngle.py +++ b/tests/test_actionAngle.py @@ -298,7 +298,7 @@ def _evaluate(self, x, t=0.0): return self._omega**2.0 * x**2.0 / 2.0 def _force(self, x, t=0.0): - return -self._omega**2.0 * x + return -(self._omega**2.0) * x ipz = HO(omega=2.23) aAH = actionAngleHarmonic(omega=ipz._omega) @@ -334,7 +334,7 @@ def _evaluate(self, x, t=0.0): return self._omega**2.0 * x**2.0 / 2.0 def _force(self, x, t=0.0): - return -self._omega**2.0 * x + return -(self._omega**2.0) * x ipz = HO(omega=2.23) aAH = actionAngleHarmonic(omega=ipz._omega) @@ -375,7 +375,7 @@ def _evaluate(self, x, t=0.0): return self._omega**2.0 * x**2.0 / 2.0 def _force(self, x, t=0.0): - return -self._omega**2.0 * x + return -(self._omega**2.0) * x ipz = HO(omega=2.236) aAH = actionAngleHarmonic(omega=ipz._omega) @@ -506,15 +506,13 @@ def test_actionAngleIsochrone_basic_freqs(): # circular orbit R, vR, vT, z, vz = 1.0, 0.0, 1.0, 0.0, 0.0 jos = aAI.actionsFreqs(R, vR, vT, z, vz) - assert ( - numpy.fabs((jos[3] - ip.epifreq(1.0)) / ip.epifreq(1.0)) < 10.0**-12.0 - ), "Circular orbit in the isochrone potential does not have Or=kappa at %g%%" % ( - 100.0 * numpy.fabs((jos[3] - ip.epifreq(1.0)) / ip.epifreq(1.0)) + assert numpy.fabs((jos[3] - ip.epifreq(1.0)) / ip.epifreq(1.0)) < 10.0**-12.0, ( + "Circular orbit in the isochrone potential does not have Or=kappa at %g%%" + % (100.0 * numpy.fabs((jos[3] - ip.epifreq(1.0)) / ip.epifreq(1.0))) ) - assert ( - numpy.fabs((jos[4] - ip.omegac(1.0)) / ip.omegac(1.0)) < 10.0**-12.0 - ), "Circular orbit in the isochrone potential does not have Op=Omega at %g%%" % ( - 100.0 * numpy.fabs((jos[4] - ip.omegac(1.0)) / ip.omegac(1.0)) + assert numpy.fabs((jos[4] - ip.omegac(1.0)) / ip.omegac(1.0)) < 10.0**-12.0, ( + "Circular orbit in the isochrone potential does not have Op=Omega at %g%%" + % (100.0 * numpy.fabs((jos[4] - ip.omegac(1.0)) / ip.omegac(1.0))) ) assert ( numpy.fabs((jos[5] - ip.verticalfreq(1.0)) / ip.verticalfreq(1.0)) < 10.0**-12.0 @@ -1202,21 +1200,18 @@ def test_actionAngleSpherical_otherIsochrone_actions(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-10.0 - ), "actionAngleSpherical applied to isochrone potential fails for Jr at %g%%" % ( - djr * 100.0 + assert djr < 10.0**-10.0, ( + "actionAngleSpherical applied to isochrone potential fails for Jr at %g%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleSpherical applied to isochrone potential fails for Lz at %g%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleSpherical applied to isochrone potential fails for Lz at %g%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-10.0 - ), "actionAngleSpherical applied to isochrone potential fails for Jz at %g%%" % ( - djz * 100.0 + assert djz < 10.0**-10.0, ( + "actionAngleSpherical applied to isochrone potential fails for Jz at %g%%" + % (djz * 100.0) ) return None @@ -1235,20 +1230,17 @@ def test_actionAngleSpherical_otherIsochrone_freqs(): dOr = numpy.fabs((jiO[3] - jiaO[3]) / jiO[3]) dOp = numpy.fabs((jiO[4] - jiaO[4]) / jiO[4]) dOz = numpy.fabs((jiO[5] - jiaO[5]) / jiO[5]) - assert ( - dOr < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Or at %g%%" % ( - dOr * 100.0 + assert dOr < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Or at %g%%" + % (dOr * 100.0) ) - assert ( - dOp < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Op at %g%%" % ( - dOp * 100.0 + assert dOp < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Op at %g%%" + % (dOp * 100.0) ) - assert ( - dOz < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Oz at %g%%" % ( - dOz * 100.0 + assert dOz < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Oz at %g%%" + % (dOz * 100.0) ) return None @@ -1268,20 +1260,17 @@ def test_actionAngleSpherical_otherIsochrone_freqs_fixed_quad(): dOr = numpy.fabs((jiO[3] - jiaO[3]) / jiO[3]) dOp = numpy.fabs((jiO[4] - jiaO[4]) / jiO[4]) dOz = numpy.fabs((jiO[5] - jiaO[5]) / jiO[5]) - assert ( - dOr < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Or at %g%%" % ( - dOr * 100.0 + assert dOr < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Or at %g%%" + % (dOr * 100.0) ) - assert ( - dOp < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Op at %g%%" % ( - dOp * 100.0 + assert dOp < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Op at %g%%" + % (dOp * 100.0) ) - assert ( - dOz < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for Oz at %g%%" % ( - dOz * 100.0 + assert dOz < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for Oz at %g%%" + % (dOz * 100.0) ) return None @@ -1300,20 +1289,17 @@ def test_actionAngleSpherical_otherIsochrone_angles(): dar = numpy.fabs((jiO[6] - jiaO[6]) / jiO[6]) dap = numpy.fabs((jiO[7] - jiaO[7]) / jiO[7]) daz = numpy.fabs((jiO[8] - jiaO[8]) / jiO[8]) - assert ( - dar < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for ar at %g%%" % ( - dar * 100.0 + assert dar < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for ar at %g%%" + % (dar * 100.0) ) - assert ( - dap < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for ap at %g%%" % ( - dap * 100.0 + assert dap < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for ap at %g%%" + % (dap * 100.0) ) - assert ( - daz < 10.0**-6.0 - ), "actionAngleSpherical applied to isochrone potential fails for az at %g%%" % ( - daz * 100.0 + assert daz < 10.0**-6.0, ( + "actionAngleSpherical applied to isochrone potential fails for az at %g%%" + % (daz * 100.0) ) return None @@ -1776,21 +1762,18 @@ def test_actionAngleAdiabatic_Isochrone_actions(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-1.2 - ), "actionAngleAdiabatic applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-1.2, ( + "actionAngleAdiabatic applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleAdiabatic applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleAdiabatic applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-1.2 - ), "actionAngleAdiabatic applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-1.2, ( + "actionAngleAdiabatic applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -1917,21 +1900,18 @@ def test_actionAngleAdiabaticGrid_Isochrone_actions(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-1.2 - ), "actionAngleAdiabatic applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-1.2, ( + "actionAngleAdiabatic applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleAdiabatic applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleAdiabatic applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-1.2 - ), "actionAngleAdiabatic applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-1.2, ( + "actionAngleAdiabatic applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -2168,11 +2148,11 @@ def test_actionAngleStaeckel_actions_order(): jrt, jpt, jzt = aAS(o, order=10000, fixed_quad=True) jr1, jp1, jz1 = aAS(o, order=5, fixed_quad=True) jr2, jp2, jz2 = aAS(o, order=50, fixed_quad=True) - assert numpy.fabs(jr1 - jrt) > numpy.fabs( - jr2 - jrt + assert ( + numpy.fabs(jr1 - jrt) > numpy.fabs(jr2 - jrt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(jz1 - jzt) > numpy.fabs( - jz2 - jzt + assert ( + numpy.fabs(jz1 - jzt) > numpy.fabs(jz2 - jzt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" return None @@ -2189,11 +2169,11 @@ def test_actionAngleStaeckel_actions_order_c(): jrt, jpt, jzt = aAS(o, order=10000) jr1, jp1, jz1 = aAS(o, order=5) jr2, jp2, jz2 = aAS(o, order=50) - assert numpy.fabs(jr1 - jrt) > numpy.fabs( - jr2 - jrt + assert ( + numpy.fabs(jr1 - jrt) > numpy.fabs(jr2 - jrt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(jz1 - jzt) > numpy.fabs( - jz2 - jzt + assert ( + numpy.fabs(jz1 - jzt) > numpy.fabs(jz2 - jzt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" return None @@ -2486,20 +2466,20 @@ def test_actionAngleStaeckel_freqs_order_c(): jrt, jpt, jzt, ort, opt, ozt = aAS.actionsFreqs(o, order=10000) jr1, jp1, jz1, or1, op1, oz1 = aAS.actionsFreqs(o, order=5) jr2, jp2, jz2, or2, op2, oz2 = aAS.actionsFreqs(o, order=50) - assert numpy.fabs(jr1 - jrt) > numpy.fabs( - jr2 - jrt + assert ( + numpy.fabs(jr1 - jrt) > numpy.fabs(jr2 - jrt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(jz1 - jzt) > numpy.fabs( - jz2 - jzt + assert ( + numpy.fabs(jz1 - jzt) > numpy.fabs(jz2 - jzt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(or1 - ort) > numpy.fabs( - or2 - ort + assert ( + numpy.fabs(or1 - ort) > numpy.fabs(or2 - ort) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(op1 - opt) > numpy.fabs( - op2 - opt + assert ( + numpy.fabs(op1 - opt) > numpy.fabs(op2 - opt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(oz1 - ozt) > numpy.fabs( - oz2 - ozt + assert ( + numpy.fabs(oz1 - ozt) > numpy.fabs(oz2 - ozt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" return None @@ -2555,29 +2535,29 @@ def test_actionAngleStaeckel_angles_order_c(): jrt, jpt, jzt, ort, opt, ozt, art, apt, azt = aAS.actionsFreqsAngles(o, order=10000) jr1, jp1, jz1, or1, op1, oz1, ar1, ap1, az1 = aAS.actionsFreqsAngles(o, order=5) jr2, jp2, jz2, or2, op2, oz2, ar2, ap2, az2 = aAS.actionsFreqsAngles(o, order=50) - assert numpy.fabs(jr1 - jrt) > numpy.fabs( - jr2 - jrt + assert ( + numpy.fabs(jr1 - jrt) > numpy.fabs(jr2 - jrt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(jz1 - jzt) > numpy.fabs( - jz2 - jzt + assert ( + numpy.fabs(jz1 - jzt) > numpy.fabs(jz2 - jzt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(or1 - ort) > numpy.fabs( - or2 - ort + assert ( + numpy.fabs(or1 - ort) > numpy.fabs(or2 - ort) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(op1 - opt) > numpy.fabs( - op2 - opt + assert ( + numpy.fabs(op1 - opt) > numpy.fabs(op2 - opt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(oz1 - ozt) > numpy.fabs( - oz2 - ozt + assert ( + numpy.fabs(oz1 - ozt) > numpy.fabs(oz2 - ozt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(ar1 - art) > numpy.fabs( - ar2 - art + assert ( + numpy.fabs(ar1 - art) > numpy.fabs(ar2 - art) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(ap1 - apt) > numpy.fabs( - ap2 - apt + assert ( + numpy.fabs(ap1 - apt) > numpy.fabs(ap2 - apt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" - assert numpy.fabs(az1 - azt) > numpy.fabs( - az2 - azt + assert ( + numpy.fabs(az1 - azt) > numpy.fabs(az2 - azt) ), "Accuracy of actionAngleStaeckel does not increase with increasing order of integration" return None @@ -3522,21 +3502,18 @@ def test_actionAngleStaeckel_otherIsochrone_actions(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -3559,21 +3536,18 @@ def test_actionAngleStaeckel_otherIsochrone_actions_fixed_quad(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0])[0] dlz = numpy.fabs((ji[1] - jia[1]) / ji[1])[0] djz = numpy.fabs((ji[2] - jia[2]) / ji[2])[0] - assert ( - djr < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -3596,21 +3570,18 @@ def test_actionAngleStaeckel_otherIsochrone_actions_c(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-3.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-3.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -3629,20 +3600,17 @@ def test_actionAngleStaeckel_otherIsochrone_freqs(): dOr = numpy.fabs((jiO[3] - jiaO[3]) / jiO[3]) dOp = numpy.fabs((jiO[4] - jiaO[4]) / jiO[4]) dOz = numpy.fabs((jiO[5] - jiaO[5]) / jiO[5]) - assert ( - dOr < 10.0**-5.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Or at %g%%" % ( - dOr * 100.0 + assert dOr < 10.0**-5.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Or at %g%%" + % (dOr * 100.0) ) - assert ( - dOp < 10.0**-5.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Op at %g%%" % ( - dOp * 100.0 + assert dOp < 10.0**-5.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Op at %g%%" + % (dOp * 100.0) ) - assert ( - dOz < 1.5 * 10.0**-4.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Oz at %g%%" % ( - dOz * 100.0 + assert dOz < 1.5 * 10.0**-4.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Oz at %g%%" + % (dOz * 100.0) ) return None @@ -3661,20 +3629,17 @@ def test_actionAngleStaeckel_otherIsochrone_angles(): dar = numpy.fabs((jiO[6] - jiaO[6]) / jiO[6]) dap = numpy.fabs((jiO[7] - jiaO[7]) / jiO[7]) daz = numpy.fabs((jiO[8] - jiaO[8]) / jiO[8]) - assert ( - dar < 10.0**-4.0 - ), "actionAngleStaeckel applied to isochrone potential fails for ar at %g%%" % ( - dar * 100.0 + assert dar < 10.0**-4.0, ( + "actionAngleStaeckel applied to isochrone potential fails for ar at %g%%" + % (dar * 100.0) ) - assert ( - dap < 10.0**-6.0 - ), "actionAngleStaeckel applied to isochrone potential fails for ap at %g%%" % ( - dap * 100.0 + assert dap < 10.0**-6.0, ( + "actionAngleStaeckel applied to isochrone potential fails for ap at %g%%" + % (dap * 100.0) ) - assert ( - daz < 10.0**-4.0 - ), "actionAngleStaeckel applied to isochrone potential fails for az at %g%%" % ( - daz * 100.0 + assert daz < 10.0**-4.0, ( + "actionAngleStaeckel applied to isochrone potential fails for az at %g%%" + % (daz * 100.0) ) return None @@ -3834,21 +3799,18 @@ def test_actionAngleStaeckelGrid_Isochrone_actions(): djr = numpy.fabs((ji[0] - jia[0]) / ji[0]) dlz = numpy.fabs((ji[1] - jia[1]) / ji[1]) djz = numpy.fabs((ji[2] - jia[2]) / ji[2]) - assert ( - djr < 10.0**-1.2 - ), "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" % ( - djr * 100.0 + assert djr < 10.0**-1.2, ( + "actionAngleStaeckel applied to isochrone potential fails for Jr at %f%%" + % (djr * 100.0) ) # Lz and Jz are easy, because ip is a spherical potential - assert ( - dlz < 10.0**-10.0 - ), "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" % ( - dlz * 100.0 + assert dlz < 10.0**-10.0, ( + "actionAngleStaeckel applied to isochrone potential fails for Lz at %f%%" + % (dlz * 100.0) ) - assert ( - djz < 10.0**-1.2 - ), "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" % ( - djz * 100.0 + assert djz < 10.0**-1.2, ( + "actionAngleStaeckel applied to isochrone potential fails for Jz at %f%%" + % (djz * 100.0) ) return None @@ -4404,20 +4366,17 @@ def test_actionAngleIsochroneApprox_bovy14(): maxdj = numpy.amax( numpy.fabs(js - numpy.tile(numpy.mean(js, axis=1), (len(times), 1)).T), axis=1 ) / numpy.mean(js, axis=1) - assert ( - maxdj[0] < 3.0 * 10.0**-2.0 - ), "Jr conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" % ( - 100.0 * maxdj[0] + assert maxdj[0] < 3.0 * 10.0**-2.0, ( + "Jr conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" + % (100.0 * maxdj[0]) ) - assert ( - maxdj[1] < 10.0**-2.0 - ), "Lz conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" % ( - 100.0 * maxdj[1] + assert maxdj[1] < 10.0**-2.0, ( + "Lz conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" + % (100.0 * maxdj[1]) ) - assert ( - maxdj[2] < 2.0 * 10.0**-2.0 - ), "Jz conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" % ( - 100.0 * maxdj[2] + assert maxdj[2] < 2.0 * 10.0**-2.0, ( + "Jz conservation for the GD-1 like orbit of Bovy (2014) fails at %f%%" + % (100.0 * maxdj[2]) ) return None @@ -4748,8 +4707,8 @@ class TwoPowerSphericalPotentialNoR2deriv(TwoPowerSphericalPotential): pytest.fail( "TwoPowerSphericalPotentialNoR2deriv appears to now have second derivatives, means that it cannot be used to test exceptions based on not having the second derivatives any longer" ) - assert "second derivatives" in str( - excinfo.value + assert ( + "second derivatives" in str(excinfo.value) ), "Estimating delta for potential lacking second derivatives should have failed with a message about the lack of second derivatives" # Generic non-axi sp = SpiralArmsPotential() @@ -4758,8 +4717,8 @@ class TwoPowerSphericalPotentialNoR2deriv(TwoPowerSphericalPotential): pytest.fail( "TwoPowerSphericalPotentialNoR2deriv appears to now have second derivatives, means that it cannot be used to test exceptions based on not having the second derivatives any longer" ) - assert "not axisymmetric" in str( - excinfo.value + assert ( + "not axisymmetric" in str(excinfo.value) ), "Estimating delta for a non-axi potential should have failed with a message about the fact that the potential is non-axisymmetric" return None @@ -4794,8 +4753,8 @@ class TwoPowerSphericalPotentialNoR2deriv(TwoPowerSphericalPotential): pytest.fail( "TwoPowerSphericalPotentialNoR2deriv appears to now have second derivatives, means that it cannot be used to test exceptions based on not having the second derivatives any longer" ) - assert "second derivatives" in str( - excinfo.value + assert ( + "second derivatives" in str(excinfo.value) ), "Estimating delta for potential lacking second derivatives should have failed with a message about the lack of second derivatives" # Generic non-axi sp = SpiralArmsPotential() @@ -4804,8 +4763,8 @@ class TwoPowerSphericalPotentialNoR2deriv(TwoPowerSphericalPotential): pytest.fail( "SpiralArms appears to now have second derivatives, means that it cannot be used to test exceptions based on not having the second derivatives any longer" ) - assert "not axisymmetric" in str( - excinfo.value + assert ( + "not axisymmetric" in str(excinfo.value) ), "Estimating delta for a non-axi potential should have failed with a message about the fact that the potential is non-axisymmetric" return None @@ -5067,9 +5026,7 @@ def test_MWPotential_warning_adiabatic(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleAdiabatic with MWPotential should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleAdiabatic with MWPotential should have thrown a warning, but didn't" # Grid with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -5085,9 +5042,7 @@ def test_MWPotential_warning_adiabatic(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleAdiabaticGrid with MWPotential should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleAdiabaticGrid with MWPotential should have thrown a warning, but didn't" return None @@ -5110,9 +5065,7 @@ def test_MWPotential_warning_staeckel(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleStaeckel with MWPotential should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleStaeckel with MWPotential should have thrown a warning, but didn't" # Grid with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -5126,9 +5079,7 @@ def test_MWPotential_warning_staeckel(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleStaeckelGrid with MWPotential should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleStaeckelGrid with MWPotential should have thrown a warning, but didn't" return None @@ -5151,9 +5102,7 @@ def test_MWPotential_warning_isochroneapprox(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleIsochroneApprox with MWPotential should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleIsochroneApprox with MWPotential should have thrown a warning, but didn't" return None @@ -6182,9 +6131,7 @@ def test_actionAngleVerticalInverse_convergence_warnings(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleVerticalInverse for large energy should have raised convergence warning, but didn't" + assert raisedWarning, "actionAngleVerticalInverse for large energy should have raised convergence warning, but didn't" for wa in w: raisedWarning = ( str(wa.message) @@ -6192,9 +6139,7 @@ def test_actionAngleVerticalInverse_convergence_warnings(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleVerticalInverse for large energy should have raised convergence warning, but didn't" + assert raisedWarning, "actionAngleVerticalInverse for large energy should have raised convergence warning, but didn't" return None @@ -6549,25 +6494,21 @@ def check_actionAngle_conserved_EccZmaxRperiRap( es, zmaxs, rperis, raps = aA.EccZmaxRperiRap( obs.R(times), obs.vR(times), obs.vT(times), obs.z(times), obs.vz(times) ) - assert ( - numpy.amax(numpy.fabs(es / numpy.mean(es) - 1)) < 10.0**tole - ), "Eccentricity conservation fails at %g%%" % ( - 100.0 * numpy.amax(numpy.fabs(es / numpy.mean(es) - 1)) + assert numpy.amax(numpy.fabs(es / numpy.mean(es) - 1)) < 10.0**tole, ( + "Eccentricity conservation fails at %g%%" + % (100.0 * numpy.amax(numpy.fabs(es / numpy.mean(es) - 1))) ) - assert ( - numpy.amax(numpy.fabs(zmaxs / numpy.mean(zmaxs) - 1)) < 10.0**tolzmax - ), "Zmax conservation fails at %g%%" % ( - 100.0 * numpy.amax(numpy.fabs(zmaxs / numpy.mean(zmaxs) - 1)) + assert numpy.amax(numpy.fabs(zmaxs / numpy.mean(zmaxs) - 1)) < 10.0**tolzmax, ( + "Zmax conservation fails at %g%%" + % (100.0 * numpy.amax(numpy.fabs(zmaxs / numpy.mean(zmaxs) - 1))) ) - assert ( - numpy.amax(numpy.fabs(rperis / numpy.mean(rperis) - 1)) < 10.0**tolrperi - ), "Rperi conservation fails at %g%%" % ( - 100.0 * numpy.amax(numpy.fabs(rperis / numpy.mean(rperis) - 1)) + assert numpy.amax(numpy.fabs(rperis / numpy.mean(rperis) - 1)) < 10.0**tolrperi, ( + "Rperi conservation fails at %g%%" + % (100.0 * numpy.amax(numpy.fabs(rperis / numpy.mean(rperis) - 1))) ) - assert ( - numpy.amax(numpy.fabs(raps / numpy.mean(raps) - 1)) < 10.0**tolrap - ), "Rap conservation fails at %g%%" % ( - 100.0 * numpy.amax(numpy.fabs(raps / numpy.mean(raps) - 1)) + assert numpy.amax(numpy.fabs(raps / numpy.mean(raps) - 1)) < 10.0**tolrap, ( + "Rap conservation fails at %g%%" + % (100.0 * numpy.amax(numpy.fabs(raps / numpy.mean(raps) - 1))) ) return None diff --git a/tests/test_actionAngleTorus.py b/tests/test_actionAngleTorus.py index 5e80df519..2e3ced5df 100644 --- a/tests/test_actionAngleTorus.py +++ b/tests/test_actionAngleTorus.py @@ -741,9 +741,7 @@ def test_actionAngleTorus_AutoFitWarning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) aAT.xvFreqs(jr, jp, jz, ar, ap, az) @@ -756,9 +754,7 @@ def test_actionAngleTorus_AutoFitWarning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) aAT.Freqs(jr, jp, jz) @@ -771,9 +767,7 @@ def test_actionAngleTorus_AutoFitWarning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) aAT.hessianFreqs(jr, jp, jz) @@ -786,9 +780,7 @@ def test_actionAngleTorus_AutoFitWarning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) aAT.xvJacobianFreqs(jr, jp, jz, ar, ap, az) @@ -801,9 +793,7 @@ def test_actionAngleTorus_AutoFitWarning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" + assert raisedWarning, "actionAngleTorus with flattened LogarithmicHaloPotential and a particular orbit should have thrown a warning, but didn't" return None diff --git a/tests/test_coords.py b/tests/test_coords.py index e22d67a79..777507631 100644 --- a/tests/test_coords.py +++ b/tests/test_coords.py @@ -2540,7 +2540,7 @@ def test_vRvz_to_pupv_oblate(): *coords.Rz_to_uv(R, z, delta=delta, oblate=True), delta=delta, oblate=True, - uv=True + uv=True, )[0] ) < 10.0**-3.0 @@ -2554,7 +2554,7 @@ def test_vRvz_to_pupv_oblate(): *coords.Rz_to_uv(R, z, delta=delta, oblate=True), delta=delta, oblate=True, - uv=True + uv=True, )[1] ) < 10.0**-3.0 diff --git a/tests/test_dynamfric.py b/tests/test_dynamfric.py index 788f4db64..03bdec6e0 100644 --- a/tests/test_dynamfric.py +++ b/tests/test_dynamfric.py @@ -263,7 +263,9 @@ def test_dynamfric_c(): normalize=0.03, sigma=4.0, b=0.8, c=1.5, pa=3.0, zvec=[1.0, 0.0, 0.0] ), potential.SCFPotential( - Acos=numpy.array([[[1.0]]]), normalize=1.0, a=3.5 # same as Hernquist + Acos=numpy.array([[[1.0]]]), + normalize=1.0, + a=3.5, # same as Hernquist ), potential.SCFPotential( Acos=numpy.array([[[1.0, 0.0], [0.3, 0.0]]]), # nonaxi @@ -387,7 +389,5 @@ def test_dynamfric_c_minr_warning(): str(rec.message.args[0]) == "Orbit integration with ChandrasekharDynamicalFrictionForce entered domain where r < minr and ChandrasekharDynamicalFrictionForce is turned off; initialize ChandrasekharDynamicalFrictionForce with a smaller minr to avoid this if you wish (but note that you want to turn it off close to the center for an object that sinks all the way to r=0, to avoid numerical instabilities)" ) - assert ( - raisedWarning - ), "Integrating an orbit that goes to r < minr with dynamical friction should have raised a warning, but didn't" + assert raisedWarning, "Integrating an orbit that goes to r < minr with dynamical friction should have raised a warning, but didn't" return None diff --git a/tests/test_evolveddiskdf.py b/tests/test_evolveddiskdf.py index 8c55aabb4..024b462c5 100644 --- a/tests/test_evolveddiskdf.py +++ b/tests/test_evolveddiskdf.py @@ -551,7 +551,9 @@ def test_mildnonaxi_vertexdev_grid(): ) assert ( numpy.fabs(vdev) < 2.0 / 180.0 * numpy.pi - ), "vertexdev of evolveddiskdf for axisymmetric potential is not close to zero" # 2 is pretty big, but the weak spiral creates that + ), ( + "vertexdev of evolveddiskdf for axisymmetric potential is not close to zero" + ) # 2 is pretty big, but the weak spiral creates that vdev = edf.vertexdev( 0.9, phi=0.2, integrate_method="rk6_c", grid=grid, gridpoints=_GRIDPOINTS ) diff --git a/tests/test_galpypaper.py b/tests/test_galpypaper.py index 43d0f4aa1..baa58260f 100644 --- a/tests/test_galpypaper.py +++ b/tests/test_galpypaper.py @@ -1,6 +1,6 @@ -""" Test that all of the examples in the galpy paper run +"""Test that all of the examples in the galpy paper run - isort:skip_file +isort:skip_file """ import os @@ -237,9 +237,7 @@ def test_potentialAPIChange_warning(): ) if raisedWarning: break - assert ( - raisedWarning - ), "Importing galpy.potential does not raise warning about evaluatePotentials API change" + assert raisedWarning, "Importing galpy.potential does not raise warning about evaluatePotentials API change" return None @@ -720,7 +718,10 @@ def test_qdf(): numpy.fabs(df.estimate_hz(0.9, 0.02) - 0.064836202345657207) < 10.0**-4.0 ), "qdf does not behave as expected" # Calculate the mean velocities - df.meanvR(0.9, 0.05), df.meanvT(0.9, 0.05), + ( + df.meanvR(0.9, 0.05), + df.meanvT(0.9, 0.05), + ) df.meanvz(0.9, 0.05) assert ( numpy.fabs(df.meanvR(0.9, 0.05) - 3.8432265354618213e-18) < 10.0**-4.0 diff --git a/tests/test_interp_potential.py b/tests/test_interp_potential.py index 815a0c158..c17f7973b 100644 --- a/tests/test_interp_potential.py +++ b/tests/test_interp_potential.py @@ -1668,7 +1668,9 @@ def test_interpolation_potential_epifreq_smalln(): / potential.epifreq(potential.MWPotential, rs) ) < 10.0**-2.0 - ), "RZPot interpolation of epifreq w/ interpRZPotential fails for vector input" # not as harsh, bc we don't have many points + ), ( + "RZPot interpolation of epifreq w/ interpRZPotential fails for vector input" + ) # not as harsh, bc we don't have many points rzpot = potential.interpRZPotential( RZPot=potential.MWPotential, rgrid=(numpy.log(1.0), numpy.log(1.3), 3), @@ -1683,7 +1685,9 @@ def test_interpolation_potential_epifreq_smalln(): / potential.epifreq(potential.MWPotential, rs) ) < 10.0**-2.0 - ), "RZPot interpolation of epifreq w/ interpRZPotential fails for vector input" # not as harsh, bc we don't have many points + ), ( + "RZPot interpolation of epifreq w/ interpRZPotential fails for vector input" + ) # not as harsh, bc we don't have many points return None diff --git a/tests/test_orbit.py b/tests/test_orbit.py index b3a3afff9..a1ca47339 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -183,34 +183,36 @@ def test_energy_jacobi_conservation(pot, ttol, tjactol, firstTest): if firstTest or "testMWPotential" in pot: # Some basic checking of the energy and Jacobi functions assert ( - o.E(pot=None) - o.E(pot=tp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" + (o.E(pot=None) - o.E(pot=tp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" assert ( - o.E() - o.E(0.0) - ) ** 2.0 < 10.0**ttol, ( - "Energy calculated with o.E() and o.E(0.) do not agree" - ) + (o.E() - o.E(0.0)) ** 2.0 < 10.0** ttol + ), "Energy calculated with o.E() and o.E(0.) do not agree" assert ( - o.Jacobi(OmegaP=None) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, ( - "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" - ) + (o.Jacobi(OmegaP=None) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=tp) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=tp)) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=[tp]) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=[the Potential the orbit was integrated with] do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=[tp])) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=[the Potential the orbit was integrated with] do not agree" if not tp.isNonAxi: assert ( - o.Jacobi(OmegaP=1.0) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" + (o.Jacobi(OmegaP=1.0) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" assert ( - o.Jacobi(OmegaP=[0.0, 0.0, 1.0]) - o.Jacobi(OmegaP=1.0) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=[0,0,1] for axisymmetric potential is not equal to o.Jacobi with OmegaP=1" + (o.Jacobi(OmegaP=[0.0, 0.0, 1.0]) - o.Jacobi(OmegaP=1.0)) ** 2.0 + < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=[0,0,1] for axisymmetric potential is not equal to o.Jacobi with OmegaP=1" assert ( - o.Jacobi(OmegaP=numpy.array([0.0, 0.0, 1.0])) - o.Jacobi(OmegaP=1.0) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=[0,0,1] for axisymmetric potential is not equal to o.Jacobi with OmegaP=1" + ( + o.Jacobi(OmegaP=numpy.array([0.0, 0.0, 1.0])) + - o.Jacobi(OmegaP=1.0) + ) + ** 2.0 + < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=[0,0,1] for axisymmetric potential is not equal to o.Jacobi with OmegaP=1" o = setup_orbit_energy(tp, axi=False, henon="Henon" in pot) try: o.E() @@ -265,27 +267,23 @@ def test_energy_jacobi_conservation(pot, ttol, tjactol, firstTest): if firstTest or "MWPotential" in pot: # Some basic checking of the energy function assert ( - o.E(pot=None) - o.E(pot=tp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" + (o.E(pot=None) - o.E(pot=tp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" assert ( - o.E() - o.E(0.0) - ) ** 2.0 < 10.0**ttol, ( - "Energy calculated with o.E() and o.E(0.) do not agree" - ) + (o.E() - o.E(0.0)) ** 2.0 < 10.0** ttol + ), "Energy calculated with o.E() and o.E(0.) do not agree" assert ( - o.Jacobi(OmegaP=None) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, ( - "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" - ) + (o.Jacobi(OmegaP=None) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=tp) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=tp)) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=[tp]) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=[tp])) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(OmegaP=1.0) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" + (o.Jacobi(OmegaP=1.0) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" o = setup_orbit_energy(tp, axi=True, henon="Henon" in pot) try: o.E() @@ -340,30 +338,26 @@ def test_energy_jacobi_conservation(pot, ttol, tjactol, firstTest): if firstTest or "MWPotential" in pot: # Some basic checking of the energy function assert ( - o.E(pot=None) - o.E(pot=ptp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the planarPotential the orbit was integrated with do not agree for planarPotential" + (o.E(pot=None) - o.E(pot=ptp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the planarPotential the orbit was integrated with do not agree for planarPotential" assert ( - o.E(pot=None) - o.E(pot=tp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree for planarPotential" + (o.E(pot=None) - o.E(pot=tp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree for planarPotential" assert ( - o.E() - o.E(0.0) - ) ** 2.0 < 10.0**ttol, ( - "Energy calculated with o.E() and o.E(0.) do not agree" - ) + (o.E() - o.E(0.0)) ** 2.0 < 10.0** ttol + ), "Energy calculated with o.E() and o.E(0.) do not agree" assert ( - o.Jacobi(OmegaP=None) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, ( - "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" - ) + (o.Jacobi(OmegaP=None) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=tp) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=tp)) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=[tp]) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=[tp])) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(OmegaP=1.0) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" + (o.Jacobi(OmegaP=1.0) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" o = setup_orbit_energy(ptp, axi=True) try: o.E() @@ -413,30 +407,26 @@ def test_energy_jacobi_conservation(pot, ttol, tjactol, firstTest): if firstTest or "MWPotential" in pot: # Some basic checking of the energy function assert ( - o.E(pot=None) - o.E(pot=ptp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the planarPotential the orbit was integrated with do not agree for planarPotential" + (o.E(pot=None) - o.E(pot=ptp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the planarPotential the orbit was integrated with do not agree for planarPotential" assert ( - o.E(pot=None) - o.E(pot=tp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree for planarPotential" + (o.E(pot=None) - o.E(pot=tp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree for planarPotential" assert ( - o.E() - o.E(0.0) - ) ** 2.0 < 10.0**ttol, ( - "Energy calculated with o.E() and o.E(0.) do not agree" - ) + (o.E() - o.E(0.0)) ** 2.0 < 10.0** ttol + ), "Energy calculated with o.E() and o.E(0.) do not agree" assert ( - o.Jacobi(OmegaP=None) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, ( - "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" - ) + (o.Jacobi(OmegaP=None) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=None is not equal to o.Jacobi" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=tp) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=tp)) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(pot=None) - o.Jacobi(pot=[tp]) - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" + (o.Jacobi(pot=None) - o.Jacobi(pot=[tp])) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with pot=None is not equal to o.Jacobi with pot=the Potential the orbit was integrated with do not agree" assert ( - o.Jacobi(OmegaP=1.0) - o.Jacobi() - ) ** 2.0 < 10.0**ttol, "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" + (o.Jacobi(OmegaP=1.0) - o.Jacobi()) ** 2.0 < 10.0** ttol + ), "o.Jacobi calculated with OmegaP=1. for axisymmetric potential is not equal to o.Jacobi (OmegaP=1 is the default for potentials without a pattern speed" o = setup_orbit_energy(ptp, axi=False) try: o.E() @@ -597,13 +587,11 @@ def test_energy_conservation_linear(pot, ttol, firstTest): if firstTest or "testMWPotential" in pot or "linearMWPotential" in pot: # Some basic checking of the energy function assert ( - o.E(pot=None) - o.E(pot=tp) - ) ** 2.0 < 10.0**ttol, "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" + (o.E(pot=None) - o.E(pot=tp)) ** 2.0 < 10.0** ttol + ), "Energy calculated with pot=None and pot=the Potential the orbit was integrated with do not agree" assert ( - o.E() - o.E(0.0) - ) ** 2.0 < 10.0**ttol, ( - "Energy calculated with o.E() and o.E(0.) do not agree" - ) + (o.E() - o.E(0.0)) ** 2.0 < 10.0** ttol + ), "Energy calculated with o.E() and o.E(0.) do not agree" o = setup_orbit_energy(tp, axi=False, henon="Henon" in pot) try: o.E() @@ -935,10 +923,8 @@ def test_integrate_SOS_3D(): o.integrate_SOS(psis, pot, method=method) Es = o.E(o.t) assert ( - numpy.std(Es) / numpy.mean(Es) - ) ** 2.0 < 10.0**-10, ( - f"Energy is not conserved by integrate_sos for method={method}" - ) + (numpy.std(Es) / numpy.mean(Es)) ** 2.0 < 10.0** -10 + ), f"Energy is not conserved by integrate_sos for method={method}" return None @@ -996,8 +982,8 @@ def test_integrate_SOS_2D(): o.integrate_SOS(psis, pot, method=method) # default is surface='x' Es = o.E(o.t) assert ( - numpy.std(Es) / numpy.mean(Es) - ) ** 2.0 < 10.0**-10, f"Energy is not conserved by integrate_sos for method={method} and surface={surface}" + (numpy.std(Es) / numpy.mean(Es)) ** 2.0 < 10.0** -10 + ), f"Energy is not conserved by integrate_sos for method={method} and surface={surface}" return None @@ -3398,17 +3384,16 @@ def test_orbit_setup(): o = Orbit([120.0, 60.0, 0.0, 10.0, 20.0, 30.0], uvw=True, lb=True, zo=0.0) assert ( numpy.fabs(o.dist() - 0.0) < 10.0**-2.0 - ), "Orbit dist setup does not agree with o.dist()" # because of tweak in the code to deal with at the Sun + ), ( + "Orbit dist setup does not agree with o.dist()" + ) # because of tweak in the code to deal with at the Sun assert ( - o.U() ** 2.0 + o.V() ** 2.0 + o.W() ** 2.0 - 10.0**2.0 - 20.0**2.0 - 30.0**2.0 - ) < 10.0**-10.0, ( - "Velocity wrt the Sun when looking at Orbit at the Sun does not agree" - ) + (o.U() ** 2.0 + o.V() ** 2.0 + o.W() ** 2.0 - 10.0**2.0 - 20.0**2.0 - 30.0**2.0) + < 10.0** -10.0 + ), "Velocity wrt the Sun when looking at Orbit at the Sun does not agree" assert ( - o.vlos() ** 2.0 - 10.0**2.0 - 20.0**2.0 - 30.0**2.0 - ) < 10.0**-10.0, ( - "Velocity wrt the Sun when looking at Orbit at the Sun does not agree" - ) + (o.vlos() ** 2.0 - 10.0**2.0 - 20.0**2.0 - 30.0**2.0) < 10.0** -10.0 + ), "Velocity wrt the Sun when looking at Orbit at the Sun does not agree" # lb w/ default and UVW o = Orbit([120.0, 60.0, 2.0, -10.0, 20.0, -25.0], lb=True, uvw=True) assert ( @@ -3733,9 +3718,7 @@ def test_orbit_setup_SkyCoord(): str(rec.message.args[0]) == "Orbit's initialization normalization ro and zo are incompatible with SkyCoord's galcen_distance (should have galcen_distance^2 = ro^2 + zo^2)" ) - assert ( - raisedWarning - ), "Orbit initialization with SkyCoord with galcen_distance incompatible with ro should have raised a warning, but didn't" + assert raisedWarning, "Orbit initialization with SkyCoord with galcen_distance incompatible with ro should have raised a warning, but didn't" # If ro and galcen_distance are both specified, don't warn if they *are* consistent (issue #370) c = apycoords.SkyCoord( ra=ra, @@ -4093,11 +4076,11 @@ def test_flip(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -4173,11 +4156,11 @@ def test_flip_inplace(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -4266,11 +4249,11 @@ def test_flip_inplace_integrated(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -4365,11 +4348,11 @@ def test_flip_inplace_integrated_evaluated(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -4512,26 +4495,18 @@ def test_newOrbit(): assert ( numpy.fabs(nos[0].phi() - o.phi(ts[-2])) < 10.0**-10.0 ), "New orbit formed from calling an old orbit does not have the correct phi" - assert not nos[ - 0 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) + assert ( + not nos[0]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" # Second t assert ( numpy.fabs(nos[1].R() - o.R(ts[-1])) < 10.0**-10.0 @@ -4551,26 +4526,18 @@ def test_newOrbit(): assert ( numpy.fabs(nos[1].phi() - o.phi(ts[-1])) < 10.0**-10.0 ), "New orbit formed from calling an old orbit does not have the correct phi" - assert not nos[ - 1 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) + assert ( + not nos[1]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" return None @@ -4676,26 +4643,18 @@ def test_newOrbit_badinterpolation(): assert ( numpy.fabs(nos[0].phi() - o.phi(ts[-2])) < 10.0**-10.0 ), "New orbit formed from calling an old orbit does not have the correct phi" - assert not nos[ - 0 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 0 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) + assert ( + not nos[0]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[0]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" # Second t assert ( numpy.fabs(nos[1].R() - o.R(ts[-1])) < 10.0**-10.0 @@ -4715,26 +4674,18 @@ def test_newOrbit_badinterpolation(): assert ( numpy.fabs(nos[1].phi() - o.phi(ts[-1])) < 10.0**-10.0 ), "New orbit formed from calling an old orbit does not have the correct phi" - assert not nos[ - 1 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._roSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) - assert not nos[ - 1 - ]._voSet, ( - "New orbit formed from calling an old orbit does not have the correct roSet" - ) + assert ( + not nos[1]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._roSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" + assert ( + not nos[1]._voSet + ), "New orbit formed from calling an old orbit does not have the correct roSet" # Try point in between, shouldn't work try: no = o(0.6) @@ -7420,9 +7371,7 @@ def test_orbitint_dissipativefallback(): str(rec.message.args[0]) == "Cannot use symplectic integration because some of the included forces are dissipative (using non-symplectic integrator odeint instead)" ) - assert ( - raisedWarning - ), "Orbit integration with symplectic integrator for dissipative force did not raise fallback warning" + assert raisedWarning, "Orbit integration with symplectic integrator for dissipative force did not raise fallback warning" return None diff --git a/tests/test_orbits.py b/tests/test_orbits.py index 1c009502d..bedfa68a8 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -6557,11 +6557,11 @@ def test_flip(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -6637,11 +6637,11 @@ def test_flip_inplace(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -6730,11 +6730,11 @@ def test_flip_inplace_integrated(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( @@ -6829,11 +6829,11 @@ def test_flip_inplace_integrated_evaluated(): numpy.fabs(o._vo - of._vo) < 10.0**-15.0 ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" if ii == 4: - assert (o._zo is None) * ( - of._zo is None + assert ( + (o._zo is None) * (of._zo is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" - assert (o._solarmotion is None) * ( - of._solarmotion is None + assert ( + (o._solarmotion is None) * (of._solarmotion is None) ), "o.flip() did not conserve physical scales and coordinate-transformation parameters" else: assert ( diff --git a/tests/test_potential.py b/tests/test_potential.py index d019e6e38..f70257e7e 100644 --- a/tests/test_potential.py +++ b/tests/test_potential.py @@ -300,8 +300,8 @@ def test_forceAsDeriv_potential(): ), f"Calculation of the Radial force as the Radial derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRforce-mpotderivR):e}, rel. diff = {numpy.fabs((tRforce-mpotderivR)/tRforce):e}" else: assert ( - tRforce - mpotderivR - ) ** 2.0 / tRforce**2.0 < 10.0**ttol, f"Calculation of the Radial force as the Radial derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRforce-mpotderivR):e}, rel. diff = {numpy.fabs((tRforce-mpotderivR)/tRforce):e}" + (tRforce - mpotderivR) ** 2.0 / tRforce** 2.0 < 10.0** ttol + ), f"Calculation of the Radial force as the Radial derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRforce-mpotderivR):e}, rel. diff = {numpy.fabs((tRforce-mpotderivR)/tRforce):e}" # azimuthal torque, if it exists if isinstance(tp, potential.linearPotential): continue @@ -365,8 +365,8 @@ def test_forceAsDeriv_potential(): ), f"Calculation of the vertical force as the vertical derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(mpotderivz):e}, rel. diff = {numpy.fabs((tzforce-mpotderivz)/tzforce):e}" else: assert ( - tzforce - mpotderivz - ) ** 2.0 / tzforce**2.0 < 10.0**ttol, f"Calculation of the vertical force as the vertical derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(mpotderivz):e}, rel. diff = {numpy.fabs((tzforce-mpotderivz)/tzforce):e}" + (tzforce - mpotderivz) ** 2.0 / tzforce** 2.0 < 10.0** ttol + ), f"Calculation of the vertical force as the vertical derivative of the {p} potential fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(mpotderivz):e}, rel. diff = {numpy.fabs((tzforce-mpotderivz)/tzforce):e}" # Test whether the second derivative of the potential is minus the derivative of the force @@ -542,8 +542,9 @@ def test_2ndDeriv_potential(): ), f"Calculation of the second Radial derivative of the potential as the Radial derivative of the {p} Radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tR2deriv-mRforcederivR):e}, rel. diff = {numpy.fabs((tR2deriv-mRforcederivR)/tR2deriv):e}" else: assert ( - tR2deriv - mRforcederivR - ) ** 2.0 / tR2deriv**2.0 < 10.0**ttol, f"Calculation of the second Radial derivative of the potential as the Radial derivative of the {p} Radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tR2deriv-mRforcederivR):e}, rel. diff = {numpy.fabs((tR2deriv-mRforcederivR)/tR2deriv):e}" + (tR2deriv - mRforcederivR) ** 2.0 / tR2deriv** 2.0 + < 10.0** ttol + ), f"Calculation of the second Radial derivative of the potential as the Radial derivative of the {p} Radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tR2deriv-mRforcederivR):e}, rel. diff = {numpy.fabs((tR2deriv-mRforcederivR)/tR2deriv):e}" # 2nd azimuthal if not isinstance(tp, potential.linearPotential) and hasattr(tp, "_phi2deriv"): for ii in range(len(Rs)): @@ -659,8 +660,9 @@ def test_2ndDeriv_potential(): ), f"Calculation of the second vertical derivative of the potential as the vertical derivative of the {p} vertical force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tz2deriv-mzforcederivz):e}, rel. diff = {numpy.fabs((tz2deriv-mzforcederivz)/tz2deriv):e}" else: assert ( - tz2deriv - mzforcederivz - ) ** 2.0 / tz2deriv**2.0 < 10.0**ttol, f"Calculation of the second vertical derivative of the potential as the vertical derivative of the {p} vertical force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tz2deriv-mzforcederivz):e}, rel. diff = {numpy.fabs((tz2deriv-mzforcederivz)/tz2deriv):e}" + (tz2deriv - mzforcederivz) ** 2.0 / tz2deriv** 2.0 + < 10.0** ttol + ), f"Calculation of the second vertical derivative of the potential as the vertical derivative of the {p} vertical force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tz2deriv-mzforcederivz):e}, rel. diff = {numpy.fabs((tz2deriv-mzforcederivz)/tz2deriv):e}" # mixed radial vertical if ( not isinstance(tp, potential.planarPotential) @@ -687,8 +689,9 @@ def test_2ndDeriv_potential(): ), f"Calculation of the mixed radial vertical derivative of the potential as the vertical derivative of the {p} radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRzderiv-mRforcederivz):e}, rel. diff = {numpy.fabs((tRzderiv-mRforcederivz)/tRzderiv):e}" else: assert ( - tRzderiv - mRforcederivz - ) ** 2.0 / tRzderiv**2.0 < 10.0**ttol, f"Calculation of the mixed radial vertical derivative of the potential as the vertical derivative of the {p} radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRzderiv-mRforcederivz):e}, rel. diff = {numpy.fabs((tRzderiv-mRforcederivz)/tRzderiv):e}" + (tRzderiv - mRforcederivz) ** 2.0 / tRzderiv** 2.0 + < 10.0** ttol + ), f"Calculation of the mixed radial vertical derivative of the potential as the vertical derivative of the {p} radial force fails at (R,Z) = ({Rs[ii]:.3f},{Zs[jj]:.3f}); diff = {numpy.fabs(tRzderiv-mRforcederivz):e}, rel. diff = {numpy.fabs((tRzderiv-mRforcederivz)/tRzderiv):e}" # mixed radial, azimuthal if not isinstance(tp, potential.linearPotential) and hasattr(tp, "_Rphideriv"): for ii in range(len(Rs)): @@ -719,8 +722,9 @@ def test_2ndDeriv_potential(): ), f"Calculation of the mixed radial azimuthal derivative of the potential as the azimuthal derivative of the {p} radial force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tRphideriv-mRforcederivphi):e}, rel. diff = {numpy.fabs((tRphideriv-mRforcederivphi)/tRphideriv):e}" else: assert ( - tRphideriv - mRforcederivphi - ) ** 2.0 / tRphideriv**2.0 < 10.0**ttol, f"Calculation of the mixed radial azimuthal derivative of the potential as the azimuthal derivative of the {p} radial force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tRphideriv-mRforcederivphi):e}, rel. diff = {numpy.fabs((tRphideriv-mRforcederivphi)/tRphideriv):e}" + (tRphideriv - mRforcederivphi) ** 2.0 / tRphideriv** 2.0 + < 10.0** ttol + ), f"Calculation of the mixed radial azimuthal derivative of the potential as the azimuthal derivative of the {p} radial force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tRphideriv-mRforcederivphi):e}, rel. diff = {numpy.fabs((tRphideriv-mRforcederivphi)/tRphideriv):e}" # mixed azimuthal, vertical if ( not isinstance(tp, potential.planarPotential) @@ -746,8 +750,9 @@ def test_2ndDeriv_potential(): ), f"Calculation of the mixed azimuthal vertical derivative of the potential as the azimuthal derivative of the {p} vertical force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tphizderiv-mzforcederivphi):e}, rel. diff = {numpy.fabs((tphizderiv-mzforcederivphi)/tphizderiv):e}" else: assert ( - tphizderiv - mzforcederivphi - ) ** 2.0 / tphizderiv**2.0 < 10.0**ttol, f"Calculation of the mixed azimuthal vertical derivative of the potential as the azimuthal derivative of the {p} vertical force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tphizderiv-mzforcederivphi):e}, rel. diff = {numpy.fabs((tphizderiv-mzforcederivphi)/tphizderiv):e}" + (tphizderiv - mzforcederivphi) ** 2.0 / tphizderiv** 2.0 + < 10.0** ttol + ), f"Calculation of the mixed azimuthal vertical derivative of the potential as the azimuthal derivative of the {p} vertical force fails at (R,phi) = ({Rs[ii]:.3f},{phis[jj]:.3f}); diff = {numpy.fabs(tphizderiv-mzforcederivphi):e}, rel. diff = {numpy.fabs((tphizderiv-mzforcederivphi)/tphizderiv):e}" # Test whether the Poisson equation is satisfied if _dens and the relevant second derivatives are implemented @@ -890,8 +895,8 @@ def test_poisson_potential(): ), f"Poisson equation relation between the derivatives of the potential and the implemented density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" else: assert ( - tpoissondens - tdens - ) ** 2.0 / tdens**2.0 < 10.0**ttol, f"Poisson equation relation between the derivatives of the potential and the implemented density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" + (tpoissondens - tdens) ** 2.0 / tdens** 2.0 < 10.0** ttol + ), f"Poisson equation relation between the derivatives of the potential and the implemented density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" return None @@ -993,9 +998,9 @@ def test_poisson_surfdens_potential(): tol = {} tol["default"] = -8.0 tol["DoubleExponentialDiskPotential"] = -3.0 # these are more difficult - tol["SphericalShellPotential"] = ( - -0 - ) # Direct integration fails to deal with delta function! + tol[ + "SphericalShellPotential" + ] = -0 # Direct integration fails to deal with delta function! # tol['SpiralArmsPotential']= -3 #these are more difficult # tol['rotatingSpiralArmsPotential']= -3 # tol['specialSpiralArmsPotential']= -4 @@ -1046,8 +1051,8 @@ def test_poisson_surfdens_potential(): ), f"Poisson equation relation between the derivatives of the potential and the implemented surface density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" else: assert ( - tpoissondens - tdens - ) ** 2.0 / tdens**2.0 < 10.0**ttol, f"Poisson equation relation between the derivatives of the potential and the implemented surface density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" + (tpoissondens - tdens) ** 2.0 / tdens** 2.0 < 10.0** ttol + ), f"Poisson equation relation between the derivatives of the potential and the implemented surface density is not satisfied for the {p} potential at (R,Z,phi) = ({Rs[ii]:.3f},{Zs[jj]:.3f},{phis[kk]:.3f}); diff = {numpy.fabs(tdens-tpoissondens):e}, rel. diff = {numpy.fabs((tdens-tpoissondens)/tdens):e}" if p == "mockRotatedAndTiltedTriaxialLogHaloPotentialwInclination": continue # takes a long time otherwise... skip after all z at one (R,phi) return None @@ -3067,20 +3072,14 @@ def test_dvcircdR_omegac_epifreq_rl_vesc(): # Kepler potential, vc = vc_0(R/R0)^-0.5 -> dvcdR= -0.5 vc_0 (R/R0)**-1.5 kp = potential.KeplerPotential(normalize=1.0) assert ( - kp.dvcircdR(1.0) + 0.5 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's rotation curve is not what it should be at R=1" - ) + (kp.dvcircdR(1.0) + 0.5) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's rotation curve is not what it should be at R=1" assert ( - kp.dvcircdR(0.5) + 0.5**-0.5 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's rotation curve is not what it should be at R=0.5" - ) + (kp.dvcircdR(0.5) + 0.5**-0.5) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's rotation curve is not what it should be at R=0.5" assert ( - kp.dvcircdR(2.0) + 0.5**2.5 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's rotation curve is not what it should be at R=2" - ) + (kp.dvcircdR(2.0) + 0.5**2.5) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's rotation curve is not what it should be at R=2" # Rotational frequency assert ( lp.omegac(1.0) - 1.0 @@ -3092,83 +3091,64 @@ def test_dvcircdR_omegac_epifreq_rl_vesc(): lp.omegac(2.0) - 0.5 ) ** 2.0 < 10.0**-16.0, "LogarithmicHalo's rotational frequency is off at R=2" assert ( - lp.toPlanar().omegac(2.0) - 0.5 - ) ** 2.0 < 10.0**-16.0, ( - "LogarithmicHalo's rotational frequency is off at R=2 through planarPotential" - ) + (lp.toPlanar().omegac(2.0) - 0.5) ** 2.0 < 10.0** -16.0 + ), "LogarithmicHalo's rotational frequency is off at R=2 through planarPotential" # Epicycle frequency, flat rotation curve assert ( - lp.epifreq(1.0) - numpy.sqrt(2.0) * lp.omegac(1.0) - ) ** 2.0 < 10.0**-16.0, "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=1" + (lp.epifreq(1.0) - numpy.sqrt(2.0) * lp.omegac(1.0)) ** 2.0 < 10.0** -16.0 + ), "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=1" assert ( - lp.epifreq(0.5) - numpy.sqrt(2.0) * lp.omegac(0.5) - ) ** 2.0 < 10.0**-16.0, "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=0.5" + (lp.epifreq(0.5) - numpy.sqrt(2.0) * lp.omegac(0.5)) ** 2.0 < 10.0** -16.0 + ), "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=0.5" assert ( - lp.epifreq(2.0) - numpy.sqrt(2.0) * lp.omegac(2.0) - ) ** 2.0 < 10.0**-16.0, "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=2" + (lp.epifreq(2.0) - numpy.sqrt(2.0) * lp.omegac(2.0)) ** 2.0 < 10.0** -16.0 + ), "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=2" assert ( - lp.toPlanar().epifreq(2.0) - numpy.sqrt(2.0) * lp.omegac(2.0) - ) ** 2.0 < 10.0**-16.0, "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=, through planar2" + (lp.toPlanar().epifreq(2.0) - numpy.sqrt(2.0) * lp.omegac(2.0)) ** 2.0 + < 10.0** -16.0 + ), "LogarithmicHalo's epicycle and rotational frequency are inconsistent with kappa = sqrt(2) Omega at R=, through planar2" # Epicycle frequency, Kepler assert ( - kp.epifreq(1.0) - kp.omegac(1.0) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=1" + (kp.epifreq(1.0) - kp.omegac(1.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=1" assert ( - kp.epifreq(0.5) - kp.omegac(0.5) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=0.5" + (kp.epifreq(0.5) - kp.omegac(0.5)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=0.5" assert ( - kp.epifreq(2.0) - kp.omegac(2.0) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=2" + (kp.epifreq(2.0) - kp.omegac(2.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's epicycle and rotational frequency are inconsistent with kappa = Omega at R=2" # Check radius of circular orbit, Kepler assert ( - kp.rl(1.0) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's radius of a circular orbit is wrong at Lz=1." - ) + (kp.rl(1.0) - 1.0) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's radius of a circular orbit is wrong at Lz=1." assert ( - kp.rl(0.5) - 1.0 / 4.0 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's radius of a circular orbit is wrong at Lz=0.5" - ) + (kp.rl(0.5) - 1.0 / 4.0) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's radius of a circular orbit is wrong at Lz=0.5" assert ( - kp.rl(2.0) - 4.0 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's radius of a circular orbit is wrong at Lz=2." - ) + (kp.rl(2.0) - 4.0) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's radius of a circular orbit is wrong at Lz=2." # Check radius of circular orbit, PowerSphericalPotential with close-to-flat rotation curve pp = potential.PowerSphericalPotential(alpha=1.8, normalize=1.0) assert ( - pp.rl(1.0) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=1." - ) + (pp.rl(1.0) - 1.0) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=1." assert ( - pp.rl(0.5) - 0.5 ** (10.0 / 11.0) - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=0.5" - ) + (pp.rl(0.5) - 0.5 ** (10.0 / 11.0)) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=0.5" assert ( - pp.rl(2.0) - 2.0 ** (10.0 / 11.0) - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=2." - ) + (pp.rl(2.0) - 2.0 ** (10.0 / 11.0)) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=2." # Check radius of circular orbit, PowerSphericalPotential with steeper rotation curve pp = potential.PowerSphericalPotential(alpha=0.5, normalize=1.0) assert ( - pp.rl(1.0) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=1." - ) + (pp.rl(1.0) - 1.0) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=1." assert ( - pp.rl(0.0625) - 0.0625 ** (4.0 / 7.0) - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=0.0625" - ) + (pp.rl(0.0625) - 0.0625 ** (4.0 / 7.0)) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=0.0625" assert ( - pp.rl(16.0) - 16.0 ** (4.0 / 7.0) - ) ** 2.0 < 10.0**-16.0, ( - "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=16." - ) + (pp.rl(16.0) - 16.0 ** (4.0 / 7.0)) ** 2.0 < 10.0** -16.0 + ), "PowerSphericalPotential's radius of a circular orbit is wrong at Lz=16." # Check radius in MWPotential2014 at very small lz, to test small lz behavior lz = 0.000001 assert ( @@ -3192,25 +3172,22 @@ def test_dvcircdR_omegac_epifreq_rl_vesc(): kp.vesc(2.0) ** 2.0 - 2.0 * kp.vcirc(2.0) ** 2.0 ) ** 2.0 < 10.0**-16.0, "KeplerPotential's escape velocity is wrong at R=2" assert ( - kp.toPlanar().vesc(2.0) ** 2.0 - 2.0 * kp.vcirc(2.0) ** 2.0 - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's escape velocity is wrong at R=2, through planar" - ) + (kp.toPlanar().vesc(2.0) ** 2.0 - 2.0 * kp.vcirc(2.0) ** 2.0) ** 2.0 + < 10.0** -16.0 + ), "KeplerPotential's escape velocity is wrong at R=2, through planar" # W/ different interface assert ( - kp.vcirc(1.0) - potential.vcirc(kp, 1.0) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's circular velocity does not agree between kp.vcirc and vcirc(kp)" + (kp.vcirc(1.0) - potential.vcirc(kp, 1.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's circular velocity does not agree between kp.vcirc and vcirc(kp)" assert ( - kp.vcirc(1.0) - potential.vcirc(kp.toPlanar(), 1.0) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's circular velocity does not agree between kp.vcirc and vcirc(kp.toPlanar)" + (kp.vcirc(1.0) - potential.vcirc(kp.toPlanar(), 1.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's circular velocity does not agree between kp.vcirc and vcirc(kp.toPlanar)" assert ( - kp.vesc(1.0) - potential.vesc(kp, 1.0) - ) ** 2.0 < 10.0**-16.0, ( - "KeplerPotential's escape velocity does not agree between kp.vesc and vesc(kp)" - ) + (kp.vesc(1.0) - potential.vesc(kp, 1.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's escape velocity does not agree between kp.vesc and vesc(kp)" assert ( - kp.vesc(1.0) - potential.vesc(kp.toPlanar(), 1.0) - ) ** 2.0 < 10.0**-16.0, "KeplerPotential's escape velocity does not agree between kp.vesc and vesc(kp.toPlanar)" + (kp.vesc(1.0) - potential.vesc(kp.toPlanar(), 1.0)) ** 2.0 < 10.0** -16.0 + ), "KeplerPotential's escape velocity does not agree between kp.vesc and vesc(kp.toPlanar)" return None @@ -3504,22 +3481,16 @@ def test_flattening(): # Check some spherical potentials kp = potential.KeplerPotential(normalize=1.0) assert ( - kp.flattening(1.0, 0.02) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "Flattening of KeplerPotential is not equal to 1 at (R,z) = (1.,0.02)" - ) + (kp.flattening(1.0, 0.02) - 1.0) ** 2.0 < 10.0** -16.0 + ), "Flattening of KeplerPotential is not equal to 1 at (R,z) = (1.,0.02)" np = potential.NFWPotential(normalize=1.0, a=5.0) assert ( - np.flattening(1.0, 0.02) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "Flattening of NFWPotential is not equal to 1 at (R,z) = (1.,0.02)" - ) + (np.flattening(1.0, 0.02) - 1.0) ** 2.0 < 10.0** -16.0 + ), "Flattening of NFWPotential is not equal to 1 at (R,z) = (1.,0.02)" hp = potential.HernquistPotential(normalize=1.0, a=5.0) assert ( - hp.flattening(1.0, 0.02) - 1.0 - ) ** 2.0 < 10.0**-16.0, ( - "Flattening of HernquistPotential is not equal to 1 at (R,z) = (1.,0.02)" - ) + (hp.flattening(1.0, 0.02) - 1.0) ** 2.0 < 10.0** -16.0 + ), "Flattening of HernquistPotential is not equal to 1 at (R,z) = (1.,0.02)" # Disk potentials should be oblate everywhere mp = potential.MiyamotoNagaiPotential(normalize=1.0, a=0.5, b=0.05) assert ( @@ -4559,9 +4530,7 @@ def test_NFW_virialquantities_diffrovo(): * 3.0 / np.rvir(ro=ro, vo=vo, H=H, Om=Om, overdens=overdens, wrtcrit=wrtcrit) ** 3.0 - ) * ( - 10.0**6.0 / H**2.0 * 8.0 * numpy.pi / 3.0 / Om * (4.302 * 10.0**-6.0) - ) + ) * (10.0**6.0 / H**2.0 * 8.0 * numpy.pi / 3.0 / Om * (4.302 * 10.0**-6.0)) assert ( numpy.fabs(od - overdens) < 0.1 ), "NFWPotential's virial quantities computed in physical units with different (ro,vo) from setup are incorrect" @@ -4966,9 +4935,7 @@ def test_MN3ExponentialDiskPotential_inputs(): raisedWarning = "MN3ExponentialDiskPotential" in str(wa.message) if raisedWarning: break - assert ( - raisedWarning - ), "MN3ExponentialDiskPotential w/o posdens, but with b/Rd > 3 did not raise galpyWarning" + assert raisedWarning, "MN3ExponentialDiskPotential w/o posdens, but with b/Rd > 3 did not raise galpyWarning" with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) mn = MN3ExponentialDiskPotential(normalize=1.0, hr=1.0, hz=0.7727, posdens=True) @@ -4977,9 +4944,7 @@ def test_MN3ExponentialDiskPotential_inputs(): raisedWarning = "MN3ExponentialDiskPotential" in str(wa.message) if raisedWarning: break - assert ( - raisedWarning - ), "MN3ExponentialDiskPotential w/o posdens, but with b/Rd > 1.35 did not raise galpyWarning" + assert raisedWarning, "MN3ExponentialDiskPotential w/o posdens, but with b/Rd > 1.35 did not raise galpyWarning" return None @@ -5363,12 +5328,14 @@ def test_SoftenedNeedleBarPotential_density(): sbp = potential.SoftenedNeedleBarPotential( normalize=1.0, a=1.0, c=0.1, b=0.3, pa=0.0 ) - assert sbp.dens(2.0, 0.0, phi=numpy.pi / 4.0) > sbp.dens( - numpy.sqrt(2.0), numpy.sqrt(2.0), phi=0.0 + assert ( + sbp.dens(2.0, 0.0, phi=numpy.pi / 4.0) + > sbp.dens(numpy.sqrt(2.0), numpy.sqrt(2.0), phi=0.0) ), "SoftenedNeedleBarPotential with flattened softening kernel does not appear to have a consistent" # Another one - assert sbp.dens(4.0, 0.0, phi=numpy.pi / 4.0) > sbp.dens( - 2.0 * numpy.sqrt(2.0), 2.0 * numpy.sqrt(2.0), phi=0.0 + assert ( + sbp.dens(4.0, 0.0, phi=numpy.pi / 4.0) + > sbp.dens(2.0 * numpy.sqrt(2.0), 2.0 * numpy.sqrt(2.0), phi=0.0) ), "SoftenedNeedleBarPotential with flattened softening kernel does not appear to have a consistent" return None @@ -6184,14 +6151,17 @@ def test_RotateAndTiltWrapper(): ) # test against the transformed potential and a MWPotential evaluated at the transformed coords assert ( - evaluatePotentials(testpot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - - evaluatePotentials( - potential.MWPotential2014, - tRphiz_test[0], - tRphiz_test[2], - phi=tRphiz_test[1], + ( + evaluatePotentials(testpot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) + - evaluatePotentials( + potential.MWPotential2014, + tRphiz_test[0], + tRphiz_test[2], + phi=tRphiz_test[1], + ) ) - ) < 1e-6, "Evaluating potential at same relative position in a Rotated and tilted MWPotential2014 and non-Rotated does not give same result" + < 1e-6 + ), "Evaluating potential at same relative position in a Rotated and tilted MWPotential2014 and non-Rotated does not give same result" # Also a triaxial NFW NFW_wrapped = potential.RotateAndTiltWrapperPotential( zvec=zvec, @@ -6202,18 +6172,32 @@ def test_RotateAndTiltWrapper(): amp=1.0, zvec=zvec, pa=galaxy_pa, b=0.7, c=0.5 ) assert ( - evaluatePotentials(NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - - evaluatePotentials(NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - ) < 1e-6, "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" + ( + evaluatePotentials( + NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + - evaluatePotentials( + NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + ) + < 1e-6 + ), "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" # Try not specifying galaxy_pa, shouldn be =0 NFW_wrapped = potential.RotateAndTiltWrapperPotential( zvec=zvec, pot=potential.TriaxialNFWPotential(amp=1.0, b=0.7, c=0.5) ) NFW_rot = potential.TriaxialNFWPotential(amp=1.0, zvec=zvec, pa=0.0, b=0.7, c=0.5) assert ( - evaluatePotentials(NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - - evaluatePotentials(NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - ) < 1e-6, "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" + ( + evaluatePotentials( + NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + - evaluatePotentials( + NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + ) + < 1e-6 + ), "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" # Try not specifying zvec, should be =[0,0,1] NFW_wrapped = potential.RotateAndTiltWrapperPotential( galaxy_pa=galaxy_pa, pot=potential.TriaxialNFWPotential(amp=1.0, b=0.7, c=0.5) @@ -6222,9 +6206,16 @@ def test_RotateAndTiltWrapper(): amp=1.0, zvec=[0.0, 0.0, 1.0], pa=galaxy_pa, b=0.7, c=0.5 ) assert ( - evaluatePotentials(NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - - evaluatePotentials(NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1]) - ) < 1e-6, "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" + ( + evaluatePotentials( + NFW_wrapped, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + - evaluatePotentials( + NFW_rot, Rphiz_test[0], Rphiz_test[2], phi=Rphiz_test[1] + ) + ) + < 1e-6 + ), "Wrapped and Internally rotated NFW potentials do not match when evaluated at the same point" # make sure the offset works as intended # triaxial NFW at x,y,z = [20.,0.,3.] NFW_wrapped = potential.RotateAndTiltWrapperPotential( @@ -6237,9 +6228,12 @@ def test_RotateAndTiltWrapper(): amp=1.0, zvec=zvec, pa=galaxy_pa, b=0.7, c=0.5 ) assert ( - evaluatePotentials(NFW_wrapped, 0.0, 0.0, phi=0.0) - - evaluatePotentials(NFW_rot, 20.0, -3.0, phi=numpy.pi) - ) < 1e-6, "Wrapped + Offset and Internally rotated NFW potentials do not match when evaluated at the same point" + ( + evaluatePotentials(NFW_wrapped, 0.0, 0.0, phi=0.0) + - evaluatePotentials(NFW_rot, 20.0, -3.0, phi=numpy.pi) + ) + < 1e-6 + ), "Wrapped + Offset and Internally rotated NFW potentials do not match when evaluated at the same point" def test_integration_RotateAndTiltWrapper(): diff --git a/tests/test_quantity.py b/tests/test_quantity.py index a179a502e..c8ed7abec 100644 --- a/tests/test_quantity.py +++ b/tests/test_quantity.py @@ -16725,23 +16725,24 @@ def test_streamdf_method_inputAsQuantity(): tdisrupt=4.5 / conversion.time_in_Gyr(220.0, 8.0), nosetup=True, ) - assert numpy.fabs( - sdf_bovy14.subhalo_encounters( - venc=200.0 * units.km / units.s, - sigma=150.0 * units.km / units.s, - nsubhalo=38.35 / (4.0 * (25.0 * units.kpc) ** 3.0 * numpy.pi / 3.0), - bmax=1.0 * units.kpc, - yoon=False, - ) - - sdf_bovy14_nou.subhalo_encounters( - venc=200.0 / vo, - sigma=150.0 / vo, - nsubhalo=38.35 / (4.0 * 25.0**3.0 * numpy.pi / 3.0) * ro**3.0, - bmax=1.0 / ro, - yoon=False, - ) - ) < 1e-6 * _NUMPY_1_22 + 1e-8 * ( - 1 - _NUMPY_1_22 + assert ( + numpy.fabs( + sdf_bovy14.subhalo_encounters( + venc=200.0 * units.km / units.s, + sigma=150.0 * units.km / units.s, + nsubhalo=38.35 / (4.0 * (25.0 * units.kpc) ** 3.0 * numpy.pi / 3.0), + bmax=1.0 * units.kpc, + yoon=False, + ) + - sdf_bovy14_nou.subhalo_encounters( + venc=200.0 / vo, + sigma=150.0 / vo, + nsubhalo=38.35 / (4.0 * 25.0**3.0 * numpy.pi / 3.0) * ro**3.0, + bmax=1.0 / ro, + yoon=False, + ) + ) + < 1e-6 * _NUMPY_1_22 + 1e-8 * (1 - _NUMPY_1_22) ), "streamdf method subhalo_encounters with Quantity input does not return correct Quantity" assert numpy.fabs( sdf_bovy14.pOparapar(0.2 / units.Gyr, 30.0 * units.deg) @@ -17800,12 +17801,8 @@ def test_orbitmethodswunits_quantity_issue326(): o = Orbit([1.0, 0.1, 1.1, 0.1, 0.2, 0.0]) # First make sure we're testing what we want to test - assert ( - not o._roSet - ), "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _roSet is True" - assert ( - not o._voSet - ), "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _voSet is True" + assert not o._roSet, "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _roSet is True" + assert not o._voSet, "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _voSet is True" # Then test methods assert isinstance( o.ra(), units.Quantity @@ -17867,12 +17864,8 @@ def test_orbitmethodswunits_quantity_overrideusephysical_issue326(): o = Orbit([1.0, 0.1, 1.1, 0.1, 0.2, 0.0]) # First make sure we're testing what we want to test - assert ( - not o._roSet - ), "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _roSet is True" - assert ( - not o._voSet - ), "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _voSet is True" + assert not o._roSet, "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _roSet is True" + assert not o._voSet, "Test of whether or not Orbit methods that should always return a Quantity do so cannot run meaningfully when _voSet is True" # Then test methods assert isinstance( o.ra(use_physical=False), units.Quantity diff --git a/tests/test_scf.py b/tests/test_scf.py index 30ef5069b..94c64632d 100644 --- a/tests/test_scf.py +++ b/tests/test_scf.py @@ -611,9 +611,7 @@ def test_FutureWarning_multid_indexing(): ) if raisedWarning: break - assert ( - not raisedWarning - ), "SCFPotential should not raise 'FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated ...', but did" + assert not raisedWarning, "SCFPotential should not raise 'FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated ...', but did" return None diff --git a/tests/test_snapshotpotential.py b/tests/test_snapshotpotential.py index 4faf3b5db..07128e261 100644 --- a/tests/test_snapshotpotential.py +++ b/tests/test_snapshotpotential.py @@ -686,19 +686,17 @@ def test_interpsnapshotKeplerpotential_Rzderiv(): # Test, but give small |z| a less constraining for r in rs: for z in zs: - assert numpy.fabs( - (sp.Rzderiv(r, z) - kp.Rzderiv(r, z)) / kp.Rzderiv(r, z) - ) < 10.0**-4.0 * ( - 1.0 + 19.0 * (numpy.fabs(z) < 0.05) + assert ( + numpy.fabs((sp.Rzderiv(r, z) - kp.Rzderiv(r, z)) / kp.Rzderiv(r, z)) + < 10.0** -4.0 * (1.0 + 19.0 * (numpy.fabs(z) < 0.05)) ), f"RZPot interpolation of Rzderiv w/ InterpSnapShotPotential of KeplerPotential fails at (R,z) = ({r:g},{z:g}) by {numpy.fabs((sp.Rzderiv(r,z)-kp.Rzderiv(r,z))/kp.Rzderiv(r,z)):g}; value is {kp.Rzderiv(r,z):g}" # This tests within the grid rs = numpy.linspace(0.01, 2.0, 10)[1:] zs = numpy.linspace(-0.2, 0.2, 20) for r in rs: for z in zs: - assert numpy.fabs( - (sp.Rzderiv(r, z) - kp.Rzderiv(r, z)) / kp.Rzderiv(r, z) - ) < 10.0**-4.0 * ( - 1.0 + 19.0 * (numpy.fabs(z) < 0.05) + assert ( + numpy.fabs((sp.Rzderiv(r, z) - kp.Rzderiv(r, z)) / kp.Rzderiv(r, z)) + < 10.0** -4.0 * (1.0 + 19.0 * (numpy.fabs(z) < 0.05)) ), f"RZPot interpolation of Rzderiv w/ InterpSnapShotPotential of KeplerPotential fails at (R,z) = ({r:g},{z:g}) by {numpy.fabs((sp.Rzderiv(r,z)-kp.Rzderiv(r,z))/kp.Rzderiv(r,z)):g}" return None diff --git a/tests/test_sphericaldf.py b/tests/test_sphericaldf.py index a800ba6f9..b33f6fb2d 100644 --- a/tests/test_sphericaldf.py +++ b/tests/test_sphericaldf.py @@ -2634,9 +2634,7 @@ def test_eddington_interpolatedpotentials_rmin(): str(rec.message.args[0]) == "Interpolated potential grid rmin is larger than the rmin to be used for the v_vesc_interpolator grid. This may adversely affect the generated samples. Proceed with care!" ) - assert ( - raisedWarning - ), "Using an interpolated potential with rmin smaller than the rmin to be used for the v_vesc_interpolator grid should have raised a warning, but didn't" + assert raisedWarning, "Using an interpolated potential with rmin smaller than the rmin to be used for the v_vesc_interpolator grid should have raised a warning, but didn't" def test_eddington_interpolatedpotentials_rmax(): @@ -2654,9 +2652,7 @@ def test_eddington_interpolatedpotentials_rmax(): str(rec.message.args[0]) == "The interpolated potential's rmax is smaller than the DF's rmax" ) - assert ( - raisedWarning - ), "Using an interpolated potential with rmax smaller than the DF's rmax should have raised a warning, but didn't" + assert raisedWarning, "Using an interpolated potential with rmax smaller than the DF's rmax should have raised a warning, but didn't" ########################### TESTS OF ERRORS AND WARNINGS####################### @@ -2713,9 +2709,7 @@ def test_anisotropic_hernquist_negdf(): str(rec.message.args[0]) == "The DF appears to have negative regions; we'll try to ignore these for sampling the DF, but this may adversely affect the generated samples. Proceed with care!" ) - assert ( - raisedWarning - ), "Using an anisotropic Hernquist DF that has negative parts should have raised a warning, but didn't" + assert raisedWarning, "Using an anisotropic Hernquist DF that has negative parts should have raised a warning, but didn't" ############################# TESTS OF UNIT HANDLING########################### diff --git a/tests/test_streamdf.py b/tests/test_streamdf.py index a3cde9073..6faa86af5 100644 --- a/tests/test_streamdf.py +++ b/tests/test_streamdf.py @@ -143,9 +143,7 @@ def test_progenitor_coordtransformparams(): ) if raisedWarning: break - assert ( - raisedWarning - ), "streamdf setup does not raise warning when progenitor's ro is different from ro" + assert raisedWarning, "streamdf setup does not raise warning when progenitor's ro is different from ro" # Test w/ diff R0 with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -169,9 +167,7 @@ def test_progenitor_coordtransformparams(): ) if raisedWarning: break - assert ( - raisedWarning - ), "streamdf setup does not raise warning when progenitor's ro is different from R0" + assert raisedWarning, "streamdf setup does not raise warning when progenitor's ro is different from R0" # Test w/ diff Vnorm with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -197,9 +193,7 @@ def test_progenitor_coordtransformparams(): ) if raisedWarning: break - assert ( - raisedWarning - ), "streamdf setup does not raise warning when progenitor's vo is different from vo" + assert raisedWarning, "streamdf setup does not raise warning when progenitor's vo is different from vo" # Test w/ diff zo with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -226,9 +220,7 @@ def test_progenitor_coordtransformparams(): ) if raisedWarning: break - assert ( - raisedWarning - ), "streamdf setup does not raise warning when progenitor's zo is different from Zsun" + assert raisedWarning, "streamdf setup does not raise warning when progenitor's zo is different from Zsun" # Test w/ diff vsun with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", galpyWarning) @@ -256,9 +248,7 @@ def test_progenitor_coordtransformparams(): ) if raisedWarning: break - assert ( - raisedWarning - ), "streamdf setup does not raise warning when progenitor's solarmotion is different from vsun" + assert raisedWarning, "streamdf setup does not raise warning when progenitor's solarmotion is different from vsun" return None @@ -276,13 +266,11 @@ def test_bovy14_freqratio(bovy14_setup): sdf_bovy14 = bovy14_setup # Test the frequency ratio assert ( - sdf_bovy14.freqEigvalRatio() - 30.0 - ) ** 2.0 < 10.0**0.0, ( - "streamdf model from Bovy (2014) does not give a frequency ratio of about 30" - ) + (sdf_bovy14.freqEigvalRatio() - 30.0) ** 2.0 < 10.0** 0.0 + ), "streamdf model from Bovy (2014) does not give a frequency ratio of about 30" assert ( - sdf_bovy14.freqEigvalRatio(isotropic=True) - 34.0 - ) ** 2.0 < 10.0**0.0, "streamdf model from Bovy (2014) does not give an isotropic frequency ratio of about 34" + (sdf_bovy14.freqEigvalRatio(isotropic=True) - 34.0) ** 2.0 < 10.0** 0.0 + ), "streamdf model from Bovy (2014) does not give an isotropic frequency ratio of about 34" return None @@ -291,11 +279,12 @@ def test_bovy14_misalignment(bovy14_setup): sdf_bovy14 = bovy14_setup # Test the misalignment assert ( - sdf_bovy14.misalignment() / numpy.pi * 180.0 + 0.5 - ) ** 2.0 < 10.0**-2.0, "streamdf model from Bovy (2014) does not give a misalighment of about -0.5 degree" + (sdf_bovy14.misalignment() / numpy.pi * 180.0 + 0.5) ** 2.0 < 10.0** -2.0 + ), "streamdf model from Bovy (2014) does not give a misalighment of about -0.5 degree" assert ( - sdf_bovy14.misalignment(isotropic=True) / numpy.pi * 180.0 - 1.3 - ) ** 2.0 < 10.0**-2.0, "streamdf model from Bovy (2014) does not give an isotropic misalighment of about 1.3 degree" + (sdf_bovy14.misalignment(isotropic=True) / numpy.pi * 180.0 - 1.3) ** 2.0 + < 10.0** -2.0 + ), "streamdf model from Bovy (2014) does not give an isotropic misalighment of about 1.3 degree" return None @@ -2810,11 +2799,12 @@ def check_closest_trackpoint(sdf, trackp, usev=False, xy=True, interp=True): xy=xy, usev=usev, ) - assert ( - indx == trackp - ), "Closest trackpoint to close to a trackpoint is not that trackpoint (%i,%i)" % ( - indx, - trackp, + assert indx == trackp, ( + "Closest trackpoint to close to a trackpoint is not that trackpoint (%i,%i)" + % ( + indx, + trackp, + ) ) return None diff --git a/tests/test_streamspraydf.py b/tests/test_streamspraydf.py index c3028463d..59abc4af3 100644 --- a/tests/test_streamspraydf.py +++ b/tests/test_streamspraydf.py @@ -342,12 +342,8 @@ def test_sample_orbit_rovoetc(): tdisrupt=4.5 / conversion.time_in_Gyr(vo, ro), ) sam = spdf_bovy14.sample(n=10) - assert ( - obs._roSet - ), "Test requires that ro be set for the progenitor orbit, but it appears not to have been set" - assert ( - not obs._voSet - ), "Test requires that vo not be set for the progenitor orbit, but it appears to have been set" + assert obs._roSet, "Test requires that ro be set for the progenitor orbit, but it appears not to have been set" + assert not obs._voSet, "Test requires that vo not be set for the progenitor orbit, but it appears to have been set" assert ( obs._roSet is sam._roSet ), "Sampled streamspraydf orbits do not have the same roSet as the progenitor orbit" @@ -383,12 +379,8 @@ def test_sample_orbit_rovoetc(): tdisrupt=4.5 / conversion.time_in_Gyr(vo, ro), ) sam = spdf_bovy14.sample(n=10) - assert ( - obs._voSet - ), "Test requires that vo be set for the progenitor orbit, but it appears not to have been set" - assert ( - not obs._roSet - ), "Test requires that ro not be set for the progenitor orbit, but it appears to have been set" + assert obs._voSet, "Test requires that vo be set for the progenitor orbit, but it appears not to have been set" + assert not obs._roSet, "Test requires that ro not be set for the progenitor orbit, but it appears to have been set" assert ( obs._roSet is sam._roSet ), "Sampled streamspraydf orbits do not have the same roSet as the progenitor orbit"