Merge 06cf5f3 into 94f1e81

README.md

@@ -51,7 +51,6 @@ pyuvdata has three major user classes:
 ## Known Issues and Planned Improvements
 * UVData: different multiple spectral windows or multiple sources are not currently supported
 * UVData: testing against miriad package
-* UVData: replacing pyephem with astropy+NOVAS for time and phase calculations
 * UVData: add support for writing CASA measurement sets
 * UVData: phasing is tested to a part in 10^3, and assumes planar array. Improvements are tracked on Issue \#148.
 * UVCal: expand support for calibration solutions: support other formats beyond FITS
@@ -103,14 +102,13 @@ python setup.py build_ext --inplace
 The numpy and astropy versions are important, so be sure to make sure these are up to date before you install.
 
 For anaconda users, we suggest using conda to install astropy, numpy, scipy, and optionally h5py, and
-conda-forge for pyephem and optionally casacore-python and healpy (e.g. ```conda install -c conda-forge pyephem```).
+conda-forge for optionally installing python-casacore and healpy (e.g. ```conda install -c conda-forge python-casacore```).
 
 * numpy >= 1.10
 * scipy
 * astropy >= 1.2
-* pyephem
 * h5py (optional: for reading and writing uvh5 format)
-* casacore-python (optional: for CASA measurement set reading functionality)
+* python-casacore (optional: for CASA measurement set reading functionality)
 * healpy (optional: working with beams in HEALPix formats)
 
 ### For CASA measurement set functionality, install python-casacore

docs/references/phasing.tex

@@ -0,0 +1,70 @@
+\documentclass[11pt, oneside]{article}   	% use "amsart" instead of "article" for AMSLaTeX format
+\usepackage{geometry}                		% See geometry.pdf to learn the layout options. There are lots.
+\geometry{letterpaper}                   		% ... or a4paper or a5paper or ... 
+%\geometry{landscape}                		% Activate for for rotated page geometry
+%\usepackage[parfill]{parskip}    		% Activate to begin paragraphs with an empty line rather than an indent
+\usepackage{graphicx}			
+
+\usepackage{amssymb}
+
+\usepackage{hyperref} 
+\hypersetup{
+    colorlinks = true
+}
+
+\usepackage{cleveref}
+\crefformat{footnote}{#2\footnotemark[#1]#3}
+
+\title{Memo: Phasing in pyuvdata}
+\author{Bryna Hazelton, Miguel Morales and the pyuvdata team}
+\date{June 28, 2018}							% Activate to display a given date or no date
+
+\begin{document}
+\maketitle
+\section{Introduction}
+This memo discusses the phase types that are supported in pyuvdata, the implementation of phasing in pyuvdata (based on the authors' conceptual understanding of phasing in radio interferometry) and the testing that has been done on the pyuvdata phasing code.
+
+\section{Types of phasing supported in pyuvdata}
+pyuvdata supports two primary kinds of phasing, zenith drift phasing (`drift') and data that are phased to a particular RA and Dec defined at a particular epoch (`phased') with very limited support for labelling and handling other or unknown phasing situations (`unknown'). Which kind of phasing any given UVData object has is stored in the `phase\_type' attribute (allowed values: `drift', `phased', `unknown').
+
+Zenith drift phasing is the phasing type that natively comes off of many low-frequency radio telescope correlators (e.g. MWA, PAPER), meaning that the correlations for each time are calculated phased to zenith (i.e. with no baseline-dependent phase applied). The baseline location vectors, the `uvw\_array' attribute, are then just given by the antenna separations in a local East-North-Up (ENU, also called topocentric) coordinate system for all times.
+
+Phased data sets have had rotations applied to the u, v, w coordinates (the ENU baseline location vectors) and baseline-dependent phasing applied to the visibilities to make the visibilities add coherently toward a particular RA and Dec defined at a particular epoch (typically J2000). Phased data sets have thus been corrected to account for the motion of the earth (rotation, precession, nutation and aberration due to the earth's motion around the sun). This is generally required for imaging and to save the data in some file formats (e.g. uvfits).
+
+Moving between these phase types requires functionality to phase zenith drift data to any given RA and Dec (defined at an epoch) and the ability to undo the phasing to get back to zenith drift. Changing between different phase centers can be achieved by first unphasing to drift and then phasing to the new phase center.
+
+\section{Conceptual description of phasing}
+
+
+\section{How phasing is done in pyuvdata}
+pyuvdata's phasing uses astropy to do all the coordinate conversions between earth referenced and solar-system barycenter referenced coordinate systems (called `frames` in astropy). Of particular interest are the ITRS frame (also called ECEF or ITRF) which is fixed to the earth and rotates with the earth, useful for describing antenna and telescope locations, and the ICRS frame which is a non-rotating frame referenced to the solar system barycenter with axes aligned with J2000 axes, useful for defining fixed celestial sources like radio galaxies. One other frame that can be useful is the GCRS frame which is non-rotating with axes aligned with ICRS but geo-center referenced so that it moves with the earth around the solar system barycenter.
+
+Phasing and unphasing in pyuvdata works from the existing `uvw\_array' by default, but can optionally work from the antenna positions (by setting `use\_ant\_pos=True'). By default, phasing and unphasing use the ICRS frame for the phased uvws, but the code can optionally phase to and from GCRS (by setting `phase\_frame=GCRS'). This is not as accurate, but is included to support testing against other codes (e.g. MWA\_tools) and to support unphasing data sets that were originally phased to GCRS (e.g. MWA uvfits files).
+
+To phase zenith drift data, the code performs the following steps:
+\begin{itemize}
+\item{Define the phase center in an astropy ICRS frame. If the epoch is `J2000' or `2000' use the RA and Dec to define a SkyCoord in the ICRS frame. Otherwise define a SkyCoord in the `FK5' frame at the given RA, Dec and epoch and convert it to an ICRS frame for phasing. This means that all phasing is done in the ICRS (J2000) frame, but if the RA and Dec are specified at another epoch, the RA and Dec are converted to ICRS before phasing. Optionally transform the phase center to the GCRS frame if phasing to GCRS.}
+\item{For each time in the `time\_array' attribute: Define the telescope location as a SkyCoord object in the ITRS frame with `obstime' equal to the time.}
+\item{For each time in the `time\_array' attribute: Define the antenna positions or uvws as a SkyCoord object in the ITRS frame with `obstime' equal to the time. For antenna positions, this requires combining the `telescope\_location` attribute (which holds the telescope location in the ITRS frame) and `antenna\_positions` attribute (which holds the antenna positions relative to the telescope location in the ITRS frame). For uvws (which are intrinsically in an ENU frame) this requires converting from ENU to ITRS using pyuvdata.utils.ECEF\_from\_ENU using the `telescope\_location' attribute as well as the `uvw\_array' attribute.}
+\item{For each time in the `time\_array' attribute: Transform the telescope location and the antenna positions or uvws to the ICRS frame (or GCRS if phasing to GCRS).}
+\item{Calculate the relative antenna positions or uvws in the ICRS (or GCRS) frame by subtracting the telescope location in the ICRS (or GCRS) frame.}
+\item{Calculate the phased antenna positions or uvws using the pyuvdata.utils.phase\_uvw function. If using antenna positions, difference these phased positions to calculate the uvws for each baseline. One way to think about what this function does is that is mathematically identical to calculating ENU values, except that it is in the ICRS (or GCRS) frame and the center of the ENU coordinate system is given by the phase center. This makes sense in that if you choose a phase center which is aligned with zenith at the observation time, the changes to the uvws are small and are just due to the changes in the earth spin axis (due to precession and nutation) and the length contraction due to earth's motion around the sun (aberration, excluded if phasing to GCRS).}
+\end{itemize}
+
+To unphase to drift, the steps are essentially followed in reverse order and unit testing guarantees that phasing followed by unphasing gets back to the same answer. 
+
+\section{Testing pyuvdata's phasing}
+
+The primary external testing reference for pyuvdata's phasing code is the phasing code used on the MWA found in the MWA\_tools repository under CONV2UVFITS (and also found in the MWA anoko repository which holds the Cotter code). The MWA code uses the SLA C Library which is proprietary. The authors of this memo and of the pyuvdata code did not attempt to copy or replicate parts of the SLA library, they simply compared calculations done using astropy to calculations performed in the MWA code (which calls the SLA library routines in various ways).
+
+\subsection{Limitations of the MWA reference code}
+By comparing the values calculated by the MWA code at various stages to values calculated using various astropy conversions, the authors of this memo determined that the MWA code is only correcting the u, v, w coordinates for precession and nutation and not for aberration. This is supported by a few lines of evidence, one is that the uvw correction is done using a rotation matrix, which cannot account for length contraction and the other is that the corrected values closely match the values generated using astropy conversions to the GCRS frame, but not to conversions to the ICRS frame (which differ by the inclusion of the effects of earth's motion around the sun, also called aberration).
+
+In addition, there are some discrepancies in the LST calculations done in the MWA code and in astropy. The MWA code uses the mean LST to calculate the RA of zenith in the True Equator and Equinox frame (TEE, also referred to as the current epoch). The mean LST calculated by astropy differs from the mean LST calculated in the MWA code by about 4 arcseconds at the particular time used for testing (mjd = 55780.1). Furthermore, the authors believe that the LST used in this calculation should actually be the apparent LST (which corrects for nutation in addition to precession). The apparent LST calculated by astropy differs from the mean LST calculated in the MWA code by about 12 arceseconds at the same testing time. Note that the LST is not used in pyuvdata's phasing code, instead all the precession and nutation corrections are applied in the astropy frame conversions.
+
+\subsection{External phasing test}
+We received a couple uvfits files from David Kaplan containing identical MWA data phased to two different phase centers. The authors presume that the phasing was done with the MWA phasing code (likely as implemented in Cotter). To test our phasing code, we compared the u, v, w coordinates from each file after unphasing both files and after unphasing and rephasing one of the files to the phase center of the other file. The results show that the uvws in each comparison (drift and rephased) agree to better than 2 cm. We also did the unphase and rephase test with calculating the uvws based on the antenna positions (rather than working with the uvws in the file) and in that test the uvws agree to better than 5 mm. 
+
+
+
+\end{document}  

docs/tutorial.rst

@@ -255,23 +255,23 @@ Phasing/unphasing data
 ::
 
   >>> from pyuvdata import UVData
-  >>> import ephem
+  >>> from astropy.time import Time
   >>> UV = UVData()
   >>> miriad_file = 'pyuvdata/data/new.uvA'
   >>> UV.read_miriad(miriad_file)
   >>> print(UV.phase_type)
   drift
 
   # Phase the data to the zenith at first time step
-  >>> UV.phase_to_time(UV.time_array[0])
+  >>> UV.phase_to_time(Time(UV.time_array[0], format='jd'))
   >>> print(UV.phase_type)
   phased
 
   # Undo phasing to try another phase center
   >>> UV.unphase_to_drift()
 
   # Phase to a specific ra/dec/epoch (in radians)
-  >>> UV.phase(5.23368, 0.710940, ephem.J2000)
+  >>> UV.phase(5.23368, 0.710940, epoch="J2000")
 
 UVData: Plotting
 ------------------