Implementing new potentials (bulge and gas disk models)

[JovanVeljanoski]

Hello!

I've been using galpy for some time now (very happy btw), and recently I got in great need of two new potentials, one that would model a bulge and one for a gas disk. The ideas behind these are already in place in the Wiki entry under "Possible galpy extensions"

• TwoPowerSphericalPotentials with flattening in the density (to which I would also add an exponential cut-off) (halo model in Binney & Tremaine 2008)
• Double exponential disk with a central hole (ISM model in Binney & Tremaine 2008)

I am happy to give it a go and try to implement one or both of these potentials. I attempted to do so for the Double exponential disk with a central hole (assuming this would be easy to do, i.e. modify the Double exponential disk that is already implemented), but I got stuck at how the potential is evaluated.

Can you please explain in some detail how is this done, or maybe point me towards some derivation of the functional form for the evaluation of this potential. I would greatly appreciate it if you could explain the meaning of the variables in the part of the code that evaluates the DoubleExponentialDiskPotential. I attach the relevant code below.

Thank you very much! J.

class DoubleExponentialDiskPotential(Potential): """Class that implements the double exponential disk potential

.. math::

    \\rho(R,z) = \\mathrm{amp}\\,\\exp\\left(-R/h_R-|z|/h_z\\right)

"""
def __init__(self,amp=1.,hr=1./3.,hz=1./16.,
             maxiter=_MAXITER,tol=0.001,normalize=False,
             ro=None,vo=None,
             new=True,kmaxFac=2.,glorder=10):
    """
    NAME:

       __init__

    PURPOSE:

       initialize a double-exponential disk potential

    INPUT:

       amp - amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density

       hr - disk scale-length (can be Quantity)

       hz - scale-height (can be Quantity)

       tol - relative accuracy of potential-evaluations

       maxiter - scipy.integrate keyword

       normalize - if True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.

       ro=, vo= distance and velocity scales for translation into internal units (default from configuration file)

    OUTPUT:

       DoubleExponentialDiskPotential object

    HISTORY:

       2010-04-16 - Written - Bovy (NYU)

       2013-01-01 - Re-implemented using faster integration techniques - Bovy (IAS)

    """
    Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units='density')
    if _APY_LOADED and isinstance(hr,units.Quantity):
        hr= hr.to(units.kpc).value/self._ro
    if _APY_LOADED and isinstance(hz,units.Quantity):
        hz= hz.to(units.kpc).value/self._ro
    self.hasC= True
    self._kmaxFac= kmaxFac
    self._glorder= glorder
    self._hr= hr
    self._scale= self._hr
    self._hz= hz
    self._alpha= 1./self._hr
    self._beta= 1./self._hz
    self._gamma= self._alpha/self._beta
    self._maxiter= maxiter
    self._tol= tol
    self._zforceNotSetUp= True #We have not calculated a typical Kz yet
    #Setup j0 zeros etc.
    self._glx, self._glw= nu.polynomial.legendre.leggauss(self._glorder)
    self._nzeros=100
    #j0 for potential and z
    self._j0zeros= nu.zeros(self._nzeros+1)
    self._j0zeros[1:self._nzeros+1]= special.jn_zeros(0,self._nzeros)
    self._dj0zeros= self._j0zeros-nu.roll(self._j0zeros,1)
    self._dj0zeros[0]= self._j0zeros[0]
    #j1 for R
    self._j1zeros= nu.zeros(self._nzeros+1)
    self._j1zeros[1:self._nzeros+1]= special.jn_zeros(1,self._nzeros)
    self._dj1zeros= self._j1zeros-nu.roll(self._j1zeros,1)
    self._dj1zeros[0]= self._j1zeros[0]
    #j2 for R2deriv
    self._j2zeros= nu.zeros(self._nzeros+1)
    self._j2zeros[1:self._nzeros+1]= special.jn_zeros(2,self._nzeros)
    self._dj2zeros= self._j2zeros-nu.roll(self._j2zeros,1)
    self._dj2zeros[0]= self._j2zeros[0]
    if normalize or \
            (isinstance(normalize,(int,float)) \
                 and not isinstance(normalize,bool)): #pragma: no cover
        self.normalize(normalize)
    #Load Kepler potential for large R
    self._kp= KeplerPotential(normalize=4.*nu.pi/self._alpha**2./self._beta)

def _evaluate(self,R,z,phi=0.,t=0.,dR=0,dphi=0):
    """
    NAME:
       _evaluate
    PURPOSE:
       evaluate the potential at (R,z)
    INPUT:
       R - Cylindrical Galactocentric radius
       z - vertical height
       phi - azimuth
       t - time
    OUTPUT:
       potential at (R,z)
    HISTORY:
       2010-04-16 - Written - Bovy (NYU)
       2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS)
    DOCTEST:
       >>> doubleExpPot= DoubleExponentialDiskPotential()
       >>> r= doubleExpPot(1.,0) #doctest: +ELLIPSIS
       ...
       >>> assert( r+1.89595350484)**2.< 10.**-6.
    """
    if True:
        if isinstance(R,float):
            floatIn= True
            R= nu.array([R])
            z= nu.array([z])
        else:
            floatIn= False
        out= nu.empty(len(R))
        indx= (R <= 6.)
        out[True-indx]= self._kp(R[True-indx],z[True-indx])
        R4max= nu.copy(R)
        R4max[(R < 1.)]= 1.
        kmax= self._kmaxFac*self._beta
        for jj in range(len(R)):
            if not indx[jj]: continue
            maxj0zeroIndx= nu.argmin((self._j0zeros-kmax*R4max[jj])**2.) #close enough
            ks= nu.array([0.5*(self._glx+1.)*self._dj0zeros[ii+1] + self._j0zeros[ii] for ii in range(maxj0zeroIndx)]).flatten()
            weights= nu.array([self._glw*self._dj0zeros[ii+1] for ii in range(maxj0zeroIndx)]).flatten()
            evalInt= special.jn(0,ks*R[jj])*(self._alpha**2.+ks**2.)**-1.5*(self._beta*nu.exp(-ks*nu.fabs(z[jj]))-ks*nu.exp(-self._beta*nu.fabs(z[jj])))/(self._beta**2.-ks**2.)
            out[jj]= -2.*nu.pi*self._alpha*nu.sum(weights*evalInt)
        if floatIn: return out[0]
        else: return out

[jobovy]

Dear Jovan,

It would be great if you would add one of these potentials.

I actually just started working on the generalization of the TwoPowerSphericalPotentials, allowing them to be triaxial (so have flattening in z and y) in a branch twopowertri. I'm hoping to make some progress on this in the next few weeks. I wasn't planning to include an exponential cut-off, and I don't know whether that complicates the math; I'm using essentially the method from Binney & Tremaine (2008) 2.5.3 which requires the psi(m) function; I don't know whether this remains analytical when you add the exponential cut-off.

The potential with the central hole for the ISM would be a great addition. The DoubleExponentialDiskPotential implementation follows Appendix A of Kuijken & Gilmore (1989). My code performs the integrals over k in (A10-A12) by integrating between the zeros of the Bessel functions, because the integrand is oscillatory.

Going through the math of that Appendix, I think it should be straightforward to add the hole part of the potential. The evaluation of the potential and the forces will remain largely the same, except that the radial Hankel transformation [Eq. (A6-A7)] is no longer analytical and would have to be done numerically. This integral does not depend on the (R,z) pair where the potential is evaluated, so it would be possible to setup an interpolation grid in the __init__ function that can then be used in the _evaluate and _R/zforce functions. The Hankel transformation is somewhat tricky to calculate numerically, because of the oscillatory behavior of the Bessel function, but that can probably be handled by integrating between the zeros in a similar way as how the integrals in (A10-12) are handled in the code.

Alternatively, it might be good to think a little and see whether there is any surface-density profile that has a central hole that does have an analytical Hankel transformation. There is nothing special about the exp(-Rm/R) implementation of the hole [that does not exactly describe the hole], so if there is some function that is similar to this, but has an analytical Hankel transformation, that would be very useful. I spent a little bit of time going through some tables of Hankel transformations to find such a function, without any luck, but perhaps it would be possible to find one (or play with Mathematica or similar to try to find a function for which the Hankel transformation is analytical).

Anyway, I'd be happy to consult if you want to start working on this. The best way to start would be to fork the code, make a branch called dblholepot and start on an implementation of DoubleExponentialDiskwHolePotential.py along the lines of DoubleExponentialDiskPotential.py.

[jobovy]

Hi Jovan,

It's been a while since you opened this issue and I wanted to check whether you are working on this and making progress? Just so I know whether yo consider this work as something in progress in case somebody else wants to work on this. If you are not planning on working on this, could you close the issue?

Thanks!

Jo

[jobovy]

Closing this, because it's not a real issue with the code and the discussion is stale.

[jobovy]

For reference: general disk potentials are now supported through the DiskSCFPotential class; this includes the gas disk model from Dehnen & Binney (1998) / Binney & Tremaine (2008).