diff --git a/doc/source/reference/bovycoords.rst b/doc/source/reference/bovycoords.rst index fedb5aaaf..6016ceacd 100644 --- a/doc/source/reference/bovycoords.rst +++ b/doc/source/reference/bovycoords.rst @@ -16,6 +16,8 @@ self-explanatory names: custom_to_radec cyl_to_rect cyl_to_rect_vec + cyl_to_spher + cyl_to_spher_vec dl_to_rphi_2d galcencyl_to_XYZ galcencyl_to_vxvyvz @@ -36,6 +38,8 @@ self-explanatory names: Rz_to_coshucosv Rz_to_uv sphergal_to_rectgal + spher_to_cyl + spher_to_cyl_vec uv_to_Rz vrpmllpmbb_to_vxvyvz vRvz_to_pupv diff --git a/doc/source/reference/coordscyltospher.rst b/doc/source/reference/coordscyltospher.rst new file mode 100644 index 000000000..2de092e52 --- /dev/null +++ b/doc/source/reference/coordscyltospher.rst @@ -0,0 +1,4 @@ +galpy.util.bovy_coords.cyl_to_spher +==================================== + +.. autofunction:: galpy.util.bovy_coords.cyl_to_spher \ No newline at end of file diff --git a/doc/source/reference/coordscyltosphervec.rst b/doc/source/reference/coordscyltosphervec.rst new file mode 100644 index 000000000..31adc8a95 --- /dev/null +++ b/doc/source/reference/coordscyltosphervec.rst @@ -0,0 +1,4 @@ +galpy.util.bovy_coords.cyl_to_spher_vec +======================================== + +.. autofunction:: galpy.util.bovy_coords.cyl_to_spher_vec \ No newline at end of file diff --git a/doc/source/reference/coordssphertocyl.rst b/doc/source/reference/coordssphertocyl.rst new file mode 100644 index 000000000..40fd013db --- /dev/null +++ b/doc/source/reference/coordssphertocyl.rst @@ -0,0 +1,4 @@ +galpy.util.bovy_coords.spher_to_cyl +==================================== + +.. autofunction:: galpy.util.bovy_coords.spher_to_cyl \ No newline at end of file diff --git a/doc/source/reference/coordssphertocylvec.rst b/doc/source/reference/coordssphertocylvec.rst new file mode 100644 index 000000000..463525d59 --- /dev/null +++ b/doc/source/reference/coordssphertocylvec.rst @@ -0,0 +1,4 @@ +galpy.util.bovy_coords.spher_to_cyl_vec +======================================== + +.. autofunction:: galpy.util.bovy_coords.spher_to_cyl_vec \ No newline at end of file diff --git a/doc/source/reference/orbit.rst b/doc/source/reference/orbit.rst index e807c111d..ebd11f872 100644 --- a/doc/source/reference/orbit.rst +++ b/doc/source/reference/orbit.rst @@ -89,6 +89,7 @@ Methods rguiding rperi SkyCoord + theta time toLinear toPlanar @@ -105,8 +106,10 @@ Methods vll vlos vphi + vr vR vra + vtheta vT vx vy diff --git a/doc/source/reference/orbitsphvr.rst b/doc/source/reference/orbitsphvr.rst new file mode 100644 index 000000000..09664f5a5 --- /dev/null +++ b/doc/source/reference/orbitsphvr.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.vr +===================== + +.. automethod:: galpy.orbit.Orbit.vr \ No newline at end of file diff --git a/doc/source/reference/orbittheta.rst b/doc/source/reference/orbittheta.rst new file mode 100644 index 000000000..e5fd0374c --- /dev/null +++ b/doc/source/reference/orbittheta.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.theta +======================== + +.. automethod:: galpy.orbit.Orbit.theta \ No newline at end of file diff --git a/doc/source/reference/orbitvtheta.rst b/doc/source/reference/orbitvtheta.rst new file mode 100644 index 000000000..66d9d17ca --- /dev/null +++ b/doc/source/reference/orbitvtheta.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.vtheta +========================= + +.. automethod:: galpy.orbit.Orbit.vtheta \ No newline at end of file diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 7cdecd3b5..bf01e0fc9 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -3275,7 +3275,111 @@ def vphi(self,*args,**kwargs): """ thiso= self._call_internal(*args,**kwargs) return (thiso[2]/thiso[0]).T + + @physical_conversion('velocity') + @shapeDecorator + def vr(self,*args,**kwargs): + """ + NAME: + + vr + + PURPOSE: + + return spherical radial velocity. For < 3 dimensions returns vR + + INPUT: + + t - (optional) time at which to get the radial velocity + + vo= (Object-wide default) physical scale for velocities to use to convert + + use_physical= use to override Object-wide default for using a physical scale for output + + OUTPUT: + + vr(t) [*input_shape,nt] + + HISTORY: + + 2020-07-01 - Written - James (UofT) + + """ + thiso = self._call_internal(*args,**kwargs) + if self.dim() == 3: + r = numpy.sqrt( numpy.square(thiso[0]) + numpy.square(thiso[3]) ) + return ((thiso[0]*thiso[1]+thiso[3]*thiso[4])/r).T + else: + return thiso[1].T + + @physical_conversion('velocity') + @shapeDecorator + def vtheta(self,*args,**kwargs): + """ + NAME: + + vtheta + + PURPOSE: + + return spherical polar velocity + INPUT: + + t - (optional) time at which to get the theta velocity + + vo= (Object-wide default) physical scale for velocities to use to convert + + use_physical= use to override Object-wide default for using a physical scale for output + + OUTPUT: + + vtheta(t) [*input_shape,nt] + + HISTORY: + + 2020-07-01 - Written - James (UofT) + + """ + thiso = self._call_internal(*args,**kwargs) + if not self.dim() == 3: + raise AttributeError("Orbit must be 3D to use vtheta()") + else: + r = numpy.sqrt(numpy.square(thiso[0])+numpy.square(thiso[3])) + return ((thiso[1]*thiso[3]-thiso[0]*thiso[4])/r).T + + @physical_conversion('angle') + @shapeDecorator + def theta(self,*args,**kwargs): + """ + NAME: + + theta + + PURPOSE: + + return spherical polar angle + + INPUT: + + t - (optional) time at which to get the angle + + OUTPUT: + + theta(t) [*input_shape,nt] + + HISTORY: + + 2020-07-01 - Written - James (UofT) + + """ + thiso = self._call_internal(*args,**kwargs) + if self.dim() != 3: + raise AttributeError("Orbit must be 3D to use theta()") + else: + return numpy.arctan2(thiso[0],thiso[3]) + + @physical_conversion('angle_deg') @shapeDecorator def ra(self,*args,**kwargs): diff --git a/galpy/util/bovy_coords.py b/galpy/util/bovy_coords.py index a081b73a8..377997225 100644 --- a/galpy/util/bovy_coords.py +++ b/galpy/util/bovy_coords.py @@ -1397,6 +1397,75 @@ def galcencyl_to_vxvyvz(vR,vT,vZ,phi,vsun=[0.,1.,0.],Xsun=1.,Zsun=0., return galcenrect_to_vxvyvz(vXg,vYg,vZg,vsun=vsun,Xsun=Xsun,Zsun=Zsun, _extra_rot=_extra_rot) +def cyl_to_spher_vec(vR,vT,vz,R,z): + """ + NAME: + + cyl_to_spher_vec + + PURPOSE: + + transform vectors from cylindrical to spherical coordinates. vtheta is positive from pole towards equator. + + INPUT: + + vR - Galactocentric cylindrical radial velocity + + vT - Galactocentric cylindrical tangential velocity + + vz - Galactocentric cylindrical vertical velocity + + R - Galactocentric cylindrical radius + + z - Galactocentric cylindrical height + + OUTPUT: + + vr,vT,vtheta + + HISTORY: + + 2020-07-01 - Written - James Lane (UofT) + + """ + r = numpy.sqrt(numpy.square(R) + numpy.square(z)) + vr = (R*vR + z*vz)/r + vtheta = (z*vR - R*vz)/r + return (vr,vT,vtheta) + +def spher_to_cyl_vec(vr,vT,vtheta,theta): + """ + NAME: + + spher_to_cyl_vec + + PURPOSE: + + transform vectors from spherical polar to cylindrical coordinates. vtheta is positive from pole towards equator, theta is 0 at pole + + INPUT: + + vr - Galactocentric spherical radial velocity + + vT - Galactocentric spherical azimuthal velocity + + vtheta - Galactocentric spherical polar velocity + + theta - Galactocentric spherical polar angle + + OUTPUT: + + vR,vT,vz + + HISTORY: + + 2020-07-01 - Written - James Lane (UofT) + + """ + vR = vr*numpy.sin(theta) + vtheta*numpy.cos(theta) + vz = vr*numpy.cos(theta) - vtheta*numpy.sin(theta) + return (vR,vT,vz) + def rect_to_cyl_vec(vx,vy,vz,X,Y,Z,cyl=False): """ NAME: diff --git a/tests/test_coords.py b/tests/test_coords.py index eee250759..491c80497 100644 --- a/tests/test_coords.py +++ b/tests/test_coords.py @@ -1412,6 +1412,38 @@ def test_lbd_to_XYZ_jac(): assert numpy.fabs(jac[5,5]-numpy.sqrt(3.)/2.*d*4.740470463496208) < 10.**-10., 'lbd_to_XYZ_jac calculation did not work as expected' return None +def test_cyl_to_spher_vec(): + # Test 45 degrees, disk plane, & polar location + vr,vT,vtheta = bovy_coords.cyl_to_spher_vec(0.6,1.3,0.6,1.,1.) + assert numpy.fabs(vr-0.6*2**0.5) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vtheta-0) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vT-1.3) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + vr,vT,vtheta = bovy_coords.cyl_to_spher_vec(-1.2,-0.7,-0.8,1.,0.) + assert numpy.fabs(vr+1.2) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vtheta-0.8) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vT+0.7) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + vr,vT,vtheta = bovy_coords.cyl_to_spher_vec(-1.2,-0.7,-0.8,0.,1.) + assert numpy.fabs(vr+0.8) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vtheta+1.2) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + assert numpy.fabs(vT+0.7) < 10.**-8, 'cyl_to_spher_vec does not work as expected' + return None + +def test_spher_to_cyl_vec(): + # Test 45 degrees, disk plane, & polar location + vR,vT,vz = bovy_coords.spher_to_cyl_vec(0.7,1.4,0.7,numpy.pi/4.) + assert numpy.fabs(vR-0.7*2**0.5) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vT-1.4) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vz-0.) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + vR,vT,vz = bovy_coords.spher_to_cyl_vec(0.5,-1.3,0.7,0.) + assert numpy.fabs(vR-0.7) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vT+1.3) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vz-0.5) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + vR,vT,vz = bovy_coords.spher_to_cyl_vec(0.5,-1.3,0.7,numpy.pi/2.) + assert numpy.fabs(vR-0.5) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vT+1.3) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + assert numpy.fabs(vz+0.7) < 10.**-8, 'spher_to_cyl_vec does not work as expected' + return None + def test_cyl_to_spher(): # Just a few quick tests r,t,p= bovy_coords.cyl_to_spher(1.2,3.2,1.) diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 9e4177804..69e96ba96 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -2192,6 +2192,11 @@ def test_orbit_setup(): assert numpy.fabs(o.pmll(obs=obs)-0.5) < 10.**-10., 'Orbit pmll setup does not agree with o.pmll()' assert numpy.fabs(o.pmbb(obs=obs)-0.) < 10.**-10., 'Orbit pmbb setup does not agree with o.pmbb()' assert numpy.fabs(o.vlos(obs=obs)-30.) < 10.**-10., 'Orbit vlos setup does not agree with o.vlos()' + #test galactocentric spherical coordinates + o= Orbit([2.,2.**0.5,1.,2.,2.**0.5,0.5]) + assert numpy.fabs(o.vr()-2.) < 10.**-10., 'Orbit galactocentric spherical coordinates are not correct' + assert numpy.fabs(o.vtheta()-0.) < 10.**-10., 'Orbit galactocentric spherical coordinates are not correct' + assert numpy.fabs(o.theta()-numpy.pi/4) < 10.**-10., 'Orbit galactocentric spherical coordinates are not correct' return None def test_orbit_setup_SkyCoord(): @@ -2409,6 +2414,18 @@ def test_attributeerrors(): pass else: raise AssertionError("o.vz() for planarOrbit should have raised AttributeError, but did not") + try: + o.theta() + except AttributeError: + pass + else: + raise AssertionError("o.theta() for planarROrbit should have raise AttributeError, but did not") + try: + o.vtheta() + except AttributeError: + pass + else: + raise AssertionError("o.vtheta() for planarROrbit should have raise AttributeError, but did not") #phi, x, y, vx, vy for Orbits that don't track phi o= Orbit([1.,0.1,1.1,0.1,0.2]) try: @@ -2472,6 +2489,18 @@ def test_attributeerrors(): pass else: raise AssertionError("o.vy() for planarROrbit should have raised AttributeError, but did not") + try: + o.theta() + except AttributeError: + pass + else: + raise AssertionError("o.theta() for planarROrbit should have raise AttributeError, but did not") + try: + o.vtheta() + except AttributeError: + pass + else: + raise AssertionError("o.vtheta() for planarROrbit should have raise AttributeError, but did not") return None # Test reversing an orbit @@ -3300,6 +3329,9 @@ def test_scalar_all(): assert isinstance(o.y(),float), 'Orbit.y() does not return a scalar' assert isinstance(o.vx(),float), 'Orbit.vx() does not return a scalar' assert isinstance(o.vy(),float), 'Orbit.vy() does not return a scalar' + assert isinstance(o.theta(),float), 'Orbit.theta() does not return a scalar' + assert isinstance(o.vtheta(),float), 'Orbit.vtheta() does not return a scalar' + assert isinstance(o.vr(),float), 'Orbit.vr() does not return a scalar' assert isinstance(o.ra(),float), 'Orbit.ra() does not return a scalar' assert isinstance(o.dec(),float), 'Orbit.dec() does not return a scalar' assert isinstance(o.ll(),float), 'Orbit.ll() does not return a scalar' @@ -3338,6 +3370,9 @@ def test_scalar_all(): assert isinstance(o.y(5.),float), 'Orbit.y() does not return a scalar' assert isinstance(o.vx(5.),float), 'Orbit.vx() does not return a scalar' assert isinstance(o.vy(5.),float), 'Orbit.vy() does not return a scalar' + assert isinstance(o.theta(5.),float), 'Orbit.theta() does not return a scalar' + assert isinstance(o.vtheta(5.),float), 'Orbit.vtheta() does not return a scalar' + assert isinstance(o.vr(5.),float), 'Orbit.vr() does not return a scalar' assert isinstance(o.ra(5.),float), 'Orbit.ra() does not return a scalar' assert isinstance(o.dec(5.),float), 'Orbit.dec() does not return a scalar' assert isinstance(o.ll(5.),float), 'Orbit.ll() does not return a scalar' @@ -3723,6 +3758,9 @@ def test_orbit_method_integrate_t_asQuantity_warning(): check_integrate_t_asQuantity_warning(o,'y') check_integrate_t_asQuantity_warning(o,'vx') check_integrate_t_asQuantity_warning(o,'vy') + check_integrate_t_asQuantity_warning(o,'theta') + check_integrate_t_asQuantity_warning(o,'vtheta') + check_integrate_t_asQuantity_warning(o,'vr') check_integrate_t_asQuantity_warning(o,'ra') check_integrate_t_asQuantity_warning(o,'dec') check_integrate_t_asQuantity_warning(o,'ll') @@ -3790,6 +3828,9 @@ def test_orbit_method_inputro_quantity(): assert numpy.fabs(o.y(ro=ro*units.kpc)-o.y(ro=ro)) < 10.**-8., 'Orbit method y does not return the correct value when input ro is Quantity' assert numpy.fabs(o.vx(ro=ro*units.kpc)-o.vx(ro=ro)) < 10.**-8., 'Orbit method vx does not return the correct value when input ro is Quantity' assert numpy.fabs(o.vy(ro=ro*units.kpc)-o.vy(ro=ro)) < 10.**-8., 'Orbit method vy does not return the correct value when input ro is Quantity' + assert numpy.fabs(o.theta(ro=ro*units.kpc)-o.theta(ro=ro)) < 10.**-8., 'Orbit method theta does not return the correct value when input ro is Quantity' + assert numpy.fabs(o.vtheta(ro=ro*units.kpc)-o.vtheta(ro=ro)) < 10.**-8., 'Orbit method vtheta does not return the correct value when input ro is Quantity' + assert numpy.fabs(o.vr(ro=ro*units.kpc)-o.vr(ro=ro)) < 10.**-8., 'Orbit method vr does not return the correct value when input ro is Quantity' assert numpy.fabs(o.ra(ro=ro*units.kpc)-o.ra(ro=ro)) < 10.**-8., 'Orbit method ra does not return the correct value when input ro is Quantity' assert numpy.fabs(o.dec(ro=ro*units.kpc)-o.dec(ro=ro)) < 10.**-8., 'Orbit method dec does not return the correct value when input ro is Quantity' assert numpy.fabs(o.ll(ro=ro*units.kpc)-o.ll(ro=ro)) < 10.**-8., 'Orbit method ll does not return the correct value when input ro is Quantity' @@ -3852,6 +3893,9 @@ def test_orbit_method_inputvo_quantity(): assert numpy.fabs(o.y(vo=vo*units.km/units.s)-o.y(vo=vo)) < 10.**-8., 'Orbit method y does not return the correct value when input vo is Quantity' assert numpy.fabs(o.vx(vo=vo*units.km/units.s)-o.vx(vo=vo)) < 10.**-8., 'Orbit method vx does not return the correct value when input vo is Quantity' assert numpy.fabs(o.vy(vo=vo*units.km/units.s)-o.vy(vo=vo)) < 10.**-8., 'Orbit method vy does not return the correct value when input vo is Quantity' + assert numpy.fabs(o.theta(vo=vo*units.km/units.s)-o.theta(vo=vo)) < 10.**-8., 'Orbit method theta does not return the correct value when input vo is Quantity' + assert numpy.fabs(o.vtheta(vo=vo*units.km/units.s)-o.vtheta(vo=vo)) < 10.**-8., 'Orbit method vtheta does not return the correct value when input vo is Quantity' + assert numpy.fabs(o.vr(vo=vo*units.km/units.s)-o.vr(vo=vo)) < 10.**-8., 'Orbit method vr does not return the correct value when input vo is Quantity' assert numpy.fabs(o.ra(vo=vo*units.km/units.s)-o.ra(vo=vo)) < 10.**-8., 'Orbit method ra does not return the correct value when input vo is Quantity' assert numpy.fabs(o.dec(vo=vo*units.km/units.s)-o.dec(vo=vo)) < 10.**-8., 'Orbit method dec does not return the correct value when input vo is Quantity' assert numpy.fabs(o.ll(vo=vo*units.km/units.s)-o.ll(vo=vo)) < 10.**-8., 'Orbit method ll does not return the correct value when input vo is Quantity'