Source code for galpy.potential.plotRotcurve

from __future__ import division, print_function

import os
import pickle
import numpy
from ..util import plot, conversion
from ..util.conversion import physical_conversion, potential_physical_input
[docs]def plotRotcurve(Pot,*args,**kwargs): """ NAME: plotRotcurve PURPOSE: plot the rotation curve for this potential (in the z=0 plane for non-spherical potentials) INPUT: Pot - Potential or list of Potential instances Rrange - Range in R to consider (needs to be in the units that you are plotting; can be Quantity) grid= grid in R phi= (None) azimuth to use for non-axisymmetric potentials savefilename= save to or restore from this savefile (pickle) +galpy.util.plot.plot args and kwargs OUTPUT: plot to output device HISTORY: 2010-07-10 - Written - Bovy (NYU) 2016-06-15 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT) """ # Using physical units or not? if isinstance(Pot,list): potro= Pot[0]._ro roSet= Pot[0]._roSet potvo= Pot[0]._vo voSet= Pot[0]._voSet else: potro= Pot._ro roSet= Pot._roSet potvo= Pot._vo voSet= Pot._voSet if (kwargs.get('use_physical',False) \ and kwargs.get('ro',roSet) and kwargs.get('vo',voSet)) or \ (not 'use_physical' in kwargs \ and kwargs.get('ro',roSet) and kwargs.get('vo',voSet)): use_physical= True potro= kwargs.get('ro',potro) potvo= kwargs.get('vo',potvo) xlabel= r'$R\,(\mathrm{kpc})$' ylabel= r"$v_c(R)\,(\mathrm{km\,s}^{-1})$" Rrange= kwargs.pop('Rrange',[0.01*potro,5.*potro]) else: use_physical= False xlabel= r"$R/R_0$" ylabel= r"$v_c(R)/v_c(R_0)$" Rrange= kwargs.pop('Rrange',[0.01,5.]) # Parse ro potro= conversion.parse_length_kpc(potro) potvo= conversion.parse_velocity_kms(potvo) Rrange[0]= conversion.parse_length_kpc(Rrange[0]) Rrange[1]= conversion.parse_length_kpc(Rrange[1]) if use_physical: Rrange[0]/= potro Rrange[1]/= potro grid= kwargs.pop('grid',1001) savefilename= kwargs.pop('savefilename',None) phi= kwargs.pop('phi',None) if not savefilename is None and os.path.exists(savefilename): print("Restoring savefile "+savefilename+" ...") savefile= open(savefilename,'rb') rotcurve= pickle.load(savefile) Rs= pickle.load(savefile) savefile.close() else: Rs= numpy.linspace(Rrange[0],Rrange[1],grid) rotcurve= calcRotcurve(Pot,Rs,phi=phi) if not savefilename == None: print("Writing savefile "+savefilename+" ...") savefile= open(savefilename,'wb') pickle.dump(rotcurve,savefile) pickle.dump(Rs,savefile) savefile.close() if use_physical: Rs*= potro rotcurve*= potvo Rrange[0]*= potro Rrange[1]*= potro if not 'xlabel' in kwargs: kwargs['xlabel']= xlabel if not 'ylabel' in kwargs: kwargs['ylabel']= ylabel if not 'xrange' in kwargs: kwargs['xrange']= Rrange if not 'yrange' in kwargs: kwargs['yrange']= [0.,1.2*numpy.amax(rotcurve)] kwargs.pop('ro',None) kwargs.pop('vo',None) kwargs.pop('use_physical',None) return plot.plot(Rs,rotcurve,*args,**kwargs)
def calcRotcurve(Pot,Rs,phi=None,t=0.): """ NAME: calcRotcurve PURPOSE: calculate the rotation curve for this potential (in the z=0 plane for non-spherical potentials) INPUT: Pot - Potential or list of Potential instances Rs - (array of) radius(i) phi= (None) azimuth to use for non-axisymmetric potentials ts - instantaneous time (optional) OUTPUT: array of vc HISTORY: 2011-04-13 - Written - Bovy (NYU) 2016-06-15 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT) """ try: grid= len(Rs) except TypeError: grid=1 Rs= numpy.array([Rs]) rotcurve= numpy.zeros(grid) for ii in range(grid): rotcurve[ii]= vcirc(Pot,Rs[ii],phi=phi,t=t,use_physical=False) return rotcurve
[docs]@potential_physical_input @physical_conversion('velocity',pop=True) def vcirc(Pot,R,phi=None,t=0.): """ NAME: vcirc PURPOSE: calculate the circular velocity at R in potential Pot INPUT: Pot - Potential instance or list of such instances R - Galactocentric radius (can be Quantity) phi= (None) azimuth to use for non-axisymmetric potentials t= time (optional; can be Quantity) OUTPUT: circular rotation velocity HISTORY: 2011-10-09 - Written - Bovy (IAS) 2016-06-15 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT) """ from ..potential import evaluateplanarRforces from ..potential import PotentialError try: return numpy.sqrt(-R*evaluateplanarRforces(Pot,R,phi=phi,t=t, use_physical=False)) except PotentialError: from ..potential import toPlanarPotential Pot= toPlanarPotential(Pot) return numpy.sqrt(-R*evaluateplanarRforces(Pot,R,phi=phi,t=t, use_physical=False))
[docs]@potential_physical_input @physical_conversion('frequency',pop=True) def dvcircdR(Pot,R,phi=None,t=0.): """ NAME: dvcircdR PURPOSE: calculate the derivative of the circular velocity wrt R at R in potential Pot INPUT: Pot - Potential instance or list of such instances R - Galactocentric radius (can be Quantity) phi= (None) azimuth to use for non-axisymmetric potentials t= time (optional; can be Quantity) OUTPUT: derivative of the circular rotation velocity wrt R HISTORY: 2013-01-08 - Written - Bovy (IAS) 2016-06-28 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT) """ from ..potential import evaluateplanarRforces, evaluateplanarR2derivs from ..potential import PotentialError tvc= vcirc(Pot,R,phi=phi,t=t,use_physical=False) try: return 0.5*(-evaluateplanarRforces(Pot,R,phi=phi,t=t,use_physical=False)+R*evaluateplanarR2derivs(Pot,R,phi=phi,t=t,use_physical=False))/tvc except PotentialError: from ..potential import RZToplanarPotential Pot= RZToplanarPotential(Pot) return 0.5*(-evaluateplanarRforces(Pot,R,phi=phi,t=t,use_physical=False)+R*evaluateplanarR2derivs(Pot,R,phi=phi,t=t,use_physical=False))/tvc