Add full set of celestial conversions; add some documentation

.gitignore

@@ -2,4 +2,5 @@
 *.project
 *.pydevproject
 *.html
+*.pyc
 /plots/

README.md

@@ -1,47 +1,31 @@
+Coordinates Benchmark
+=====================
+
+Here we compare astropy coordinate conversion results against other Python packages.
+
+For now we only compare precision, in the future we also want to compare speed.
+
+Links
+-----
+
+Repository: https://github.com/astropy/coordinates-benchmark
+Results summary: http://mpia.de/~robitaille/astropy/coordinates/summary.html
+
+Mailing list discussions:
+* https://groups.google.com/forum/?fromgroups#!topic/astropy-dev/tbcpMQ1rOcY
+* https://groups.google.com/forum/?fromgroups#!topic/astropy-dev/-ulzSY9qJKk
+* https://groups.google.com/forum/?fromgroups#!searchin/astropy-dev/benchmark
+
+Documentation is in the `docs` folder:
+* Specification.rst
+* Notes.rst
+* Tools.rst
+* HOWTO.rst
+
 History
 -------
 
 - Fixed kapteyn transforms to explicitly specify the epoch of observation for
   the B1950 conversion, and the equinox for the ecliptic conversion
 
 - Explicitly precess astropy.coordinates fk4 result to B1950
-
-Notes
------
-
-From David Berry:
-
-> 1) The differences between kapteyn and (idl,pyast,tpm) on the J2000 to
-> B1950 tests is reduced seriously if you specify the epoch of
-> observation explicitly in kapteyn/convert.py. That is, change
->
-> transform_to(skyout='fk4', tag='b1950')
->
-> to
->
-> transform_to(skyout='fk4,J2000_OBS', tag='b1950')
->
-> This reduces the median difference from 0.217 arc-sec, to 0.0044
-> arc-sec (this must mean that kapteyn is using a different default
-> epoch of observation to the other systems).
->
-> 2) Improved agreement in J2000->ecliptic can be achieved by
-> specifying the equinox explicitly in kapteyn/convert.py. Change:
->
-> transform_to(skyout='ecliptic', tag='ecliptic')
->
-> to
->
-> transform_to(skyout='ecliptic,J2000', tag='ecliptic')
->
-> This causes kapteyn and IDL have zero difference. The kaptyn->pyast
-> median difference goes down from 0.04 to 0.036. The remaining
-> difference between kapteyn and pyast seems to be down to:
->
-> a) DIfferent obliquity models. kaptyn uses the IAU 2000 model, and
-> pyast uses the IAU 2006 model.
-> b) Different precession models. kaptyn uses IAU 1976 and pyast uses IAU 2006.
->
-> Reverting the models used by pyast gives zero difference between
-> kapteyn and pyast for the J2000->Ecliptic test
-

astrolib/convert.py

@@ -0,0 +1,54 @@
+"""
+Coordinate conversions with the astrolib coords package.
+astrolib coords is a Python wrapper around the TPM C library [2].
+
+I couldn't find a repository for TPM, but there is another
+Python wrapper for it called pytpm which has TPM bundled
+and is on github, so you can look at the code there.
+
+[1] http://www.scipy.org/AstroLibCoordsHome
+[2] http://www.sal.wisc.edu/~jwp/astro/tpm/tpm.html
+[3] https://github.com/phn/pytpm/tree/master/src/tpm
+"""
+import numpy as np
+from coords import Position 
+
+SUPPORTED_SYSTEMS = 'fk5 fk4 galactic ecliptic'.split()
+
+def convert(coords, systems):
+
+    if not set(systems.values()).issubset(SUPPORTED_SYSTEMS):
+        return None
+
+    lons = np.zeros_like(coords['lon'])
+    lats = np.zeros_like(coords['lat'])
+
+    for ii, (lon, lat) in enumerate(zip(coords['lon'], coords['lat'])):
+
+        # Create coordinate in input system
+        if systems['in'] == 'fk5':
+            coord = Position((lon, lat))
+        elif systems['in'] == 'fk4':
+            coord = Position((lon, lat), system='celestial', equinox='B1950')
+        elif systems['in'] == 'galactic':
+            coord = Position((lon, lat), system='galactic')
+        elif systems['in'] == 'ecliptic':
+            coord = Position((lon, lat), system='ecliptic')
+        else:
+            raise ValueError()
+
+        # Convert to appropriate output system
+        if systems['out'] == 'fk5':
+            lon, lat = coord.j2000()
+        elif systems['out'] == 'fk4':
+            lon, lat = coord.b1950()
+        elif systems['out'] == 'galactic':
+            lon, lat = coord.galactic()
+        elif systems['out'] == 'ecliptic':
+            lon, lat = coord.ecliptic()
+        else:
+            raise ValueError()
+
+        lons[ii], lats[ii] = lon, lat
+
+    return dict(lon=lons, lat=lats)