Merge c1bb5bd into 549508a

doc/source/actionAngle.rst

@@ -1,3 +1,5 @@
+.. _actionangle:
+
 Action-angle coordinates
 =========================
 

doc/source/examples/dierickx_eccentricities.py

@@ -0,0 +1,192 @@
+import os, os.path
+import pexpect
+import subprocess
+from astropy.io import fits, ascii
+from astropy import units
+import numpy as np
+from optparse import OptionParser
+import sys
+import tempfile
+from ftplib import FTP
+import shutil
+import pickle
+from tqdm import tqdm
+from galpy.potential import LogarithmicHaloPotential
+from galpy.potential import evaluatePotentials as evalPot
+from galpy.orbit import Orbit
+from galpy.actionAngle import estimateDeltaStaeckel, actionAngleStaeckel, UnboundError
+from galpy.util import bovy_coords
+import matplotlib.pyplot as plt
+import numpy
+
+_ERASESTR= "                                                                                "
+
+def calc_eccentricity(args, options):
+    table = os.path.join(args[0],'table2.dat')
+    readme = os.path.join(args[0],'ReadMe')
+    dierickx = ascii.read(table, readme=readme)
+    vxvv = np.dstack([dierickx['RAdeg'], dierickx['DEdeg'], dierickx['Dist']/1e3, dierickx['pmRA'], dierickx['pmDE'], dierickx['HRV']])[0]
+    ro, vo, zo = 8., 220., 0.025
+    ra, dec= vxvv[:,0], vxvv[:,1]
+    lb= bovy_coords.radec_to_lb(ra,dec,degree=True)
+    pmra, pmdec= vxvv[:,3], vxvv[:,4]
+    pmllpmbb= bovy_coords.pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec,degree=True)
+    d, vlos= vxvv[:,2], vxvv[:,5]
+    rectgal= bovy_coords.sphergal_to_rectgal(lb[:,0],lb[:,1],d,vlos,pmllpmbb[:,0], pmllpmbb[:,1],degree=True)
+    vsolar= np.array([-10.1,4.0,6.7])
+    vsun= np.array([0.,1.,0.,])+vsolar/vo
+    X = rectgal[:,0]/ro
+    Y = rectgal[:,1]/ro
+    Z = rectgal[:,2]/ro
+    vx = rectgal[:,3]/vo
+    vy = rectgal[:,4]/vo
+    vz = rectgal[:,5]/vo
+    vsun= np.array([0.,1.,0.,])+vsolar/vo
+    Rphiz= bovy_coords.XYZ_to_galcencyl(X,Y,Z,Zsun=zo/ro)
+    vRvTvz= bovy_coords.vxvyvz_to_galcencyl(vx,vy,vz,Rphiz[:,0],Rphiz[:,1],Rphiz[:,2],vsun=vsun,Xsun=1.,Zsun=zo/ro,galcen=True)
+    #do the integration and individual analytic estimate for each object
+    ts= np.linspace(0.,20.,10000)
+    lp= LogarithmicHaloPotential(normalize=1.)
+    e_ana = numpy.zeros(len(vxvv))
+    e_int = numpy.zeros(len(vxvv))
+    print('Performing orbit integration and analytic parameter estimates for Dierickx et al. sample...')
+    for i in tqdm(range(len(vxvv))):
+        try:
+            orbit = Orbit(vxvv[i], radec=True, vo=220., ro=8.)
+            e_ana[i] = orbit.e(analytic=True, pot=lp, c=True)
+        except UnboundError:
+            e_ana[i] = np.nan
+        orbit.integrate(ts, lp)
+        e_int[i] = orbit.e(analytic=False)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(e_int, e_ana,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{galpy\ integrated}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ analytic}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-integratedeanalytice.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.hist(e_int, bins=30)
+    plt.xlim(0.,1.)
+    plt.xlabel(r'$\mathrm{galpy}\ e$')
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0], 'dierickx-integratedehist.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(dierickx['e'], e_int,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ integrated}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-integratedee.png'), format='png', dpi=200)
+    fig = plt.figure()
+    fig.set_size_inches(1.5*columnwidth, 1.5*columnwidth)
+    plt.scatter(dierickx['e'], e_ana,  s=1, color='Black', lw=0.)
+    plt.xlabel(r'$\mathrm{Dierickx\ et\ al.}\ e$')
+    plt.ylabel(r'$\mathrm{galpy\ estimated}\ e$')
+    plt.xlim(0.,1.)
+    plt.ylim(0.,1.)
+    fig.tight_layout()
+    plt.savefig(os.path.join(args[0],'dierickx-analyticee.png'), format='png', dpi=200)
+    arr = numpy.recarray(len(e_ana), dtype=[('analytic_e', float), ('integrated_e', float)])
+    arr['analytic_e'] = e_ana
+    arr['integrated_e'] = e_int
+    with open(os.path.join(args[0],'eccentricities.dat'), 'w') as file:
+        pickle.dump(arr, file)
+        file.close()
+
+def get_table(args,options):
+    cat = 'J/ApJ/725/L186/'
+    tab2name = 'table2.dat.gz'
+    tab2readme = 'ReadMe'
+    out = args[0]
+    ensure_dir(os.path.join(out,tab2name))
+    vizier(cat, os.path.join(out,tab2name), os.path.join(out,tab2readme), catalogname=tab2name, readmename=tab2readme)
+    subprocess.call(['gunzip', os.path.join(out,tab2name)])
+
+def vizier(cat,filePath,ReadMePath,
+           catalogname='catalog.dat',readmename='ReadMe'):
+    """
+    NAME:
+       vizier
+    PURPOSE:
+       download a catalog and its associated ReadMe from Vizier
+    INPUT:
+       cat - name of the catalog (e.g., 'III/272' for RAVE, or J/A+A/... for journal-specific catalogs)
+       filePath - path of the file where you want to store the catalog (note: you need to keep the name of the file the same as the catalogname to be able to read the file with astropy.io.ascii)
+       ReadMePath - path of the file where you want to store the ReadMe file
+       catalogname= (catalog.dat) name of the catalog on the Vizier server
+       readmename= (ReadMe) name of the ReadMe file on the Vizier server
+    OUTPUT:
+       (nothing, just downloads)
+    HISTORY:
+       2016-09-12 - Written - Bovy (UofT)
+    """
+    _download_file_vizier(cat,filePath,catalogname=catalogname)
+    _download_file_vizier(cat,ReadMePath,catalogname=readmename)
+    return None
+
+
+def _download_file_vizier(cat,filePath,catalogname='catalog.dat'):
+    '''
+    Stolen from Jo Bovy's gaia_tools package!
+    '''
+    sys.stdout.write('\r'+"Downloading file %s ...\r" \
+                         % (os.path.basename(filePath)))
+    sys.stdout.flush()
+    try:
+        # make all intermediate directories
+        os.makedirs(os.path.dirname(filePath)) 
+    except OSError: pass
+    # Safe way of downloading
+    downloading= True
+    interrupted= False
+    file, tmp_savefilename= tempfile.mkstemp()
+    os.close(file) #Easier this way
+    ntries= 1
+    while downloading:
+        try:
+            ftp= FTP('cdsarc.u-strasbg.fr')
+            ftp.login('anonymous', 'test')
+            ftp.cwd(os.path.join('pub','cats',cat))
+            with open(tmp_savefilename,'wb') as savefile:
+                ftp.retrbinary('RETR %s' % catalogname,savefile.write)
+            shutil.move(tmp_savefilename,filePath)
+            downloading= False
+            if interrupted:
+                raise KeyboardInterrupt
+        except:
+            raise
+            if not downloading: #Assume KeyboardInterrupt
+                raise
+            elif ntries > _MAX_NTRIES:
+                raise IOError('File %s does not appear to exist on the server ...' % (os.path.basename(filePath)))
+        finally:
+            if os.path.exists(tmp_savefilename):
+                os.remove(tmp_savefilename)
+        ntries+= 1
+    sys.stdout.write('\r'+_ERASESTR+'\r')
+    sys.stdout.flush()        
+    return None
+
+def ensure_dir(f):
+    """ Ensure a a file exists and if not make the relevant path """
+    d = os.path.dirname(f)
+    if not os.path.exists(d):
+        os.makedirs(d)
+
+def get_options():
+    #no options yet - probably none needed?
+    usage = "usage: %prog [options] <outpath>"
+    parser = OptionParser(usage=usage)
+    return parser
+
+if __name__ == '__main__':
+    parser = get_options()
+    options, args= parser.parse_args()
+    get_table(args,options)
+    calc_eccentricity(args, options)

doc/source/orbit.rst

@@ -346,6 +346,52 @@ behavior
 
 .. image:: images/lp-orbit-integration-EzJz.png
 
+.. _fastchar:
+
+Fast orbit characterization
+-----------------------
+
+It is also possible to use galpy for the fast estimation of orbit parameters via the Staeckel 
+approximation (originally used by `Binney (2012) <http://adsabs.harvard.edu/abs/2012MNRAS.426.1324B>`_
+for the appoximation of actions in axisymmetric potentials), without performing any orbit integration. 
+The method uses the geometry of the orbit tori to estimate the orbit parameters. After initialising 
+an ``Orbit`` instance, this is done (as in the previous section) specifying ``analytic=True`` and 
+selecting ``type='staeckel'`` (default in vX.X of galpy).
+
+>>> o.e(analytic=True, type='staeckel')
+
+if running the above without integrating the orbit, the potential should also be specified
+in the usual way
+
+>>> o.e(analytic=True, type='staeckel', pot=mp)
+
+again, where ``mp`` is the Miyamoto-Nagai potential of :ref:`Introduction:
+Rotation curves <rotcurves>`. This interface automatically computes the necessary Delta 
+parameter based on the initial condition of the ``Orbit`` object.
+
+While this is useful and fast for individual ``Orbit`` objects, it is likely that users will
+want to rapidly evaluate the orbit parameters of large numbers of objects. It is possible
+to perform the orbital parameter estimation above through the :ref:`actionAngle <actionangle>` 
+interface. To do this, we need arrays of the phase-space points ``R``, ``vR``, ``vT``, ``z``, ``vz``, and 
+``phi`` for the objects. We can then optionally estimate the individual Delta parameter 
+for these phase-space points using
+
+>>> from galpy.actionAngle import estimateDeltaStaeckel
+>>> deltas = estimateDeltaStaeckel(mp, R, z, no_median=True)
+
+where ``no_median=True`` specifies that the function return the delta parameter at each given point
+rather than the median of the calculated deltas. The orbit parameters are then calculated by first 
+specifying an ``actionAngle`` instance (with an arbitrary delta parameter), then using the 
+``EccZmaxRperiRap`` method with the data points and the estimated delta array:
+
+>>> aAS = actionAngleStaeckel(pot=mp, delta=0.4)
+>>> e, Zmax, rperi, rap = aAS.EccZmaxRperiRap(R, vR, vT, z, vz, phi, delta=deltas)
+
+Here, delta can alternatively be either a single scalar value for all points, negating the need for
+the estimation at each point. This method is also applicable in the ``actionAngleIsochrone``, 
+``actionAngleSpherical``, and ``actionAngleAdiabatic`` modules. This method returns parameters at
+speeds as fast as 3 microseconds per object. 
+
 
 Accessing the raw orbit
 -----------------------
@@ -507,32 +553,70 @@ a set of thick disk stars. We start by downloading the sample of SDSS
 SEGUE (`2009AJ....137.4377Y
 <http://adsabs.harvard.edu/abs/2009AJ....137.4377Y>`_) thick disk
 stars compiled by Dierickx et al. (`2010arXiv1009.1616D
-<http://adsabs.harvard.edu/abs/2010arXiv1009.1616D>`_) at
-
-http://www.mpia-hd.mpg.de/homes/rix/Data/Dierickx-etal-tab2.txt
+<http://adsabs.harvard.edu/abs/2010arXiv1009.1616D>`_) from CDS at `this 
+link <http://vizier.cfa.harvard.edu/viz-bin/Cat?cat=J%2FApJ%2F725%2FL186&target=http&>`_.
+Downloading the table and the ReadMe will allow you to read in the data using ``astropy.io.ascii``
+like so
+
+>>> from astropy.io import ascii
+>>> dierickx = ascii.read('table2.dat', readme='ReadMe')
+>>> vxvv = numpy.dstack([dierickx['RAdeg'], dierickx['DEdeg'], dierickx['Dist']/1e3, dierickx['pmRA'], dierickx['pmDE'], dierickx['HRV']])[0]
 
 After reading in the data (RA,Dec,distance,pmRA,pmDec,vlos; see above)
 as a vector ``vxvv`` with dimensions [6,ndata] we (a) define the
 potential in which we want to integrate the orbits, and (b) integrate
-each orbit and save its eccentricity (running this for all 30,000-ish
+each orbit and save its eccentricity as calculated analytically following the :ref:`Staeckel 
+approximation method <fastchar>` and by orbit integration (running this for all 30,000-ish
 stars will take about half an hour)
 
+>>> from galpy.actionAngle import UnboundError
+>>> ts= np.linspace(0.,20.,10000)
 >>> lp= LogarithmicHaloPotential(normalize=1.)
->>> ts= nu.linspace(0.,20.,10000)
->>> mye= nu.zeros(ndata)
->>> for ii in range(len(e)):
-...         o= Orbit(vxvv[ii,:],radec=True,vo=220.,ro=8.) #Initialize
-...         o.integrate(ts,lp) #Integrate
-...         mye[ii]= o.e() #Calculate eccentricity
+>>> e_ana = numpy.zeros(len(vxvv))
+>>> e_int = numpy.zeros(len(vxvv))
+>>> for i in range(len(vxvv)):
+...	#calculate analytic e estimate, catch any 'unbound' orbits
+...     try:
+...         orbit = Orbit(vxvv[i], radec=True, vo=220., ro=8.)
+...         e_ana[i] = orbit.e(analytic=True, pot=lp, c=True)
+...     except UnboundError:
+...         #parameters cannot be estimated analytically
+...         e_ana[i] = np.nan
+...     #integrate the orbit and return the numerical e value
+...     orbit.integrate(ts, lp)
+...     e_int[i] = orbit.e(analytic=False)
+
+We then find the following eccentricity distribution (from the numerical eccentricities)
+
+.. image:: images/dierickx-integratedehist.png
+
+The eccentricity calculated by integration in galpy compare well with those
+calculated by Dierickx et al., except for a few objects
 
-We then find the following eccentricity distribution
+.. image:: images/dierickx-integratedee.png
 
-.. image:: images/dierickx-myehist.png
+and the analytical estimates are equally as good:
 
-The eccentricity calculated by galpy compare well with those
-calculated by Dierickx et al., except for a few objects
+.. image:: images/dierickx-analyticee.png
+
+In comparing the analytic and integrated eccentricity estimates - one can see that in this case
+the estimation is almost exact, due to the spherical symmetry of the chosen potential:
+
+.. image:: images/dierickx-integratedeanalytice.png
+
+A script that calculates and plots everything can be downloaded
+:download:`here <examples/dierickx_eccentricities.py>`. To generate the plots just run::
+
+    python dierickx_eccentricities.py ../path/to/folder
+
+specifiying the location you want to put the plots and data.
+
+Alternatively - one can transform the observed coordinates into spherical coordinates and perform 
+the estimations in one batch using the ``actionAngle`` interface, which takes considerably less time:
 
-.. image:: images/dierickx-myee.png
+>>> from galpy import actionAngle
+>>> deltas = actionAngle.estimateDeltaStaeckel(lp, Rphiz[:,0], Rphiz[:,2], no_median=True)
+>>> aAS = actionAngleStaeckel(pot=lp, delta=0.)
+>>> par = aAS.EccZmaxRperiRap(Rphiz[:,0], vRvTvz[:,0], vRvTvz[:,1], Rphiz[:,2], vRvTvz[:,2], Rphiz[:,1], delta=deltas)
 
-The script that calculates and plots everything can be downloaded
-:download:`here <examples/dierickx-edist.py>`.
+The above code calculates the parameters in roughly 100ms on a single core.

galpy/actionAngle_src/actionAngleStaeckel.py

@@ -1045,7 +1045,7 @@ def _vminFindStart(v,E,Lz,I3V,delta,u0,cosh2u0,sinh2u0,
 
 @potential_physical_input
 @physical_conversion('position',pop=True)
-def estimateDeltaStaeckel(pot,R,z):
+def estimateDeltaStaeckel(pot,R,z, no_median=False):
     """
     NAME:
        estimateDeltaStaeckel
@@ -1054,6 +1054,8 @@ def estimateDeltaStaeckel(pot,R,z):
     INPUT:
        pot - Potential instance or list thereof
        R,z- coordinates (if these are arrays, the median estimated delta is returned, i.e., if this is an orbit)
+       no_median - (False) if True, and input is array, return all calculated values of delta (useful for quickly 
+       estimating delta for many phase space points)
     OUTPUT:
        delta
     HISTORY:
@@ -1070,7 +1072,8 @@ def estimateDeltaStaeckel(pot,R,z):
                                                              use_physical=False)))/evaluateRzderivs(pot,R[ii],z[ii],use_physical=False)) for ii in range(len(R))])
         indx= (delta2 < 0.)*(delta2 > -10.**-10.)
         delta2[indx]= 0.
-        delta2= nu.median(delta2[True^nu.isnan(delta2)])
+        if not no_median:
+        	delta2= nu.median(delta2[True^nu.isnan(delta2)])
     else:
         delta2= (z**2.-R**2. #eqn. (9) has a sign error
                  +(3.*R*_evaluatezforces(pot,R,z)