def __init__(self,pot=None,delta=None,Rmax=5.,
                 nE=25,npsi=25,nLz=30,numcores=1,
                 interpecc=False,
                 **kwargs):
        """
        NAME:
           __init__
        PURPOSE:
           initialize an actionAngleStaeckelGrid object
        INPUT:
           pot= potential or list of potentials

           delta= focus of prolate confocal coordinate system (can be Quantity)

           Rmax = Rmax for building grids (natural units)

           nE=, npsi=, nLz= grid size

           interpecc= (False) if True, also interpolate the approximate eccentricity, zmax, rperi, and rapo

           numcores= number of cpus to use to parallellize

           ro= distance from vantage point to GC (kpc; can be Quantity)

           vo= circular velocity at ro (km/s; can be Quantity)

        OUTPUT:
         
           instance

        HISTORY:

            2012-11-29 - Written - Bovy (IAS)

            2017-12-15 - Written - Bovy (UofT)

        """
        actionAngle.__init__(self,
                             ro=kwargs.get('ro',None),vo=kwargs.get('vo',None))
        if pot is None:
            raise IOError("Must specify pot= for actionAngleStaeckelGrid")
        self._pot= flatten_potential(pot)
        if delta is None:
            raise IOError("Must specify delta= for actionAngleStaeckelGrid")
        if ext_loaded and 'c' in kwargs and kwargs['c']:
            self._c= True
        else:
            self._c= False
        self._delta= delta
        if _APY_LOADED and isinstance(self._delta,units.Quantity):
            self._delta= self._delta.to(units.kpc).value/self._ro
        self._Rmax= Rmax
        self._Rmin= 0.01
        #Set up the actionAngleStaeckel object that we will use to interpolate
        self._aA= actionAngleStaeckel.actionAngleStaeckel(pot=self._pot,delta=self._delta,c=self._c)
        #Build grid
        self._Lzmin= 0.01
        self._Lzs= numpy.linspace(self._Lzmin,
                                  self._Rmax\
                                      *galpy.potential.vcirc(self._pot,
                                                             self._Rmax),
                                  nLz)
        self._Lzmax= self._Lzs[-1]
        self._nLz= nLz
        #Calculate E_c(R=RL), energy of circular orbit
        self._RL= numpy.array([galpy.potential.rl(self._pot,l) for l in self._Lzs])
        self._RLInterp= interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                                 self._RL,k=3)
        self._ERL= _evaluatePotentials(self._pot,self._RL,
                                       numpy.zeros(self._nLz))\
                                       +self._Lzs**2./2./self._RL**2.
        self._ERLmax= numpy.amax(self._ERL)+1.
        self._ERLInterp= interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                                  numpy.log(-(self._ERL-self._ERLmax)),k=3)
        self._Ramax= 200./8.
        self._ERa= _evaluatePotentials(self._pot,self._Ramax,0.) +self._Lzs**2./2./self._Ramax**2.
        #self._EEsc= numpy.array([self._ERL[ii]+galpy.potential.vesc(self._pot,self._RL[ii])**2./4. for ii in range(nLz)])
        self._ERamax= numpy.amax(self._ERa)+1.
        self._ERaInterp= interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                                  numpy.log(-(self._ERa-self._ERamax)),k=3)
        y= numpy.linspace(0.,1.,nE)
        self._nE= nE
        psis= numpy.linspace(0.,1.,npsi)*numpy.pi/2.
        self._npsi= npsi
        jr= numpy.zeros((nLz,nE,npsi))
        jz= numpy.zeros((nLz,nE,npsi))
        u0= numpy.zeros((nLz,nE))
        jrLzE= numpy.zeros((nLz))
        jzLzE= numpy.zeros((nLz))
        #First calculate u0
        thisLzs= (numpy.tile(self._Lzs,(nE,1)).T).flatten()
        thisERL= (numpy.tile(self._ERL,(nE,1)).T).flatten()
        thisERa= (numpy.tile(self._ERa,(nE,1)).T).flatten()
        thisy= (numpy.tile(y,(nLz,1))).flatten()
        thisE= _invEfunc(_Efunc(thisERa,thisERL)+thisy*(_Efunc(thisERL,thisERL)-_Efunc(thisERa,thisERL)),thisERL)
        if isinstance(self._pot,galpy.potential.interpRZPotential) and hasattr(self._pot,'_origPot'):
            u0pot= self._pot._origPot
        else:
            u0pot= self._pot
        if self._c:
            mu0= actionAngleStaeckel_c.actionAngleStaeckel_calcu0(thisE,thisLzs,
                                                                  u0pot,
                                                                  self._delta)[0]
        else:
            if numcores > 1:
                mu0= multi.parallel_map((lambda x: self.calcu0(thisE[x],
                                                               thisLzs[x])),
                                        range(nE*nLz),
                                        numcores=numcores)
            else:
                mu0= list(map((lambda x: self.calcu0(thisE[x],
                                                     thisLzs[x])),
                              range(nE*nLz)))
        u0= numpy.reshape(mu0,(nLz,nE))
        thisR= self._delta*numpy.sinh(u0)
        thisv= numpy.reshape(self.vatu0(thisE.flatten(),thisLzs.flatten(),
                                        u0.flatten(),
                                        thisR.flatten()),(nLz,nE))
        self.thisv= thisv
        #reshape
        thisLzs= numpy.reshape(thisLzs,(nLz,nE))
        thispsi= numpy.tile(psis,(nLz,nE,1)).flatten()
        thisLzs= numpy.tile(thisLzs.T,(npsi,1,1)).T.flatten()
        thisR= numpy.tile(thisR.T,(npsi,1,1)).T.flatten()
        thisv= numpy.tile(thisv.T,(npsi,1,1)).T.flatten()
        mjr, mlz, mjz= self._aA(thisR, #R
                                thisv*numpy.cos(thispsi), #vR
                                thisLzs/thisR, #vT
                                numpy.zeros(len(thisR)), #z
                                thisv*numpy.sin(thispsi), #vz
                                fixed_quad=True) 
        if interpecc:
            mecc, mzmax, mrperi, mrap=\
                self._aA.EccZmaxRperiRap(thisR, #R
                                         thisv*numpy.cos(thispsi), #vR
                                         thisLzs/thisR, #vT
                                         numpy.zeros(len(thisR)), #z
                                         thisv*numpy.sin(thispsi)) #vz
        if isinstance(self._pot,galpy.potential.interpRZPotential) and hasattr(self._pot,'_origPot'):
            #Interpolated potentials have problems with extreme orbits
            indx= (mjr == 9999.99)
            indx+= (mjz == 9999.99)
            #Re-calculate these using the original potential, hopefully not too slow
            tmpaA= actionAngleStaeckel.actionAngleStaeckel(pot=self._pot._origPot,delta=self._delta,c=self._c)
            mjr[indx], dum, mjz[indx]= tmpaA(thisR[indx], #R
                                             thisv[indx]*numpy.cos(thispsi[indx]), #vR
                                             thisLzs[indx]/thisR[indx], #vT
                                             numpy.zeros(numpy.sum(indx)), #z
                                             thisv[indx]*numpy.sin(thispsi[indx]), #vz
                                             fixed_quad=True)
            if interpecc:
                mecc[indx], mzmax[indx], mrperi[indx], mrap[indx]=\
                    self._aA.EccZmaxRperiRap(thisR[indx], #R
                                             thisv[indx]*numpy.cos(thispsi[indx]), #vR
                                             thisLzs[indx]/thisR[indx], #vT
                                             numpy.zeros(numpy.sum(indx)), #z
                                             thisv[indx]*numpy.sin(thispsi[indx])) #vz
        jr= numpy.reshape(mjr,(nLz,nE,npsi))
        jz= numpy.reshape(mjz,(nLz,nE,npsi))
        if interpecc:
            ecc= numpy.reshape(mecc,(nLz,nE,npsi))
            zmax= numpy.reshape(mzmax,(nLz,nE,npsi))
            rperi= numpy.reshape(mrperi,(nLz,nE,npsi))
            rap= numpy.reshape(mrap,(nLz,nE,npsi))
            zmaxLzE= numpy.zeros((nLz))
            rperiLzE= numpy.zeros((nLz))
            rapLzE= numpy.zeros((nLz))
        for ii in range(nLz):
            jrLzE[ii]= numpy.nanmax(jr[ii,(jr[ii,:,:] != 9999.99)])#:,:])
            jzLzE[ii]= numpy.nanmax(jz[ii,(jz[ii,:,:] != 9999.99)])#:,:])
            if interpecc:
                zmaxLzE[ii]= numpy.amax(zmax[ii,numpy.isfinite(zmax[ii])])
                rperiLzE[ii]= numpy.amax(rperi[ii,numpy.isfinite(rperi[ii])])
                rapLzE[ii]= numpy.amax(rap[ii,numpy.isfinite(rap[ii])])
        jrLzE[(jrLzE == 0.)]= numpy.nanmin(jrLzE[(jrLzE > 0.)])
        jzLzE[(jzLzE == 0.)]= numpy.nanmin(jzLzE[(jzLzE > 0.)])
        if interpecc:
            zmaxLzE[(zmaxLzE == 0.)]= numpy.nanmin(zmaxLzE[(zmaxLzE > 0.)])
            rperiLzE[(rperiLzE == 0.)]= numpy.nanmin(rperiLzE[(rperiLzE > 0.)])
            rapLzE[(rapLzE == 0.)]= numpy.nanmin(rapLzE[(rapLzE > 0.)])
        for ii in range(nLz):
            jr[ii,:,:]/= jrLzE[ii]
            jz[ii,:,:]/= jzLzE[ii]
            if interpecc:
                zmax[ii,:,:]/= zmaxLzE[ii]
                rperi[ii,:,:]/= rperiLzE[ii]
                rap[ii,:,:]/= rapLzE[ii]
        #Deal w/ 9999.99
        jr[(jr > 1.)]= 1.
        jz[(jz > 1.)]= 1.
        #Deal w/ NaN
        jr[numpy.isnan(jr)]= 0.
        jz[numpy.isnan(jz)]= 0.
        if interpecc:
            ecc[(ecc < 0.)]= 0.
            ecc[(ecc > 1.)]= 1.
            ecc[numpy.isnan(ecc)]= 0.
            ecc[numpy.isinf(ecc)]= 1.
            zmax[(zmax > 1.)]= 1.
            zmax[numpy.isnan(zmax)]= 0.
            zmax[numpy.isinf(zmax)]= 1.
            rperi[(rperi > 1.)]= 1.
            rperi[numpy.isnan(rperi)]= 0.
            rperi[numpy.isinf(rperi)]= 0. # typically orbits that can reach 0
            rap[(rap > 1.)]= 1.
            rap[numpy.isnan(rap)]= 0.
            rap[numpy.isinf(rap)]= 1.
        #First interpolate the maxima
        self._jr= jr
        self._jz= jz
        self._u0= u0
        self._jrLzE= jrLzE
        self._jzLzE= jzLzE
        self._jrLzInterp= interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                                   numpy.log(jrLzE+10.**-5.),k=3)
        self._jzLzInterp= interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                                   numpy.log(jzLzE+10.**-5.),k=3)
        if interpecc:
            self._ecc= ecc
            self._zmax= zmax
            self._rperi= rperi
            self._rap= rap
            self._zmaxLzE= zmaxLzE
            self._rperiLzE= rperiLzE
            self._rapLzE= rapLzE
            self._zmaxLzInterp=\
                interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                         numpy.log(zmaxLzE+10.**-5.),k=3)
            self._rperiLzInterp=\
                interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                         numpy.log(rperiLzE+10.**-5.),k=3)
            self._rapLzInterp=\
                interpolate.InterpolatedUnivariateSpline(self._Lzs,
                                                         numpy.log(rapLzE+10.**-5.),k=3)
        #Interpolate u0
        self._logu0Interp= interpolate.RectBivariateSpline(self._Lzs,
                                                           y,
                                                           numpy.log(u0),
                                                           kx=3,ky=3,s=0.)
        #spline filter jr and jz, such that they can be used with ndimage.map_coordinates
        self._jrFiltered= ndimage.spline_filter(numpy.log(self._jr+10.**-10.),order=3)
        self._jzFiltered= ndimage.spline_filter(numpy.log(self._jz+10.**-10.),order=3)
        if interpecc:
            self._eccFiltered= ndimage.spline_filter(numpy.log(self._ecc+10.**-10.),order=3)
            self._zmaxFiltered= ndimage.spline_filter(numpy.log(self._zmax+10.**-10.),order=3)
            self._rperiFiltered= ndimage.spline_filter(numpy.log(self._rperi+10.**-10.),order=3)
            self._rapFiltered= ndimage.spline_filter(numpy.log(self._rap+10.**-10.),order=3)
        # Check the units
        self._check_consistent_units()
        return None
def __init__ ( self , pot = None , delta = None , Rmax = 5. , nE = 25 , npsi = 25 , nLz = 30 , numcores = 1 , interpecc = False , ** kwargs ) : actionAngle . __init__ ( self , ro = kwargs . get ( 'ro' , None ) , vo = kwargs . get ( 'vo' , None ) ) if pot is None : raise IOError ( "Must specify pot= for actionAngleStaeckelGrid" ) self . _pot = flatten_potential ( pot ) if delta is None : raise IOError ( "Must specify delta= for actionAngleStaeckelGrid" ) if ext_loaded and 'c' in kwargs and kwargs [ 'c' ] : self . _c = True else : self . _c = False self . _delta = delta if _APY_LOADED and isinstance ( self . _delta , units . Quantity ) : self . _delta = self . _delta . to ( units . kpc ) . value / self . _ro self . _Rmax = Rmax self . _Rmin = 0.01 self . _aA = actionAngleStaeckel . actionAngleStaeckel ( pot = self . _pot , delta = self . _delta , c = self . _c ) self . _Lzmin = 0.01 self . _Lzs = numpy . linspace ( self . _Lzmin , self . _Rmax * galpy . potential . vcirc ( self . _pot , self . _Rmax ) , nLz ) self . _Lzmax = self . _Lzs [ - 1 ] self . _nLz = nLz self . _RL = numpy . array ( [ galpy . potential . rl ( self . _pot , l ) for l in self . _Lzs ] ) self . _RLInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , self . _RL , k = 3 ) self . _ERL = _evaluatePotentials ( self . _pot , self . _RL , numpy . zeros ( self . _nLz ) ) + self . _Lzs ** 2. / 2. / self . _RL ** 2. self . _ERLmax = numpy . amax ( self . _ERL ) + 1. self . _ERLInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( - ( self . _ERL - self . _ERLmax ) ) , k = 3 ) self . _Ramax = 200. / 8. self . _ERa = _evaluatePotentials ( self . _pot , self . _Ramax , 0. ) + self . _Lzs ** 2. / 2. / self . _Ramax ** 2. self . _ERamax = numpy . amax ( self . _ERa ) + 1. self . _ERaInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( - ( self . _ERa - self . _ERamax ) ) , k = 3 ) y = numpy . linspace ( 0. , 1. , nE ) self . _nE = nE psis = numpy . linspace ( 0. , 1. , npsi ) * numpy . pi / 2. self . _npsi = npsi jr = numpy . zeros ( ( nLz , nE , npsi ) ) jz = numpy . zeros ( ( nLz , nE , npsi ) ) u0 = numpy . zeros ( ( nLz , nE ) ) jrLzE = numpy . zeros ( ( nLz ) ) jzLzE = numpy . zeros ( ( nLz ) ) thisLzs = ( numpy . tile ( self . _Lzs , ( nE , 1 ) ) . T ) . flatten ( ) thisERL = ( numpy . tile ( self . _ERL , ( nE , 1 ) ) . T ) . flatten ( ) thisERa = ( numpy . tile ( self . _ERa , ( nE , 1 ) ) . T ) . flatten ( ) thisy = ( numpy . tile ( y , ( nLz , 1 ) ) ) . flatten ( ) thisE = _invEfunc ( _Efunc ( thisERa , thisERL ) + thisy * ( _Efunc ( thisERL , thisERL ) - _Efunc ( thisERa , thisERL ) ) , thisERL ) if isinstance ( self . _pot , galpy . potential . interpRZPotential ) and hasattr ( self . _pot , '_origPot' ) : u0pot = self . _pot . _origPot else : u0pot = self . _pot if self . _c : mu0 = actionAngleStaeckel_c . actionAngleStaeckel_calcu0 ( thisE , thisLzs , u0pot , self . _delta ) [ 0 ] else : if numcores > 1 : mu0 = multi . parallel_map ( ( lambda x : self . calcu0 ( thisE [ x ] , thisLzs [ x ] ) ) , range ( nE * nLz ) , numcores = numcores ) else : mu0 = list ( map ( ( lambda x : self . calcu0 ( thisE [ x ] , thisLzs [ x ] ) ) , range ( nE * nLz ) ) ) u0 = numpy . reshape ( mu0 , ( nLz , nE ) ) thisR = self . _delta * numpy . sinh ( u0 ) thisv = numpy . reshape ( self . vatu0 ( thisE . flatten ( ) , thisLzs . flatten ( ) , u0 . flatten ( ) , thisR . flatten ( ) ) , ( nLz , nE ) ) self . thisv = thisv thisLzs = numpy . reshape ( thisLzs , ( nLz , nE ) ) thispsi = numpy . tile ( psis , ( nLz , nE , 1 ) ) . flatten ( ) thisLzs = numpy . tile ( thisLzs . T , ( npsi , 1 , 1 ) ) . T . flatten ( ) thisR = numpy . tile ( thisR . T , ( npsi , 1 , 1 ) ) . T . flatten ( ) thisv = numpy . tile ( thisv . T , ( npsi , 1 , 1 ) ) . T . flatten ( ) mjr , mlz , mjz = self . _aA ( thisR , thisv * numpy . cos ( thispsi ) , thisLzs / thisR , numpy . zeros ( len ( thisR ) ) , thisv * numpy . sin ( thispsi ) , fixed_quad = True ) if interpecc : mecc , mzmax , mrperi , mrap = self . _aA . EccZmaxRperiRap ( thisR , thisv * numpy . cos ( thispsi ) , thisLzs / thisR , numpy . zeros ( len ( thisR ) ) , thisv * numpy . sin ( thispsi ) ) if isinstance ( self . _pot , galpy . potential . interpRZPotential ) and hasattr ( self . _pot , '_origPot' ) : indx = ( mjr == 9999.99 ) indx += ( mjz == 9999.99 ) tmpaA = actionAngleStaeckel . actionAngleStaeckel ( pot = self . _pot . _origPot , delta = self . _delta , c = self . _c ) mjr [ indx ] , dum , mjz [ indx ] = tmpaA ( thisR [ indx ] , thisv [ indx ] * numpy . cos ( thispsi [ indx ] ) , thisLzs [ indx ] / thisR [ indx ] , numpy . zeros ( numpy . sum ( indx ) ) , thisv [ indx ] * numpy . sin ( thispsi [ indx ] ) , fixed_quad = True ) if interpecc : mecc [ indx ] , mzmax [ indx ] , mrperi [ indx ] , mrap [ indx ] = self . _aA . EccZmaxRperiRap ( thisR [ indx ] , thisv [ indx ] * numpy . cos ( thispsi [ indx ] ) , thisLzs [ indx ] / thisR [ indx ] , numpy . zeros ( numpy . sum ( indx ) ) , thisv [ indx ] * numpy . sin ( thispsi [ indx ] ) ) jr = numpy . reshape ( mjr , ( nLz , nE , npsi ) ) jz = numpy . reshape ( mjz , ( nLz , nE , npsi ) ) if interpecc : ecc = numpy . reshape ( mecc , ( nLz , nE , npsi ) ) zmax = numpy . reshape ( mzmax , ( nLz , nE , npsi ) ) rperi = numpy . reshape ( mrperi , ( nLz , nE , npsi ) ) rap = numpy . reshape ( mrap , ( nLz , nE , npsi ) ) zmaxLzE = numpy . zeros ( ( nLz ) ) rperiLzE = numpy . zeros ( ( nLz ) ) rapLzE = numpy . zeros ( ( nLz ) ) for ii in range ( nLz ) : jrLzE [ ii ] = numpy . nanmax ( jr [ ii , ( jr [ ii , : , : ] != 9999.99 ) ] ) jzLzE [ ii ] = numpy . nanmax ( jz [ ii , ( jz [ ii , : , : ] != 9999.99 ) ] ) if interpecc : zmaxLzE [ ii ] = numpy . amax ( zmax [ ii , numpy . isfinite ( zmax [ ii ] ) ] ) rperiLzE [ ii ] = numpy . amax ( rperi [ ii , numpy . isfinite ( rperi [ ii ] ) ] ) rapLzE [ ii ] = numpy . amax ( rap [ ii , numpy . isfinite ( rap [ ii ] ) ] ) jrLzE [ ( jrLzE == 0. ) ] = numpy . nanmin ( jrLzE [ ( jrLzE > 0. ) ] ) jzLzE [ ( jzLzE == 0. ) ] = numpy . nanmin ( jzLzE [ ( jzLzE > 0. ) ] ) if interpecc : zmaxLzE [ ( zmaxLzE == 0. ) ] = numpy . nanmin ( zmaxLzE [ ( zmaxLzE > 0. ) ] ) rperiLzE [ ( rperiLzE == 0. ) ] = numpy . nanmin ( rperiLzE [ ( rperiLzE > 0. ) ] ) rapLzE [ ( rapLzE == 0. ) ] = numpy . nanmin ( rapLzE [ ( rapLzE > 0. ) ] ) for ii in range ( nLz ) : jr [ ii , : , : ] /= jrLzE [ ii ] jz [ ii , : , : ] /= jzLzE [ ii ] if interpecc : zmax [ ii , : , : ] /= zmaxLzE [ ii ] rperi [ ii , : , : ] /= rperiLzE [ ii ] rap [ ii , : , : ] /= rapLzE [ ii ] jr [ ( jr > 1. ) ] = 1. jz [ ( jz > 1. ) ] = 1. jr [ numpy . isnan ( jr ) ] = 0. jz [ numpy . isnan ( jz ) ] = 0. if interpecc : ecc [ ( ecc < 0. ) ] = 0. ecc [ ( ecc > 1. ) ] = 1. ecc [ numpy . isnan ( ecc ) ] = 0. ecc [ numpy . isinf ( ecc ) ] = 1. zmax [ ( zmax > 1. ) ] = 1. zmax [ numpy . isnan ( zmax ) ] = 0. zmax [ numpy . isinf ( zmax ) ] = 1. rperi [ ( rperi > 1. ) ] = 1. rperi [ numpy . isnan ( rperi ) ] = 0. rperi [ numpy . isinf ( rperi ) ] = 0. rap [ ( rap > 1. ) ] = 1. rap [ numpy . isnan ( rap ) ] = 0. rap [ numpy . isinf ( rap ) ] = 1. self . _jr = jr self . _jz = jz self . _u0 = u0 self . _jrLzE = jrLzE self . _jzLzE = jzLzE self . _jrLzInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( jrLzE + 10. ** - 5. ) , k = 3 ) self . _jzLzInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( jzLzE + 10. ** - 5. ) , k = 3 ) if interpecc : self . _ecc = ecc self . _zmax = zmax self . _rperi = rperi self . _rap = rap self . _zmaxLzE = zmaxLzE self . _rperiLzE = rperiLzE self . _rapLzE = rapLzE self . _zmaxLzInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( zmaxLzE + 10. ** - 5. ) , k = 3 ) self . _rperiLzInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( rperiLzE + 10. ** - 5. ) , k = 3 ) self . _rapLzInterp = interpolate . InterpolatedUnivariateSpline ( self . _Lzs , numpy . log ( rapLzE + 10. ** - 5. ) , k = 3 ) self . _logu0Interp = interpolate . RectBivariateSpline ( self . _Lzs , y , numpy . log ( u0 ) , kx = 3 , ky = 3 , s = 0. ) self . _jrFiltered = ndimage . spline_filter ( numpy . log ( self . _jr + 10. ** - 10. ) , order = 3 ) self . _jzFiltered = ndimage . spline_filter ( numpy . log ( self . _jz + 10. ** - 10. ) , order = 3 ) if interpecc : self . _eccFiltered = ndimage . spline_filter ( numpy . log ( self . _ecc + 10. ** - 10. ) , order = 3 ) self . _zmaxFiltered = ndimage . spline_filter ( numpy . log ( self . _zmax + 10. ** - 10. ) , order = 3 ) self . _rperiFiltered = ndimage . spline_filter ( numpy . log ( self . _rperi + 10. ** - 10. ) , order = 3 ) self . _rapFiltered = ndimage . spline_filter ( numpy . log ( self . _rap + 10. ** - 10. ) , order = 3 ) self . _check_consistent_units ( ) return None
NAME:
           __init__
        PURPOSE:
           initialize an actionAngleStaeckelGrid object
        INPUT:
           pot= potential or list of potentials

           delta= focus of prolate confocal coordinate system (can be Quantity)

           Rmax = Rmax for building grids (natural units)

           nE=, npsi=, nLz= grid size

           interpecc= (False) if True, also interpolate the approximate eccentricity, zmax, rperi, and rapo

           numcores= number of cpus to use to parallellize

           ro= distance from vantage point to GC (kpc; can be Quantity)

           vo= circular velocity at ro (km/s; can be Quantity)

        OUTPUT:
         
           instance

        HISTORY:

            2012-11-29 - Written - Bovy (IAS)

            2017-12-15 - Written - Bovy (UofT)