diff --git a/Euler/Euler-methods.tex b/Euler/Euler-methods.tex
index 5a2a4ec..4065021 100644
--- a/Euler/Euler-methods.tex
+++ b/Euler/Euler-methods.tex
@@ -253,14 +253,14 @@ \subsection{Piecewise linear}
 expansion is used for $\Rb\Lb \overline{\Delta \qb}$.  In fact, any vector
 can be decomposed in this fashion:
 \begin{equation}
-{\boldsymbol{\chi}} = {\bf I} {\boldsymbol{\chi}} = \Rb \Lb {\boldsymbol{\chi}} = 
+{\boldsymbol{\chi}} = {\bf I} {\boldsymbol{\chi}} = \Rb \Lb {\boldsymbol{\chi}} =
    \sum_\nu (\lb^\enu \cdot {\boldsymbol{\chi}}) \rb^\enu
 \end{equation}
 And then it is easy to see that the above manipulations for $\Ab \Delta \qb$
 can be expressed as:
 \begin{equation}
-\Ab \Delta \qb =  \Ab \sum_\nu (\lb^\enu \cdot \overline{\Delta \qb}) \rb^\enu = 
-  \sum_\nu (\lb^\enu \cdot \overline{\Delta \qb}) \Ab \rb^\enu = 
+\Ab \Delta \qb =  \Ab \sum_\nu (\lb^\enu \cdot \overline{\Delta \qb}) \rb^\enu =
+  \sum_\nu (\lb^\enu \cdot \overline{\Delta \qb}) \Ab \rb^\enu =
   \sum_\nu (\lb^\enu \cdot \overline{\Delta \qb}) \lambda^\enu \rb^\enu
 \end{equation}
 where we used $\Ab \rb^\enu = \lambda^\enu \rb^\enu$.  The quantity $(\lb^\enu
@@ -735,18 +735,18 @@ \subsection{Limiting on characteristic variables}
 limiting on the characteristic variables rather than the primitive
 variables.  The characteristic slopes for the quantity carried by the
 wave $\nu$ can be found from the primitive variables
-as: 
+as:
 %
-\begin{equation} 
-\Delta w^{(\nu)} = \lb^{(\nu)} \cdot \Delta \qb 
-\end{equation} 
+\begin{equation}
+\Delta w^{(\nu)} = \lb^{(\nu)} \cdot \Delta \qb
+\end{equation}
 %
 any limiting would then be done to $\Delta w^{(\nu)}$ and the limited
 primitive variables would be recovered as:
 \begin{equation}
   \overline{\Delta \qb} = \sum_\nu \overline{\Delta w}^{(\nu)}
-  \rb^{(\nu)} 
-\end{equation} 
+  \rb^{(\nu)}
+\end{equation}
 (here we use an overline to indicate limiting).
 
 This is attractive because it is more in the spirit of the linear
@@ -775,7 +775,7 @@ \section{Riemann solvers}
 \includegraphics[width=\linewidth]{multiple_interfaces}
 \caption[Riemann wave structure at each interface]{\label{fig:euler:multiple_interfaces}
   The Riemann wave structure resulting from the separate Riemann problems at each interface.  For each
-  interface, we find the state on the interface and use this to evaluate the flux through the interface.} 
+  interface, we find the state on the interface and use this to evaluate the flux through the interface.}
 \end{figure}
 
 As discussed in \S~\ref{Euler:riemann:solution}, we need to determine
@@ -817,8 +817,8 @@ \section{Riemann solvers}
   approximation does a reasonable job near the intersection and only
   diverges significantly for small $p$ (which is where the solution
   should really be a
-  rarefaction).\\ 
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/riemann-2shock.py}{riemann-2shock.py}}}
+  rarefaction).\\
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/riemann-2shock.py}{riemann-2shock.py}}}
 \end{figure}
 
 An additional approximation concerns rarefactions.  Recall that a
@@ -856,7 +856,7 @@ \section{Conservative update}
 
 Once we have the fluxes, the conservative update is done as
 \begin{equation}
-\Uc^{n+1}_i = \Uc^n_i + \frac{\Delta t}{\Delta x} 
+\Uc^{n+1}_i = \Uc^n_i + \frac{\Delta t}{\Delta x}
    \left ( \Fb_{i-\myhalf}^{n+\myhalf} - \Fb_{i+\myhalf}^{n+\myhalf} \right )
 \end{equation}
 The timestep, $\Delta t$ is determined by the time it takes for the
@@ -1152,7 +1152,7 @@ \section{Source terms}
 
 Note that the source here is cell-centered.  This expansion is
 second-order accurate.  This is the approach outlined in Miller
-\& Colella \cite{millercolella:2002}.  Also notice the 
+\& Colella \cite{millercolella:2002}.  Also notice the
 similarity of this source term to a second-order Euler method
 for integrating ODEs.
 
@@ -1215,22 +1215,22 @@ \section{Source terms}
 \end{align}
 As written, this appears to be an implicit update (since $\Uc^{n+1}$
 depends on $\Hb^{n+1}$), but often, the form of the source terms allows
-you to update the equations in sequence explicitly.  
+you to update the equations in sequence explicitly.
 
 Again, using a constant gravitational acceleration as an example, and
 looking in 1-d for simplicity, $\Uc = (\rho, \rho u, \rho E)^\intercal$
 and $\Hb = (0, \rho g, \rho u g)^\intercal$, so our update sequence is:
 \begin{align}
-\rho^{n+1}_i = \rho^n_i &+ \frac{\Delta t}{\Delta x} 
-   \left [\rho_{i-\myhalf}^{n+\myhalf} u_{i-\myhalf}^{n+\myhalf} - 
+\rho^{n+1}_i = \rho^n_i &+ \frac{\Delta t}{\Delta x}
+   \left [\rho_{i-\myhalf}^{n+\myhalf} u_{i-\myhalf}^{n+\myhalf} -
           \rho_{i+\myhalf}^{n+\myhalf} u_{i+\myhalf}^{n+\myhalf} \right ] \\
-(\rho u)^{n+1}_i = (\rho u)^n_i &+ \frac{\Delta t}{\Delta x} 
-   \left [ \rho_{i-\myhalf}^{n+\myhalf} (u_{i-\myhalf}^{n+\myhalf})^2  - 
+(\rho u)^{n+1}_i = (\rho u)^n_i &+ \frac{\Delta t}{\Delta x}
+   \left [ \rho_{i-\myhalf}^{n+\myhalf} (u_{i-\myhalf}^{n+\myhalf})^2  -
            \rho_{i+\myhalf}^{n+\myhalf} (u_{i+\myhalf}^{n+\myhalf})^2  \right ]
-         + \frac{\Delta t}{\Delta x} \left ( p_{i-\myhalf}^{n+\myhalf} - 
+         + \frac{\Delta t}{\Delta x} \left ( p_{i-\myhalf}^{n+\myhalf} -
                                               p_{i+\myhalf}^{n+\myhalf} \right ) \nonumber \\
        &+ \frac{\Delta t}{2}(\rho^n_i + \rho^{n+1}_i) g \\
-(\rho E)^{n+1}_i = (\rho E)^n_i &+ \frac{\Delta t}{\Delta x} 
+(\rho E)^{n+1}_i = (\rho E)^n_i &+ \frac{\Delta t}{\Delta x}
    \left [ \left (\rho_{i-\myhalf}^{n+\myhalf} E_{i-\myhalf}^{n+\myhalf}  + p_{i-\myhalf}^{n+\myhalf} \right ) u_{i-\myhalf}^{n+\myhalf} - \right . \nonumber \\
    &\phantom{+ \frac{\Delta t}{\Delta x} \left[ \right.}   \left .   \left (\rho_{i+\myhalf}^{n+\myhalf} E_{i+\myhalf}^{n+\myhalf}  + p_{i+\myhalf}^{n+\myhalf} \right ) u_{i+\myhalf}^{n+\myhalf} \right ]
     + \frac{\Delta t}{2} \left [ (\rho u)^n + (\rho u)^{n+1} \right] g
@@ -1282,19 +1282,19 @@ \section{Simple geometries}
 the same ideas as in \S~\ref{euler:sec:sourceterms}.
 
 The conservative update now needs to include the geometry factors.  However,
-it is complicated by the fact that in the momentum equation, the pressure 
-term is a gradient, not a divergence, and therefore has different 
+it is complicated by the fact that in the momentum equation, the pressure
+term is a gradient, not a divergence, and therefore has different
 geometic factors.  The update of the system appears as:
 \begin{align}
-\rho^{n+1}_i = \rho^n_i &+ \frac{\Delta t}{V_i} 
-   \left [A_{i-\myhalf} \rho_{i-\myhalf}^{n+\myhalf} u_{i-\myhalf}^{n+\myhalf} - 
+\rho^{n+1}_i = \rho^n_i &+ \frac{\Delta t}{V_i}
+   \left [A_{i-\myhalf} \rho_{i-\myhalf}^{n+\myhalf} u_{i-\myhalf}^{n+\myhalf} -
           A_{i+\myhalf} \rho_{i+\myhalf}^{n+\myhalf} u_{i+\myhalf}^{n+\myhalf} \right ] \\
-(\rho u)^{n+1}_i = (\rho u)^n_i &+ \frac{\Delta t}{V_i} 
-   \left [ A_{i-\myhalf} \rho_{i-\myhalf}^{n+\myhalf} (u_{i-\myhalf}^{n+\myhalf})^2  - 
+(\rho u)^{n+1}_i = (\rho u)^n_i &+ \frac{\Delta t}{V_i}
+   \left [ A_{i-\myhalf} \rho_{i-\myhalf}^{n+\myhalf} (u_{i-\myhalf}^{n+\myhalf})^2  -
            A_{i+\myhalf} \rho_{i+\myhalf}^{n+\myhalf} (u_{i+\myhalf}^{n+\myhalf})^2  \right ]\nonumber \\
-        & + \frac{\Delta t}{\Delta r} \left ( p_{i-\myhalf}^{n+\myhalf} - 
+        & + \frac{\Delta t}{\Delta r} \left ( p_{i-\myhalf}^{n+\myhalf} -
                                               p_{i+\myhalf}^{n+\myhalf} \right ) \\
-(\rho E)^{n+1}_i = (\rho E)^n_i &+ \frac{\Delta t}{V_i} 
+(\rho E)^{n+1}_i = (\rho E)^n_i &+ \frac{\Delta t}{V_i}
    \left [ A_{i-\myhalf} \left (\rho_{i-\myhalf}^{n+\myhalf} E_{i-\myhalf}^{n+\myhalf}  + p_{i-\myhalf}^{n+\myhalf} \right ) u_{i-\myhalf}^{n+\myhalf} - \right . \nonumber \\
    &\phantom{+ \frac{\Delta t}{V_i} \left[ \right.}   \left .   A_{i+\myhalf} \left (\rho_{i+\myhalf}^{n+\myhalf} E_{i+\myhalf}^{n+\myhalf}  + p_{i+\myhalf}^{n+\myhalf} \right ) u_{i+\myhalf}^{n+\myhalf} \right ]
 \end{align}
@@ -1349,7 +1349,7 @@ \section{Simple geometries}
 Just as with the 1-d spherical case, the pressure term in the momentum
 equation needs to be treated separately from the flux, since it enters
 as a gradient and not a divergence.\footnote{It is common to see the divergence
-term expressed as 
+term expressed as
 \begin{equation}
 \nabla \cdot {\boldsymbol{\phi}} = \frac{1}{r^\alpha} \ddr{(r^\alpha \phi^{(r)})} + \ldots
 \end{equation}
@@ -1541,7 +1541,7 @@ \subsection{Sedov blast wave}
 \begin{figure}[t]
 \centering
 \includegraphics[width=0.75\linewidth]{sedov_compare}
-\caption[2-d cylindrical Sedov problem]{\label{fig:Euler:sedov2d_compare} Angle-average 
+\caption[2-d cylindrical Sedov problem]{\label{fig:Euler:sedov2d_compare} Angle-average
   profile for the 2-d Sedov explosion from Figure~\ref{fig:Euler:sedov2d} shown
   with the analytic solution.  This was constructed using the {\tt sedov\_compare.py}
   script in \pyro.}
@@ -1552,7 +1552,7 @@ \subsection{Advection}
 
 We can run a simple advection test analogous to the tests we used in
 Ch.~\ref{ch:advection}.  However, because we are now doing hydrodynamics,
-we need to suppress the dynamics.  This is accomplished by putting the 
+we need to suppress the dynamics.  This is accomplished by putting the
 profile we want to advect in the density field and then put it in
 pressure equilibrium by adjusting the internal energy.  For example, consider
 a Gaussian profile.  We initialize the density as:
@@ -1598,9 +1598,9 @@ \subsection{Slow moving shock}
 Slow moving (or stationary) shocks can be difficult to model, as
 oscillations can be setup behind the shock (this is discussed a little
 in \cite{colellawoodward:1984,leveque:2002}).  We can produce a slow
-moving shock as a shock tube, and we can use the jump conditions 
+moving shock as a shock tube, and we can use the jump conditions
 across a shock that were derived for the Riemann problem to find the
-conditions to setup a stationary (or slow-moving) shock.  
+conditions to setup a stationary (or slow-moving) shock.
 
 The speed of a right-moving shock was found (see
 Eq.~\ref{eq:euler:shockspeedjump}) as:
@@ -1619,7 +1619,7 @@ \subsection{Slow moving shock}
 (which was the star state when we discussed the Riemann problem) using
 the jump conditions, Eqs.~\ref{eq:euler:shockrhojump} and
 \ref{eq:euler:shockujump}.\footnote{The script
-  \href{https://github.com/zingale/hydro_examples/blob/master/compressible/slow_shock.py}{\tt
+  \href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/slow_shock.py}{\tt
     slow\_shock.py} will find the initial conditions to generate a
   stationary shock.}
 
@@ -1660,7 +1660,7 @@ \subsection{Two-dimensional Riemann problems}
 problems initialize the 4 quadrants of the domain with different
 states, and watch the ensuing evolution.  There are some
 analytic estimates that can be compared to, but also these
-tests can provide a means of assessing the symmetry of 
+tests can provide a means of assessing the symmetry of
 a code in the presence of complex flows.  We use the setup
 corresponding to {\em configuration 3} in that paper (this
 same setup is used in \cite{leveque:1997}).
@@ -1811,12 +1811,12 @@ \subsection{General equation of state}
 
 The above methods were formulated with a constant gamma equation of
 state.  A general equation of state (such as degenerate electrons)
-requires a more complex method.  Most methods are designed to 
+requires a more complex method.  Most methods are designed to
 reduce the need to call a complex equation of state frequently,
-and work by augmenting the vector of primitive variables with 
-additional thermodynamic information.  There are two parts of the 
-adaption to a general equation of state: the interface states and 
-the Riemann problem. 
+and work by augmenting the vector of primitive variables with
+additional thermodynamic information.  There are two parts of the
+adaption to a general equation of state: the interface states and
+the Riemann problem.
 
 \subsubsection{Carrying $\gamma_e$}
 
@@ -1851,8 +1851,8 @@ \subsubsection{Carrying $\gamma_e$}
 in the tracing in the construction of interface states%
 \footnote{A {\sf Jupyter} notebook using {\sf SymPy} that derives these
 eigenvectors is available here:
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/euler-generaleos.ipynb}{euler-generaleos.ipynb}}}.
-If we write 
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/euler-generaleos.ipynb}{euler-generaleos.ipynb}}}.
+If we write
 our system as
 \begin{equation}
 \qb = \left ( \begin{array}{c} \tau \\ u \\ p \\ \gamma_e \end{array} \right )
@@ -1996,6 +1996,3 @@ \subsubsection{Carrying $(\rho e)$}
 explicitly adds the transverse terms found in the multi-dimensional form
 of Eq.~\ref{eq:euler:pgeneral}
 to the normal states of $p$.
-
-
-
diff --git a/Euler/Euler-theory.tex b/Euler/Euler-theory.tex
index 88bad40..8a1e94f 100644
--- a/Euler/Euler-theory.tex
+++ b/Euler/Euler-theory.tex
@@ -10,13 +10,13 @@ \section{Euler equation properties}
 dissipative terms.  However, for astrophysical flows, the scales on
 which these dissipative terms operate are usually much smaller than
 the system of interest (equivalently, Reynolds numbers of
-astrophysical flows are very large).} describe conservation of 
+astrophysical flows are very large).} describe conservation of
 mass, momentum, and energy in the fluid approximation.  Their general
 form, without any source terms, is:  \MarginPar{need to add a discussion of the N-S equations and dimensionless numbers}
 \begin{align}
 \ddt{\rho} + \nabla \cdot (\rho \Ub) &= 0 \\
 \ddt{(\rho \Ub)} + \nabla \cdot (\rho \Ub \Ub) + \nabla p &= 0 \\
-\ddt{(\rho E)} + \nabla \cdot (\rho E \Ub + p \Ub ) &= 0 
+\ddt{(\rho E)} + \nabla \cdot (\rho E \Ub + p \Ub ) &= 0
 \end{align}
 Here $\rho$ is the density, $\Ub$ is the velocity vector, $\Ub =
 u\hat{x} + v\hat{y}$, $p$ is the pressure, and $E$ is the total energy
@@ -80,7 +80,7 @@ \section{Euler equation properties}
 \equiv \rho u$, $\mathcal{E} \equiv \rho E$, and assuming a gamma-law
 EOS%
 \footnote{we can relax this assumption by writing $p = p(\rho, e)$, and then
-taking the derivatives of this as needed:  $\partial p/\partial \rho$, 
+taking the derivatives of this as needed:  $\partial p/\partial \rho$,
 $\partial p/\partial m = \partial p /\partial e|_\rho \partial e/\partial m$,
 and $\partial p/\partial \mathcal{E} = \partial p/\partial e|_\rho \partial e/\partial \mathcal{E}$,
 with $e = (\mathcal{E} - \myhalf m^2/\rho)/\rho$.  But as we'll see, there are
@@ -99,14 +99,14 @@ \section{Euler equation properties}
 The Jacobian%
 \footnote{
 The Jacobian, ${\bf J}$ of a vector $\Fb(\Uc)$ with $\Fb = (f_1, f_2, \ldots, f_n)^\intercal$ and
-$\Uc = (u_1, u_2, \ldots, u_n)^\intercal$ is 
+$\Uc = (u_1, u_2, \ldots, u_n)^\intercal$ is
 \begin{equation*}
 {\bf J} \equiv \frac{\partial \Fb}{\partial \Uc} = \left (
   \begin{array}{cccc}
      \partial f_1/\partial u_1 & \partial f_1/\partial u_2 & \ldots & \partial f_1/\partial u_n \\
      \partial f_2/\partial u_1 & \partial f_2/\partial u_2 & \ldots & \partial f_2/\partial u_n \\
      \vdots                    & \vdots                    & \ddots & \vdots \\
-     \partial f_n/\partial u_1 & \partial f_n/\partial u_2 & \ldots & \partial f_n/\partial u_n 
+     \partial f_n/\partial u_1 & \partial f_n/\partial u_2 & \ldots & \partial f_n/\partial u_n
   \end{array} \right )
 \end{equation*}
 }
@@ -142,7 +142,7 @@ \section{Euler equation properties}
 \end{align}
 }
 \end{exercise}
-Notice that the velocity equation looks like Burgers' equation, and 
+Notice that the velocity equation looks like Burgers' equation, and
 is nonlinear.  This nonlinearity will admit shock and rarefaction
 solutions like we saw with Burgers' equation.
 
@@ -250,7 +250,7 @@ \section{Euler equation properties}
 \end{exercise}
 A {\sf Jupyter} notebook using {\sf SymPy} that derives these
 eigenvectors is available:
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/euler.ipynb}{euler.ipynb}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/euler.ipynb}{euler.ipynb}}
 \footnote{if you don't have {\sf Jupyter}, you can view this rendered online via the github link.}
 
 A final form of the equations is called the {\em characteristic form}.  Here,
@@ -283,7 +283,7 @@ \section{Euler equation properties}
 \end{equation}
 }
 \end{exercise}
-To transform our primitive variable system into the characteristic form, we 
+To transform our primitive variable system into the characteristic form, we
 start by multiplying by $\Lb$:
 \begin{equation}
 \Lb\qb_t + \Lb \Ab \qb_x = 0
@@ -304,7 +304,7 @@ \section{Euler equation properties}
 \begin{align}
 w_t^\evm + \lambda^\evm w_x^\evm &= 0 \\
 w_t^\evz + \lambda^\evz w_x^\evz &= 0 \\
-w_t^\evp + \lambda^\evp w_x^\evp &= 0 
+w_t^\evp + \lambda^\evp w_x^\evp &= 0
 \end{align}
 
 If the system were linear, then the solution to each would simply be
@@ -357,7 +357,7 @@ \section{Euler equation properties}
   the {\em Sod problem}, and in general, a setup with two states
   separated by a discontinuity is called a shock-tube problem.  We see three
   waves carrying changes in the solution. \\
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/riemann-sod.py}{riemann-sod.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/riemann-sod.py}{riemann-sod.py}}}
 \end{figure}
 
 
@@ -444,7 +444,7 @@ \subsection{Rarefactions}
 \end{equation}
 The eigenvalues of $\Ab_s$ are again $u$, $u-c$, and $u+c$.  The right eigenvectors
 are\footnote{Again, see the {\sf SymPy} {\sf Jupyter} notebook:
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/euler.ipynb}{euler.ipynb}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/euler.ipynb}{euler.ipynb}}}
 \begin{align}
 \rb_s^\evm = \left ( \begin{array}{c} 1 \\ -c/\rho \\ 0 \end{array} \right )
 %
@@ -466,9 +466,9 @@ \subsection{Rarefactions}
 wave gives us the conditions (using our original primitive variable system):
 \begin{align}
 \lb^\evp \cdot d\qb &= 0 \\
-\lb^\evz \cdot d\qb &= 0 
+\lb^\evz \cdot d\qb &= 0
 \end{align}
-or 
+or
 \begin{equation}
 \left ( \begin{array}{ccc} 0 & \frac{\rho}{2c} & \frac{1}{2c^2} \end{array} \right)
    \left ( \begin{array}{c} d\rho \\ du \\ dp \end{array} \right ) = 0
@@ -480,9 +480,9 @@ \subsection{Rarefactions}
 Expanding these, we have the system:
 \begin{align}
 du + \frac{1}{\rho c} dp &= 0 \\
-d\rho - \frac{1}{c^2} dp &= 0 
+d\rho - \frac{1}{c^2} dp &= 0
 \end{align}
-Defining the {\em Lagrangian sound speed}, $C \equiv \rho c$, and the 
+Defining the {\em Lagrangian sound speed}, $C \equiv \rho c$, and the
 specific volume, $\tau = 1/\rho$, we can rewrite this system as:
 \begin{equation}
 du = -\frac{dp}{C} , \,\,\, d\tau = -\frac{dp}{C^2} \quad \mbox{across the left wave}
@@ -509,7 +509,7 @@ \subsection{Rarefactions}
 \begin{equation}
 p = K \rho^\gamma
 \end{equation}
-where $K$ is a constant that depends on the entropy, and do the 
+where $K$ is a constant that depends on the entropy, and do the
 integrals analytically.
 
 \begin{exercise}[Riemann invariants for gamma-law gas]
@@ -558,7 +558,7 @@ \subsection{Shocks}
 Rankine-Hugoniot jump conditions.  There will be one condition for
 each of our conservation laws, and together they tell us the speed of
 the shock and how density and pressure jump across
-it.  
+it.
 
 The Rankine-Hugoniot conditions are:
 \begin{equation}
@@ -579,10 +579,10 @@ \subsection{Shocks}
 \begin{align}
 \rho_\star \hat{u}_\star &= \rho_s \hat{u}_s \\
 \rho_\star \hat{u}_\star^2 + p_\star &= \rho_s \hat{u}_s^2 + p_s \\
-\rho_\star \hat{u}_\star e_\star + \frac{1}{2} \rho_\star \hat{u}_\star^3 + \hat{u}_\star p_\star &= 
+\rho_\star \hat{u}_\star e_\star + \frac{1}{2} \rho_\star \hat{u}_\star^3 + \hat{u}_\star p_\star &=
   \rho_s \hat{u}_s e_s + \frac{1}{2} \rho_s \hat{u}_s^3 + \hat{u}_s p_s
 \end{align}
-Our goal is to find how each variable jumps across the shock.  We'll 
+Our goal is to find how each variable jumps across the shock.  We'll
 work this out for the general EOS case, following the ideas in \cite{colellaglaz:1985}.
 
 Starting with the mass flux, we can express:
@@ -641,9 +641,9 @@ \subsection{Shocks}
 [u] = \frac{[p]}{W_R}
 \end{equation}
 
-The last jump condition is for energy.  Since all the terms in the 
+The last jump condition is for energy.  Since all the terms in the
 energy equation are proportional to velocity, there will be no sign
-difference between the left and right shock jump conditions.  
+difference between the left and right shock jump conditions.
 We start by introducing the mass flux:
 \begin{equation}
 W_s e_\star + \frac{1}{2} W_s \hat{u}_\star^2 + \frac{W_s}{\rho_\star} p_\star =
@@ -653,7 +653,7 @@ \subsection{Shocks}
 \begin{equation}
 [e] + \frac{p_\star}{\rho_\star} - \frac{p_s}{\rho_s} + \frac{1}{2} \left ( \hat{u}_\star^2 - \hat{u}_s^2 \right ) = 0
 \end{equation}
-getting rid of the velocities using $\hat{u}_\star^2 = W_s^2/\rho_\star^2$ and 
+getting rid of the velocities using $\hat{u}_\star^2 = W_s^2/\rho_\star^2$ and
 $\hat{u}_s^2 = W_s^2/\rho_s^2$, we have:
 \begin{equation}
 [e] + \frac{p_\star}{\rho_\star} - \frac{p_s}{\rho_s} + \frac{1}{2} W_s^2 \left ( \frac{1}{\rho_\star^2} - \frac{1}{\rho_s^2} \right ) = 0
@@ -682,7 +682,7 @@ \subsection{Shocks}
   \item use Newton's method (or another technique) with $[e] = -\bar{p} [\tau]$
      to find a correction to $\rho_\star$
   \end{enumerate}
-\item compute 
+\item compute
   \begin{equation}
     \frac{1}{W_s^2} = - \frac{[\tau]}{[p]}
   \end{equation}
@@ -708,7 +708,7 @@ \subsection{Shocks}
 
 \begin{exercise}[Shock jump conditions for $\gamma$-law EOS]
 {
-Introducing 
+Introducing
 \begin{equation}
 e = \frac{p}{\rho} \frac{1}{\gamma -1}
 \end{equation}
@@ -744,7 +744,7 @@ \subsection{Shocks}
 \end{equation}
 with the `$-$' or the left shock and the `+' for the right shock.
 \end{exercise}
-  
+
 \subsection{Finding the Star State}
 
 \label{Euler:riemann:starstate}
@@ -757,12 +757,12 @@ \subsection{Finding the Star State}
 \begin{equation}
 u_{\star,L}(p) = \begin{cases}
    u_{\star,L}^\mathrm{shock}(p) & p > p_L \\
-   u_{\star,L}^\mathrm{rare}(p) & p \le p_L 
+   u_{\star,L}^\mathrm{rare}(p) & p \le p_L
    \end{cases}
 \qquad
 u_{\star,r}(p) = \begin{cases}
    u_{\star,R}^\mathrm{shock}(p) & p > p_R \\
-   u_{\star,R}^\mathrm{rare}(p) & p \le p_R 
+   u_{\star,R}^\mathrm{rare}(p) & p \le p_R
    \end{cases}
 \end{equation}
 The solution
@@ -811,7 +811,7 @@ \subsection{Finding the Star State}
 curves are shown.  The solution to the Riemann problem is the point
 where the curves intersect.  The line style of the curve indicates
 where we are a shock or rarefaction.  Where $p > p_s$, where $s \in {L,R}$, we have a shock.\\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/riemann-phase.py}{riemann-phase.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/riemann-phase.py}{riemann-phase.py}}}
 \end{figure}
 
 
@@ -832,7 +832,7 @@ \subsection{Complete Solution}
 \includegraphics[width=0.49\linewidth]{riemann_waves_ifc_R}
 \includegraphics[width=0.49\linewidth]{riemann_waves_ifc_Rstar} \\
 \includegraphics[width=0.49\linewidth]{riemann_waves_ifc_Lstar}
-\includegraphics[width=0.49\linewidth]{riemann_waves_ifc_L} 
+\includegraphics[width=0.49\linewidth]{riemann_waves_ifc_L}
 \caption[Wave configuration for the Riemann problem]
         {\label{fig:euler:riemann_sample} An illustration of the 3
           waves emanating from an initial discontinuity (the origin of
@@ -849,7 +849,7 @@ \subsection{Complete Solution}
 configurations of the waves.  Note that the middle wave is always a
 contact but the left (1) and right (3) waves can be either a shock or
 rarefaction.  To find out which region the interface falls in, we simply
-look at the speeds.  
+look at the speeds.
 
 The first speed to consider is the contact wave, that has a speed of
 simply $S_c = u_\star$.  If $S_c < \xi$, then we are choosing between
@@ -865,7 +865,7 @@ \subsection{Complete Solution}
 between just two regions, either $L$-$L_\star$ or $R$-$R_\star$.  This
 leaves just a single wave to consider: the right wave for cases a and
 b; and the left wave for cases c and d.  We need the wave speed for
-this---it will depend on whether it is a shock or a rarefaction.  
+this---it will depend on whether it is a shock or a rarefaction.
 
 For a shock, the wave speed is given by
 Eq.~\ref{eq:euler-theory:shockgeneral}.  We do the same procedure as
@@ -882,13 +882,13 @@ \subsection{Complete Solution}
 corresponding wave speeds are:
 \begin{itemize}
 \item left (1) rarefaction:
-  \begin{itemize} 
+  \begin{itemize}
   \item $\lambda_\mathrm{head} = u_L - c_L$
   \item $\lambda_\mathrm{tail} = u_\star - c_\star$
   \end{itemize}
-  
+
 \item right (3) rarefaction:
-  \begin{itemize} 
+  \begin{itemize}
   \item $\lambda_\mathrm{head} = u_R + c_R$
   \item $\lambda_\mathrm{tail} = u_\star + c_\star$
   \end{itemize}
@@ -921,15 +921,15 @@ \subsection{Complete Solution}
 \includegraphics[width=0.6\linewidth]{rarefaction_center} \\
 \includegraphics[width=0.6\linewidth]{rarefaction_right}
 \caption[Rarefaction configuration for the Riemann problem]
-        {\label{fig:euler:rarefaction_sample} An illustration of the 
+        {\label{fig:euler:rarefaction_sample} An illustration of the
           structure of the left rarefaction wave, for the case
           where the contact and right wave are to the right of the
           interface.  Here, we are choosing between states $L$ and $L_\star$.
-          In (a), both the head and tail of the rarefaction are to the 
+          In (a), both the head and tail of the rarefaction are to the
           left of the interface, so state $L_\star$ is on the interface.
-          In (c), both the head and tail of the rarefaction are to the 
+          In (c), both the head and tail of the rarefaction are to the
           right of the interface, so state $L$ is on the interface.
-          Inbetween, case (b), shows a rarefaction that spans the 
+          Inbetween, case (b), shows a rarefaction that spans the
           interface.  For this case, we need to integrate the Riemann
           invariants to find the state on the interface.}
 \end{figure}
@@ -962,7 +962,7 @@ \section{Other thermodynamic equations}
 
 Notice that internal energy, $(\rho e)$ is not a conserved quantity
 (in particular, the $p\nabla \cdot \Ub$ term is not in conservative
-form).  
+form).
 
 Another energy-like quantity that we can consider is specific enthalpy,
 \begin{equation}
@@ -1192,4 +1192,3 @@ \subsection{Eigensystem with temperature}
 The eigensystem will change with this addition.  We don't explore this here at this
 time.
 \fi
-
diff --git a/advection/advection-basics.tex b/advection/advection-basics.tex
index eb43fcf..52a8eae 100644
--- a/advection/advection-basics.tex
+++ b/advection/advection-basics.tex
@@ -75,7 +75,7 @@ \section{First-order advection in 1-d and finite-differences}
 \label{eq:fo}
 \end{equation}
 This is an {\em explicit} method, since the new solution, $a_i^{n+1}$,
-depends only on information at the old time level, $n$.  
+depends only on information at the old time level, $n$.
 
 Finally, we also need to specify a boundary condition for this.  Our
 choice of spatial derivative is one-sided---it uses information from
@@ -112,10 +112,10 @@ \section{First-order advection in 1-d and finite-differences}
 \Delta t \le \frac{\Delta x}{u} \enskip .
 \end{equation}
 This is called the {\em Courant-Friedrichs-Lewy} or {\em CFL}
-condition.  A dimensionless quantity called the {\em CFL number} is 
-defined as 
+condition.  A dimensionless quantity called the {\em CFL number} is
+defined as
 \begin{equation}
-\cfl = \frac{\Delta t u }{\Delta x} 
+\cfl = \frac{\Delta t u }{\Delta x}
 \end{equation}
 Stability requires $\cfl \le 1$.
 %
@@ -133,7 +133,7 @@ \section{First-order advection in 1-d and finite-differences}
 \end{exercise}
 
 Keep in mind that, in general, we will be solving a non-linear
-system of equations, so it is not possible to run with $\cfl=1$, 
+system of equations, so it is not possible to run with $\cfl=1$,
 since $u$ (and therefore $\cfl$) will change from zone to zone.
 Instead, one looks at the most restrictive timestep over all the
 zones and uses that for the entire system.
@@ -174,10 +174,10 @@ \section{First-order advection in 1-d and finite-differences}
 {\label{fig:fdadvect} Finite-difference solution to the first-order
 finite-difference upwind method for advection, using 65 points and
 a variety of CFL numbers. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/fdadvect.py}{fdadvect.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/fdadvect.py}{fdadvect.py}}}
 \end{figure}
 %
-This method is first-order accurate.  
+This method is first-order accurate.
 
 
 Ultimately we will want higher-order accurate methods.  The most obvious
@@ -201,7 +201,7 @@ \section{First-order advection in 1-d and finite-differences}
 finite-difference method for advection using 65 points, modeling for
 only 1/10$^\mathrm{th}$ of a period.  The panel of the left is
 with $\cfl = 0.1$ and the panel on the right is $\cfl = 0.5$.\\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/fdadvect.py}{fdadvect.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/fdadvect.py}{fdadvect.py}}}
 \end{figure}
 
 You will find that no matter what value of $\cfl$ you choose with the
@@ -216,7 +216,7 @@ \section{Stability}
 
 \MarginPar{add a discussion of domain of dependence}
 
-The classic method for understanding stability is to consider the growth 
+The classic method for understanding stability is to consider the growth
 of a single Fourier mode in our discretization.  That is, substitute in
 $a_i^n = A^n e^{Ii\theta}$, where $I = \sqrt{-1}$\footnote{our use of $i$ and $j$ as spatial indices presents an unfortunate clash of notation here, hence the use of $I$ for the imaginary unit}, and $\theta$ represents a
 phase.  A method is stable if $|A^{n+1}/A^n| \le 1$.  FTCS appears as:
@@ -237,14 +237,14 @@ \section{Stability}
 We see that there is no value of $\cfl$ that can make the method stable
 ($|A^{n+1}/A^n| > 1$ always).  Doing the same analysis for
 Eq.~\ref{eq:fo} would show that the upwind method is stable for $0\le
-\cfl \le 1$. 
+\cfl \le 1$.
 
 \begin{exercise}[Stability of the upwind method]
 {Using the above stability analysis, considering the amplitude of a
   single Fourier mode, show that the growth of a mode for the upwind
   method (Eq.~\ref{eq:fo}) is:
   \begin{equation}
-    \left | \frac{A^{n+1}}{A^n} \right |^2 = 
+    \left | \frac{A^{n+1}}{A^n} \right |^2 =
            1 - 2\cfl(1-\cfl)(1-\cos\theta)
   \end{equation}
   and stability requires $2\cfl(1-\cfl) \ge 0$ or $0 \le \cfl \le 1$.
@@ -270,8 +270,8 @@ \section{Stability}
 and replace $a_i^{n+1}$ with a Taylor expansion in time, and replace
 $a_{i-1}^n$ with a Taylor expansion in space, keeping terms to
 $O(\Delta t^3)$ and $O(\Delta x^3)$.  Replacing $\partial a/\partial t$
-with $-u \partial a/ \partial x$ in the higher-order terms, show 
-that our difference equation more closely corresponds to 
+with $-u \partial a/ \partial x$ in the higher-order terms, show
+that our difference equation more closely corresponds to
 \begin{eqnarray}
 \label{eq:advect_trunc_analysis}
 a_t + u a_x &=& \frac{u \Delta x}{2} \left ( 1 - \frac{\Delta t u}{\Delta x} \right ) \frac{\partial^2 a}{\partial x^2} \\
@@ -341,7 +341,7 @@ \section{Implicit-in-time}
 identical, so we only need to update one of these.  Taking $a_0^{n+1} = a_{N-1}^{n+1}$, our system in matrix form appears as:
 \begin{equation}
 \renewcommand{\arraystretch}{1.5}
-\left ( \begin{array}{ccccccc} 
+\left ( \begin{array}{ccccccc}
 1+\cfl & & & & & & -\cfl \\
 -\cfl  & 1+\cfl &  \\
  &  -\cfl & 1+\cfl &  \\
@@ -378,7 +378,7 @@ \section{Implicit-in-time}
 expensive than explicit methods.  However, stability analysis would show
 that this implicit discretization is stable for any choice of $\cfl$. (But
 one must not confuse stability with accuracy---the most accurate solutions
-with this method will still have a small $\cfl$).  Also note that the form of 
+with this method will still have a small $\cfl$).  Also note that the form of
 the matrix will change depending on the choice of boundary conditions.
 Figure~\ref{fig:fdadvect-implicit} shows the result of solving this
 implicit system.
@@ -391,7 +391,7 @@ \section{Implicit-in-time}
 {\label{fig:fdadvect-implicit} Finite-difference solution to the implicit first-order
 finite-difference upwind method for advection, using 65 points and
 a variety of CFL numbers. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/fdadvect_implicit.py}{fdadvect\_implicit.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/fdadvect_implicit.py}{fdadvect\_implicit.py}}}
 \end{figure}
 
 \begin{exercise}[Implicit advection]
@@ -412,10 +412,10 @@ \section{Eulerian vs.\ Lagrangian frames}
    \frac{dx}{dt}
 \end{equation}
 So the value of this derivative depends on the path, $x(t)$, that we choose
-to follow.  
+to follow.
 
 Consider an observer who is stationary.  They will watch the flow move
-past them, so $dx/dt = 0$, and $da/dt = \partial a/\partial t$.  
+past them, so $dx/dt = 0$, and $da/dt = \partial a/\partial t$.
 This fixed frame is called {\em Eulerian frame}.
 
 Instead imagine an observer who moves with the flow, at the velocity $u$.
@@ -449,7 +449,7 @@ \section{Errors and convergence rate}
   definition of the norms} as:
 \begin{equation}
 \epsilon^\mathrm{abs} = \| a^\mathrm{final} - a^\mathrm{init} \|_2 \equiv
-   \left [ \frac{1}{N} \sum_{i=1}^N 
+   \left [ \frac{1}{N} \sum_{i=1}^N
    ( a_i^\mathrm{final} - a_i^\mathrm{init} )^2
   \right ]^{\myhalf}
 \end{equation}
diff --git a/advection/advection-higherorder.tex b/advection/advection-higherorder.tex
index 1d76aa7..8900099 100644
--- a/advection/advection-higherorder.tex
+++ b/advection/advection-higherorder.tex
@@ -4,8 +4,8 @@ \section{Advection and the finite-volume method}
 
 \label{ch:adv:sndorder}
 
-In these notes, we will typically use a {\em finite-volume} discretization.  Here we 
-explore this method for the 
+In these notes, we will typically use a {\em finite-volume} discretization.  Here we
+explore this method for the
 advection equation.  First we rewrite the advection equation in {\em
   conservation form}:
 \begin{equation}
@@ -13,7 +13,7 @@ \section{Advection and the finite-volume method}
 \label{eq:advect-cons}
 \end{equation}
 where $f(a) = ua$ is the flux of the quantity $a$.  In conservation form,
-the time derivative of a quantity is related to the divergence of 
+the time derivative of a quantity is related to the divergence of
 its flux.
 
 % figure created with figures/advection/fv-ghost.py
@@ -36,9 +36,9 @@ \section{Advection and the finite-volume method}
 from ${x_{i-\myhalf}}$ to ${x_{i+\myhalf}}$, normalizing by the zone width,
 $\Delta x$:
 \begin{align}
-\frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} a_t \, dx &= 
+\frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} a_t \, dx &=
    - \frac{1}{\Delta x} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \frac{\partial}{\partial x} f(a) \, dx \\
-\frac{\partial}{\partial t} a_i &= 
+\frac{\partial}{\partial t} a_i &=
    - \frac{1}{\Delta x} \left \{ \left [f(a) \right ]_{i+\myhalf} - \left [f(a) \right ]_{i-\myhalf} \right \} \label{adv:eq:fvadv}
 \end{align}
 This is an evolution equation for the zone-average of $a$, and shows
@@ -209,8 +209,8 @@ \section{Second-order predictor-corrector scheme}
 {\label{fig:fvadvect} Second-order finite volume advection showing the
 result of advecting a tophat profile through five periods with both
 unlimited and limited slopes.  This calculation used 64 zones and
-$\cfl=0.7$. \\ 
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
+$\cfl=0.7$. \\
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
 \end{figure}
 
 \subsection{Limiting}
@@ -235,9 +235,9 @@ \subsection{Limiting}
   \frac{a_i - a_{i-1}}{\Delta x}, \frac{a_{i+1} - a_i}{\Delta x} \right )
 \end{equation}
 instead of Eq.~\ref{eq:slopecentered}.
-with 
+with
 \begin{equation}
-\mathtt{minmod}(a,b) = \left \{ 
+\mathtt{minmod}(a,b) = \left \{
     \begin{array}{ll}
     a & \mathit{if~} |a| < |b| \mathrm{~and~} a\cdot b > 0 \\
     b & \mathit{if~} |b| < |a| \mathrm{~and~} a\cdot b > 0 \\
@@ -277,7 +277,7 @@ \subsection{Limiting}
 \includegraphics[width=0.75\linewidth]{rea-nolimit-start_003} \\
 \includegraphics[width=0.75\linewidth]{rea-nolimit-start_004} \\
 \includegraphics[width=0.75\linewidth]{rea-nolimit-start_005} \\
-\includegraphics[width=0.75\linewidth]{rea-nolimit-start_006} 
+\includegraphics[width=0.75\linewidth]{rea-nolimit-start_006}
 %\includegraphics[width=0.55\linewidth]{rea-nolimit-start_007} \\
 %\includegraphics[width=0.55\linewidth]{rea-nolimit-start_008} \\
 \caption[The effect of no limiting on initially discontinuous data]{\label{fig:limitingex}Initially discontinuous data evolved for several steps with
@@ -310,7 +310,7 @@ \subsection{Limiting}
 \end{equation}
 Then the limited difference is
 \begin{equation}
-\left . \frac{\partial a}{\partial x} \right |_i = 
+\left . \frac{\partial a}{\partial x} \right |_i =
  \left \{
 \begin{array}{ll}
 \min \left [ \frac{| a_{i+1} - a_{i-1} |}{2 \Delta x},
@@ -324,7 +324,7 @@ \subsection{Limiting}
 %
 Note that a slightly different form of this limiter is presented in
 \cite{leveque:2002}, where all quantities are in a {\tt minmod}, which appears to limit a bit less.
-This is second-order accurate for smooth flows.  
+This is second-order accurate for smooth flows.
 
 The main goal of a limiter is to reduce the slope near extrema.
 Figure~\ref{fig:limit} shows a finite-volume grid with the original
@@ -372,18 +372,18 @@ \subsection{Limiting}
 \subsection{Reconstruct-evolve-average}
 
 Another way to think about these methods is as a reconstruction,
-evolve, and average (R-E-A) process (see Figure~\ref{fig:rea}).  
+evolve, and average (R-E-A) process (see Figure~\ref{fig:rea}).
 
 We can write the conservative update as:
 \begin{align}
-a_i^{n+1} &= a_i^n + \frac{\Delta t}{\Delta x} 
+a_i^{n+1} &= a_i^n + \frac{\Delta t}{\Delta x}
     (u a^{n+\myhalf}_{i-\myhalf} - u a^{n+\myhalf}_{i+\myhalf} ) \\
-          &= a_i^n + \cfl (a^{n+\myhalf}_{i-\myhalf} - a^{n+\myhalf}_{i+\myhalf} ) 
+          &= a_i^n + \cfl (a^{n+\myhalf}_{i-\myhalf} - a^{n+\myhalf}_{i+\myhalf} )
 \end{align}
 If we take $u > 0$, then the Riemann problem will always choose the
 left state, so we can write this as:
 \begin{equation}
-a_i^{n+1} = a_i^n + 
+a_i^{n+1} = a_i^n +
      \cfl \biggl [\underbrace{\left (a_{i-1}^n + \frac{1}{2} (1-\cfl) \Delta a_{i-1} \right )}_{a_{i-\myhalf,L}} -
               \underbrace{\left (a_{i}^n + \frac{1}{2} (1-\cfl) \Delta a_{i} \right )}_{a_{i+\myhalf,L}}
        \biggr ] \label{eq:rea_orig}
@@ -394,20 +394,20 @@ \subsection{Reconstruct-evolve-average}
 \begin{equation}
 a_i(x) = a_i + \frac{\Delta a_i}{\Delta x} (x - x_i)
 \end{equation}
-Consider zone $i$.  
+Consider zone $i$.
 If we are advecting with a CFL number $\cfl$, then that means that the fraction
 $\cfl$ of the zone immediately to the left of the $i-\myhalf$ interface will advect
 into zone $i$ over the timestep.  And only the fraction $1-\cfl$ in zone $i$
-immediately to the right of the interface will stay in that zone.  This 
-is indicated by the shaded regions in Figure~\ref{fig:rea}. 
+immediately to the right of the interface will stay in that zone.  This
+is indicated by the shaded regions in Figure~\ref{fig:rea}.
 
 The average of the quantity $a$ from zone $i-1$ that will advect into
-zone $i$ is 
+zone $i$ is
 \begin{eqnarray}
-\mathcal{I}_< &=& \frac{1}{\cfl \Delta x} 
+\mathcal{I}_< &=& \frac{1}{\cfl \Delta x}
    \int_{x_{i-\myhalf} - \cfl\Delta x}^{x_{i-\myhalf}} a_{i-1}(x) dx \\
 %
- &=& \frac{1}{\cfl \Delta x} 
+ &=& \frac{1}{\cfl \Delta x}
    \int_{x_{i-\myhalf} - \cfl\Delta x}^{x_{i-\myhalf}}
         \left [ a_{i-1} + \frac{\Delta a_{i-1}}{\Delta x} (x - x_{i-1} ) \right ] dx  \\
  &=& a_{i-1} + \frac{1}{2} \Delta a_{i-1} (1-\cfl)
@@ -416,16 +416,16 @@ \subsection{Reconstruct-evolve-average}
 And the average of the quantity $a$ in zone $i$ that will remain in zone $i$
 is
 \begin{eqnarray}
-\mathcal{I}_> &=& \frac{1}{(1-\cfl) \Delta x} 
+\mathcal{I}_> &=& \frac{1}{(1-\cfl) \Delta x}
    \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl) \Delta x} a_{i}(x) dx \\
 %
- &=& \frac{1}{(1-\cfl) \Delta x} 
-   \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl)\Delta x} 
+ &=& \frac{1}{(1-\cfl) \Delta x}
+   \int_{x_{i-\myhalf}}^{x_{i-\myhalf} + (1-\cfl)\Delta x}
         \left [ a_{i} + \frac{\Delta a_{i}}{\Delta x} (x - x_{i} ) \right ] dx  \\
  &=& a_{i} - \frac{1}{2} \Delta a_{i} \cfl
 \end{eqnarray}
 
-The final part of the R-E-A procedure is to average the over the 
+The final part of the R-E-A procedure is to average the over the
 advected profiles in the new cell.  The weighted average of the
 amount brought in from the left of the interface and that that remains
 in the cell is
@@ -435,10 +435,10 @@ \subsection{Reconstruct-evolve-average}
    + (1-\cfl) \left [ a^n_i - \frac{1}{2} \Delta a_i \cfl \right ] \\
           &= a^n_i + \cfl \left [a^n_{i-1} + \frac{1}{2}(1 - \cfl) \Delta a_{i-1} \right ]
            - \cfl \left [ a^n_i + \frac{1}{2} (1-\cfl) \Delta a_i \right ]
-\end{align}          
+\end{align}
 This is identical to Eq.~\ref{eq:rea_orig}.  This demonstrates that the
 R-E-A procedure is equivalent to our reconstruction, prediction of the
-interface states, solving the Riemann problem, and doing the 
+interface states, solving the Riemann problem, and doing the
 conservative flux update.
 
 % this figure sequence is created by figures/advection/rea.py
@@ -484,18 +484,18 @@ \subsection{Reconstruct-evolve-average}
         {\label{fig:advnorm} Convergence for the second-order
           finite-volume method with no limiting, MC, and minmod limiting advecting a
           Gaussian initial profile with $\cfl = 0.8$ for 5 periods. \\
-        \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
+        \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
 \end{figure}
 
 \begin{figure}[ht]
 \centering
 \includegraphics[width=0.8\linewidth]{fv-gaussian-limiters} \\[1em]
-\includegraphics[width=0.8\linewidth]{fv-tophat-limiters} 
+\includegraphics[width=0.8\linewidth]{fv-tophat-limiters}
 \caption[Effect of different limiters on evolution]
         {\label{fig:limiter_panel} The effect of different limiters on the
           evolution of a Gaussian initial profile (top) and a tophat initial
           profile (bottom), using $\cfl = 0.8$ and 5 periods. \\
-        \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
+        \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/advection.py}{advection.py}}}
 \end{figure}
 
 As seen in figure~\ref{fig:advnorm}, the choice of limiters can have a great
@@ -566,7 +566,7 @@ \section{Multi-dimensional advection}
 \caption[A 2-d grid with zone-centered indexes]{\label{fig:2dgrid} A
   simple 2-d grid with the zone-centered indexes.  The $\times$s mark
   the left and right interface states at the upper edge of the $i,j$ zone in each
-  coordinate direction.  For a finite-volume discretization, $a_{i,j}$ 
+  coordinate direction.  For a finite-volume discretization, $a_{i,j}$
   represents the average of $a(x,y)$ over the zone.}
 \end{figure}
 
@@ -581,21 +581,21 @@ \section{Multi-dimensional advection}
 define the average of $a$ in a zone by integrating it over the
 volume:
 \begin{equation}
-a_{i,j} = \frac{1}{\Delta x \Delta y} 
-   \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} 
+a_{i,j} = \frac{1}{\Delta x \Delta y}
+   \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}}
    a(x,y,t) \, dx \, dy
 \end{equation}
 Integrating Eq.~\ref{eq:advect2d-cons} over $x$ and $y$, we have:
 \begin{align}
-\frac{1}{\Delta x \Delta y} 
-  \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} 
-  \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} a_t \, dx \, dy =  
+\frac{1}{\Delta x \Delta y}
+  \int_{x_{i-\myhalf}}^{x_{i+\myhalf}}
+  \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} a_t \, dx \, dy =
   &- \frac{1}{\Delta x \Delta y}
        \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}}
       (u a)_x \, dx \, dy \nonumber \\
   &- \frac{1}{\Delta x \Delta y}
        \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}}
-      (v a)_y \, dx \, dy 
+      (v a)_y \, dx \, dy
 \end{align}
 or using the divergence theorem,
 \begin{align}
@@ -610,7 +610,7 @@ \section{Multi-dimensional advection}
 Now we integrate over time---the left side of our expression becomes
 just the difference between the new and old state.
 \begin{align}
- a_{i,j}^{n+1} - a_{i,j}^n = 
+ a_{i,j}^{n+1} - a_{i,j}^n =
   &- \frac{1}{\Delta x\Delta y} \int_{t^n}^{t^{n+1}} \int_{y_{j-\myhalf}}^{y_{j+\myhalf}}
      \left \{ (u a)_{i+\myhalf,j} - (u a)_{i-\myhalf,j} \right \} dy dt \nonumber \\
   &- \frac{1}{\Delta x\Delta y} \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}}
@@ -628,13 +628,13 @@ \section{Multi-dimensional advection}
 \item y-face
 \begin{equation}
 \langle (va)_{i,j+\myhalf}\rangle_{(t)} = \frac{1}{\Delta x \Delta t}
-    \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx dt 
+    \int_{t^n}^{t^{n+1}} \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx dt
 \end{equation}
 \end{itemize}
 where $\langle . \rangle_{(t)}$ denotes the time-average over the face.
 For a second-order accurate method in time, we replace the average in
 time with the flux at the midpoint in time and the average over the face
-with the flux at the center of the face: 
+with the flux at the center of the face:
 \begin{equation}
 \langle (ua)_{i+\myhalf,j} \rangle_{(t)} \approx (ua)_{i+\myhalf,j}^{n+\myhalf}
 \end{equation}
@@ -647,7 +647,7 @@ \section{Multi-dimensional advection}
 
 For the advection problem, since $u$ and $v$ are constant, we need
 to find the interface states of $a$ alone.
-There are two methods for computing these interface states, 
+There are two methods for computing these interface states,
 $a_{i\pm\myhalf,j}^{n+\myhalf}$ on $x$-interfaces and  $a_{i,j\pm\myhalf}^{n+\myhalf}$ on $y$-interfaces:
 dimensionally split and unsplit.  Dimensionally split methods are
 easier to code, since each dimension is operated on independent of the
@@ -691,17 +691,17 @@ \subsection{Dimensionally split}
 case (Eq.~\ref{eq:statel} and \ref{eq:stater}).  For example, the
 $a_{i+\myhalf,j,L}^{n+\myhalf}$ state is:
 \begin{eqnarray}
-a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + 
-  \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + 
-  \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + 
+a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n +
+  \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} +
+  \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} +
   \mathcal{O}(\Delta x^2) + \mathcal{O}(\Delta t^2) \nonumber \\
- &=& a_{i,j}^n + 
-   \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + 
-   \frac{\Delta t}{2} \left ( 
+ &=& a_{i,j}^n +
+   \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} +
+   \frac{\Delta t}{2} \left (
    - u \left .\frac{\partial a}{\partial x} \right |_{i,j} \right
    ) + \ldots \nonumber \\
-    &=& a_{i,j}^n + 
-   \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) 
+    &=& a_{i,j}^n +
+   \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right )
    \left .\frac{\partial a}{\partial x} \right |_{i,j} +
    \ldots \label{eq:statels}
 \end{eqnarray}
@@ -736,18 +736,18 @@ \subsection{Unsplit multi-dimensional advection}
 
 The construction of the $a_{i+\myhalf,j,L}^{n+\myhalf}$ interface state appears as
 \begin{eqnarray}
-a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n + 
-  \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + 
-  \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} + 
+a_{i+\myhalf,j,L}^{n+\myhalf} &=& a_{i,j}^n +
+  \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} +
+  \frac{\Delta t}{2} \left .\frac{\partial a}{\partial t} \right |_{i,j} +
   \mathcal{O}(\Delta x^2) + \mathcal{O}(\Delta t^2) \nonumber \\
- &=& a_{i,j}^n + 
-   \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} + 
-   \frac{\Delta t}{2} \left ( 
-   - u \left .\frac{\partial a}{\partial x} \right |_{i,j} 
+ &=& a_{i,j}^n +
+   \frac{\Delta x}{2} \left .\frac{\partial a}{\partial x} \right |_{i,j} +
+   \frac{\Delta t}{2} \left (
+   - u \left .\frac{\partial a}{\partial x} \right |_{i,j}
    - v \left .\frac{\partial a}{\partial y} \right |_{i,j} \right
    ) + \ldots \nonumber \\
-    &=& \underbrace{a_{i,j}^n + 
-   \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) 
+    &=& \underbrace{a_{i,j}^n +
+   \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right )
    \left .\frac{\partial a}{\partial x} \right |_{i,j}}_{\hat{a}_{i+\myhalf,j,L}^{n+\myhalf}} \underbrace{-
    \frac{\Delta t}{2} v \left .\frac{\partial a}{\partial y} \right |_{i,j}}_{\mathrm{``transverse~flux~difference''}} +
    \ldots \label{eq:statelu}
@@ -756,7 +756,7 @@ \subsection{Unsplit multi-dimensional advection}
 explicitly appearance of the ``transverse flux difference'' in the unsplit
 interface state.  We rewrite this as:
 \begin{equation}
-a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} 
+a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf}
    - \frac{\Delta t}{2} v \left .\frac{\partial a}{\partial y} \right |_{i,j}
 \end{equation}
 Here, the $\hat{a}_{i+\myhalf,j,L}^{n+\myhalf}$ term is just the normal
@@ -768,7 +768,7 @@ \subsection{Unsplit multi-dimensional advection}
 compute these one-dimensional $\hat{a}$'s at the left and right every
 interface in both coordinate directions.  Note that these states are
 still one-dimensional, since we have not used any information from the
-transverse direction in their computation.  
+transverse direction in their computation.
 
 \item Solve a Riemann problem at each of these interfaces:
 \begin{eqnarray}
@@ -778,13 +778,13 @@ \subsection{Unsplit multi-dimensional advection}
                               \hat{a}_{i,j+\myhalf,R}^{n+\myhalf}) \\
 \end{eqnarray}
 These states are given the `$^T$' superscript since they are the states
-that are used in computing the transverse flux difference.  
+that are used in computing the transverse flux difference.
 
-\item Correct the 
+\item Correct the
 previously computed normal interface states (the $\hat{a}$'s) with
 the transverse flux difference:
 \begin{equation}
-a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf} 
+a_{i+\myhalf,j,L}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,L}^{n+\myhalf}
    - \frac{\Delta t}{2} v \frac{a^T_{i,j+\myhalf} - a^T_{i,j-\myhalf}}{\Delta y}
 \end{equation}
 Notice that the
@@ -793,11 +793,11 @@ \subsection{Unsplit multi-dimensional advection}
 right state, it would be those that are transverse, but to the right
 of the interface:
 \begin{equation}
-a_{i+\myhalf,j,R}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,R}^{n+\myhalf} 
+a_{i+\myhalf,j,R}^{n+\myhalf} = \hat{a}_{i+\myhalf,j,R}^{n+\myhalf}
    - \frac{\Delta t}{2} v \frac{a^T_{i+1,j+\myhalf} - a^T_{i+1,j-\myhalf}}{\Delta y}
 \end{equation}
 \end{itemize}
-A similar procedure happens at the $y$-interfaces.  
+A similar procedure happens at the $y$-interfaces.
 
 Figure~\ref{fig:unsplitstates} illustrates the steps involved in
 the construction of the $a_{i+\myhalf,j,L}^{n+\myhalf}$ state.  \MarginPar{more panels illustrating the sequence?}
@@ -821,7 +821,7 @@ \subsection{Unsplit multi-dimensional advection}
 
 Once all of the full states (normal prediction $+$ transverse flux
 difference) are computed to the left and right of all the interfaces
-($x$ and $y$), we solve another Riemann problem to find the final 
+($x$ and $y$), we solve another Riemann problem to find the final
 state on each interface.
 \begin{equation}
 a_{i+\myhalf,j}^{n+\myhalf} = \mathcal{R}(a_{i+\myhalf,j,L}^{n+\myhalf},
@@ -829,7 +829,7 @@ \subsection{Unsplit multi-dimensional advection}
 \end{equation}
 The final conservative update is then done via Eq.~\ref{eq:update2du}.
 
-See \cite{colella:1990} for more details on this unsplit method, 
+See \cite{colella:1990} for more details on this unsplit method,
 and \cite{saltzman:1994} for details of the extension to 3-d.
 
 Figures~\ref{fig:adv:gaussian-2d} and \ref{fig:adv:tophat-2d} show the
@@ -916,12 +916,12 @@ \subsection{Method-of-lines in multi-dimensions}
 \item x-face:
 \begin{equation}
 \langle (ua)_{i+\myhalf,j} \rangle = \frac{1}{\Delta y}
-    \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} (ua)_{i+\myhalf,j}\, dy 
+    \int_{y_{j-\myhalf}}^{y_{j+\myhalf}} (ua)_{i+\myhalf,j}\, dy
 \end{equation}
 \item y-face
 \begin{equation}
 \langle (va)_{i,j+\myhalf} \rangle = \frac{1}{\Delta x}
-    \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx 
+    \int_{x_{i-\myhalf}}^{x_{i+\myhalf}} (va)_{i,j+\myhalf}\, dx
 \end{equation}
 \end{itemize}
 These are similar to the expressions above, except there is no
@@ -979,7 +979,7 @@ \subsection{Method-of-lines in multi-dimensions}
 \begin{equation}
 a^n_{i,j} = A^n e^{Ii\theta}e^{Ij\phi}
 \end{equation}
-Next, to simplify things, we'll assume that $u = v$, and $\Delta x = \Delta y$, 
+Next, to simplify things, we'll assume that $u = v$, and $\Delta x = \Delta y$,
 then $\cfl = u \Delta t/ \Delta x = v \Delta t / \Delta y$.
 Inserting this into our difference equation gives:
 \begin{equation}
@@ -991,7 +991,7 @@ \subsection{Method-of-lines in multi-dimensions}
 [0, 2\pi]$ and $\phi \in [0, 2\pi]$ as a function of $\cfl$.  In doing
 this, we find that this discretization is unstable ($| {A^{n+1}}/{A^n}
 |^2 > 1$) for $\cfl > 1/2$\footnote{See the script
-  \href{https://github.com/zingale/hydro_examples/blob/master/compressible/cfl.py}{cfl.py}}.
+  \href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/cfl.py}{cfl.py}}.
 
 This restriction is $1/d$, where $d$ is the dimensionality, and it is
 analogous to the difference between the timestep limits for stability
@@ -1004,7 +1004,7 @@ \subsection{Method-of-lines in multi-dimensions}
 method in \cite{mccorquodalecolella} allows for $\Delta t \sum_{i=1}^d (|\Ub \cdot \mathbf{e}_d|)/\Delta x_d  \lesssim 1.4$.  Total variation diminishing Runge-Kutta methods are popular
 for this application~\cite{gottliebshu:1996}.
 
-% run as: 
+% run as:
 %   ./pyro.py advection_rk tophat inputs.tophat
 %   ./plot.py -W 7.7 -H 6 -o tophat_rk4_cfl08.pdf advection_rk tophat_0040
 %
@@ -1019,7 +1019,7 @@ \subsection{Method-of-lines in multi-dimensions}
   the method-of-lines approach with 4th-order Runge-Kutta time integration.
   We use $u = v = 1$ for one period on a $32\times 32$ grid.  On the left
   is $\cfl = 0.8$---this is clearly unstable.  On the right is $\cfl = 0.4$,
-  which is stable.  
+  which is stable.
   This was run as {\tt ./pyro.py advection\_rk tophat inputs.tophat}}
 \end{figure}
 
@@ -1301,7 +1301,7 @@ \subsubsection{Implementation issues with WENO schemes}
 \includegraphics[width=0.8\linewidth]{weno-converge-gaussian-rk4}
 \caption[High order WENO convergence rates for linear advection]
 {\label{fig:weno-converge-gaussian-rk4} WENO solutions for advecting a Gaussian five periods, using two different orders. A fourth order Runge-Kutta method is used for time evolution. For $r=3$ we see the expected fifth ($2 r - 1$) order convergence. For $r=5$ we see fourth order, rather than the expected ninth order, convergence. This is as the error from the time integrator dominates. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/weno.py}{weno.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/weno.py}{weno.py}}}
 \end{figure}
 %
 
@@ -1311,7 +1311,7 @@ \subsubsection{Implementation issues with WENO schemes}
 \includegraphics[width=0.8\linewidth]{weno-converge-gaussian}
 \caption[Very high order WENO convergence rates for linear advection]
 {\label{fig:weno-converge-gaussian} WENO solutions for advecting a Gaussian one periods, using four different orders. An eighth order Dormand-Price Runge-Kutta method is used for time evolution. This minimizes the time integrator error and we see convergence at order $2 r - 2$ for all schemes, although the time integrator error eventually shows in the highest order case. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/weno.py}{weno.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/weno.py}{weno.py}}}
 \end{figure}
 %
 
@@ -1341,18 +1341,18 @@ \section{Going further}
   \end{equation}
   The solution procedure is largely the same as described above.  We write:
 \begin{align}
-\phi_{i+\myhalf,j,L}^{n+\myhalf} &= \phi_{i,j}^n + 
+\phi_{i+\myhalf,j,L}^{n+\myhalf} &= \phi_{i,j}^n +
    \frac{\Delta x}{2} \frac{\partial \phi}{\partial x} +
    \frac{\Delta t}{2} \frac{\partial \phi}{\partial t} + \ldots \nonumber \\
 %
- &= \phi_{i,j}^n + 
+ &= \phi_{i,j}^n +
     \frac{\Delta x}{2} \frac{\partial \phi}{\partial x} +
     \frac{\Delta t}{2} \left [ -(\phi u)_x -(\phi v)_y \right ]_{i,j} \nonumber\\
 %
- &= \underbrace{\phi_{i,j}^n + 
+ &= \underbrace{\phi_{i,j}^n +
    \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u_{i,j} \right )
         \frac{\partial \phi}{\partial x}}_{\hat{\phi}_{i+\myhalf,j,L}^{n+\myhalf}}
-   - \frac{\Delta t}{2} \left [\phi u_x \right ]_{i,j} 
+   - \frac{\Delta t}{2} \left [\phi u_x \right ]_{i,j}
    - \frac{\Delta t}{2} \left [ (\phi v)_y \right ]_{i,j}
 \end{align}
   and upwinding is used to resolve the Riemann problem for both the
@@ -1361,7 +1361,7 @@ \section{Going further}
   according to the velocity field in much the fashion shown here, and
   is described in \S~\ref{sec:lm:density}.
 
-  
+
   For compressible hydrodynamics, we often have density-weighted quantities
   that we advect.  This extension is described in \S~\ref{sec:euler:further}.
 
@@ -1370,7 +1370,7 @@ \section{Going further}
 
 \section{\pyro\ experimentation}
 
-To gain some experiences with these ideas, we can use the advection 
+To gain some experiences with these ideas, we can use the advection
 solver in \pyro\ (see Appendix~\ref{app:pyro} to get started).
 The \pyro\ advection solver implements the second-order unsplit
 advection algorithm described in the previous sections.  To run
@@ -1402,4 +1402,3 @@ \section{\pyro\ experimentation}
 and compare the results to the unsplit version.  Pick a suitable norm
 and measure the convergence of the methods.}
 \end{exercise}
-
diff --git a/burgers/burgers.tex b/burgers/burgers.tex
index d49d0de..d924031 100644
--- a/burgers/burgers.tex
+++ b/burgers/burgers.tex
@@ -71,7 +71,7 @@ \section{Burgers' equation}
           velocity distribution with low velocity left of high
           velocity (bottom) Approximate characteristic structure for
           Burgers' equation, using $u_0 = u(t)$.  Note that after a
-          short period of time, the characteristics diverge---this is 
+          short period of time, the characteristics diverge---this is
           a rarefaction.}
 \end{figure}
 
@@ -203,7 +203,7 @@ \section{Burgers' equation}
   This initial velocity state creates a divergent flow.  The curves
   are shown 0.02~s apart, with the darker grayscale representing later
   in
-  time. \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/burgers/burgers.py}{burgers.py}}}
+  time. \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/burgers/burgers.py}{burgers.py}}}
 \end{figure}
 
 
@@ -224,7 +224,7 @@ \section{Burgers' equation}
   The initial conditions here are sinusoidal and the solution quickly
   steepens into a shock.The curves are shown 0.02~s apart, with the
   darker grayscale representing later in
-  time. \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/burgers/burgers.py}{burgers.py}}}
+  time. \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/burgers/burgers.py}{burgers.py}}}
 \end{figure}
 
 Figure~\ref{fig:burgers-shock} shows the solution to Burgers'
@@ -276,12 +276,12 @@ \section{Burgers' equation}
 conservative and non-conservation formulations of Burgers' equation
 as:
 \begin{equation}
-\frac{u^{n+1}_i - u^n_i}{\Delta t} = 
+\frac{u^{n+1}_i - u^n_i}{\Delta t} =
    -\frac{1}{2} \frac{(u^n_i)^2 - (u^n_{i-1})^2}{\Delta x}
 \end{equation}
 and
 \begin{equation}
-  \frac{u_{i}^{n+1} - u_i^n}{\Delta t} = - 
+  \frac{u_{i}^{n+1} - u_i^n}{\Delta t} = -
       \frac{u_i^n (u_i^n - u_{i-1}^n)}{\Delta x}
 \end{equation}
 (Note: these discretizations are upwind so long as $u > 0$).
diff --git a/diffusion/diffusion.tex b/diffusion/diffusion.tex
index acb8909..2efe0cc 100644
--- a/diffusion/diffusion.tex
+++ b/diffusion/diffusion.tex
@@ -25,11 +25,11 @@ \section{Diffusion}
 will see, we can cast our differenced equations in a form reminiscent
 of an elliptic problem.
 
-In one-dimension, and assuming that the diffusion coefficient, $k$, is 
+In one-dimension, and assuming that the diffusion coefficient, $k$, is
 constant, we have
 \begin{equation}
-\frac{\partial \phi}{\partial t} = 
-  \frac{\partial }{\partial x} 
+\frac{\partial \phi}{\partial t} =
+  \frac{\partial }{\partial x}
   \left ( k \frac{\partial \phi}{\partial x} \right )
 \end{equation}
 
@@ -38,7 +38,7 @@ \section{Diffusion}
 part of the energy equation in compressible flow), species/mass
 diffusion for multi-species flows, or the viscous terms in
 incompressible flows.  In this form, the diffusion coefficient (or
-conductivity), $k$, can be a function of $x$, or even $\phi$.  
+conductivity), $k$, can be a function of $x$, or even $\phi$.
 
 
 We will consider a constant diffusion coefficient as our model problem:
@@ -52,12 +52,12 @@ \section{Diffusion}
 \section{Explicit differencing}
 
 A nice overview of the discretizations for diffusion and their stability
-is given in \cite{richtmyermorton}.  
+is given in \cite{richtmyermorton}.
 The simplest way to difference this equation is {\em explicit} in time
 (i.e., the righthand side depends only on the old state):
 \begin{equation}
 \label{eq:diff:explicitdiff}
-\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
   k \frac{\phi_{i+1}^n - 2\phi_i^n + \phi_{i-1}^n}{\Delta x^2}
 \end{equation}
 This is second-order accurate in space, but only first order accurate in
@@ -82,7 +82,7 @@ \section{Explicit differencing}
 \begin{equation}
 \Delta t = \frac{C}{2} \frac{\Delta x^2}{k}
 \end{equation}
-where $C$ is a constant, $C < 1$. 
+where $C$ is a constant, $C < 1$.
 Note the $\Delta x^2$ dependence---this constraint can become really
 restrictive.
 
@@ -142,7 +142,7 @@ \section{Explicit differencing}
 Diffusion of a Gaussian using the explicit differencing of
 Eq.~\ref{eq:diff:explicitdiff} with 64 zones and C = 0.8, shown at
 several times.  The dotted line is the analytic
-solution.  \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
+solution.  \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
 \end{figure}
 
 
@@ -153,7 +153,7 @@ \section{Explicit differencing}
 Diffusion of a Gaussian using the explicit differencing of
 Eq.~\ref{eq:diff:explicitdiff} with different resolutions.
 several times.  The dotted line is the analytic
-solution.  \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
+solution.  \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
 \end{figure}
 
 As with advection, if we exceed the timestep limit ($C > 1$), then
@@ -166,7 +166,7 @@ \section{Explicit differencing}
 Diffusion of a Gaussian using the explicit differencing of
 Eq.~\ref{eq:diff:explicitdiff} with 64 zones, but a timestep
 with $C > 1$, showing that the solution is unstable.
- \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
+ \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
 \end{figure}
 
 Our spatial order-of-accuracy is second-order.
@@ -190,9 +190,9 @@ \section{Explicit differencing}
 \centering
 \includegraphics[width=\linewidth]{diffexplicit-converge-0_8}
 \caption[Error convergence of explicit diffusion]{\label{fig:diffexplicit_converge}
-Error convergence with resolution of the explicit diffusion using the differencing of 
+Error convergence with resolution of the explicit diffusion using the differencing of
 Eq.~\ref{eq:diff:explicitdiff}
- \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
+ \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_explicit.py}{diffusion\_explicit.py}}}
 \end{figure}
 
 \section{Implicit with direct solve}
@@ -202,7 +202,7 @@ \section{Implicit with direct solve}
 diffusion.  A backward-Euler implicit discretization would be:
 \begin{equation}
 \label{eq:diff:foimplicit}
-\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
   k \frac{\phi_{i+1}^{n+1} - 2\phi_i^{n+1} + \phi_{i-1}^{n+1}}{\Delta x^2}
 \end{equation}
 The only difference with Eq.~\ref{eq:diff:explicitdiff} is the
@@ -226,7 +226,7 @@ \section{Implicit with direct solve}
 \end{equation}
 This is a set of coupled algebraic equations.  We can write this in
 matrix form.  Using a cell-centered grid, we will solve for the values
-$[\mathrm{lo},\mathrm{hi}]$.  
+$[\mathrm{lo},\mathrm{hi}]$.
 The implicit method can use any $C > 0$.
 
 We specify boundary conditions by modifying the stencil
@@ -239,7 +239,7 @@ \section{Implicit with direct solve}
 \end{equation}
 and substituting this into Eq~\ref{eq:diff:implicit}, the update for the leftmost cell is:
 \begin{equation}
- (1 + \alpha) \phi_\mathrm{lo}^{n+1} -\alpha \phi_\mathrm{lo+1}^{n+1}  = 
+ (1 + \alpha) \phi_\mathrm{lo}^{n+1} -\alpha \phi_\mathrm{lo+1}^{n+1}  =
   \phi_\mathrm{lo}^n
 \end{equation}
 If we choose Dirichlet BCs on the right ($\phi |_{x=x_l} = A$), then:
@@ -299,7 +299,7 @@ \subsubsection{Crank-Nicolson time discretization}
 A second-order in time discretization requires us to center the
 righthand side in time.  We do this as:
 \begin{equation}
-\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+\frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
   \frac{k}{2} \left ( \frac{\phi_{i+1}^{n} - 2\phi_i^{n} + \phi_{i-1}^{n}}{\Delta x^2} +
                       \frac{\phi_{i+1}^{n+1} - 2\phi_i^{n+1} + \phi_{i-1}^{n+1}}{\Delta x^2} \right )
 \end{equation}
@@ -367,7 +367,7 @@ \subsubsection{Crank-Nicolson time discretization}
 \right )
 \end{equation}
 
-Figure~\ref{fig:diffuse} shows the result of using $\alpha = 0.8$ and $\alpha = 8.0$.  We 
+Figure~\ref{fig:diffuse} shows the result of using $\alpha = 0.8$ and $\alpha = 8.0$.  We
 see that they are both stable, but that the smaller timestep is closer to the analytic
 solution (especially at early times).
 
@@ -379,7 +379,7 @@ \subsubsection{Crank-Nicolson time discretization}
   Implicit diffusion of a Gaussian (with Crank-Nicolson
   discretization) with $C = 0.8$ and $C = 8.0$.  The exact solution at
   each time is shown as the dotted
-  line. \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_implicit.py}{diffusion\_implicit.py}}}
+  line. \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_implicit.py}{diffusion\_implicit.py}}}
 \end{figure}
 
 \begin{exercise}[1-d implicit diffusion]
@@ -389,7 +389,7 @@ \subsubsection{Crank-Nicolson time discretization}
   that above.  Use a timestep close to the explicit step, a grid with
   N = 128 zones.
 
-  Use a Gaussian for your initial conditions (as you did for the 
+  Use a Gaussian for your initial conditions (as you did for the
   explicit problem.}
 \end{exercise}
 
@@ -397,17 +397,17 @@ \subsubsection{Crank-Nicolson time discretization}
 \section{Implicit multi-dimensional diffusion via multigrid}
 \label{diff:sec:implicit_mg}
 
-Instead of doing a direct solve of the matrix form of the system, we 
-can use multigrid techniques.  Consider the Crank-Nicolson system we just 
+Instead of doing a direct solve of the matrix form of the system, we
+can use multigrid techniques.  Consider the Crank-Nicolson system we just
 looked at:
 \begin{equation}
-\frac{\phi^{n+1}_i - \phi^n_i}{\Delta t} = 
+\frac{\phi^{n+1}_i - \phi^n_i}{\Delta t} =
    \frac{1}{2} \left ( k \nabla^2 \phi^n_i + k \nabla^2 \phi^{n+1}_i \right )
 \end{equation}
-Grouping all the 
+Grouping all the
 $n+1$ terms on the left, we find:
 \begin{equation}
-\phi^{n+1}_i - \frac{\Delta t}{2} k \nabla^2 \phi^{n+1}_i = 
+\phi^{n+1}_i - \frac{\Delta t}{2} k \nabla^2 \phi^{n+1}_i =
     \phi^n_i + \frac{\Delta t}{2} k \nabla^2 \phi^n_i
 \end{equation}
 This is in the form of a constant-coefficient Helmholtz equation,
@@ -433,7 +433,7 @@ \section{Implicit multi-dimensional diffusion via multigrid}
 \label{eq:diff:discrete1d}
 \phi_{i} \leftarrow
  \left .    \left ( f_{i} + \frac{\beta}{\Delta x^2} [\phi_{i+1}
-                             + \phi_{i-1}] \right ) \middle / 
+                             + \phi_{i-1}] \right ) \middle /
 \left ( \alpha + \frac{2 \beta}{\Delta x^2}  \right )  \right .
 \end{equation}
 Note that when you take $\alpha = 0$ and $\beta = -1$ in
@@ -466,7 +466,7 @@ \section{Implicit multi-dimensional diffusion via multigrid}
 \caption[Comparison of 2-d implicit diffusion with analytic
   solution]{\label{fig:diffusion:multid} A comparison of 2-d implicit
   diffusion from Figure~\ref{fig:diff:twodcompare} with the analytic solution at several times.
-  The 2-d solution was averaged over angles to yield a profile as a 
+  The 2-d solution was averaged over angles to yield a profile as a
   function of radius, using the {\tt
     gauss\_diffusion\_compare.py} in \pyro.}
 \end{figure}
@@ -480,10 +480,10 @@ \section{Implicit multi-dimensional diffusion via multigrid}
 
 
 \begin{exercise}[Implicit multi-dimensional diffusion]
-{The diffusion solver in \pyro\ uses Crank-Nicolson differencing in time.  
+{The diffusion solver in \pyro\ uses Crank-Nicolson differencing in time.
 Modify the solver to do first-order backward Euler.  This will change
 the form of the linear system (coefficients and righthand side), but
-should not require any changes to the multigrid solver itself. 
+should not require any changes to the multigrid solver itself.
 
 Compare the solution with the Crank-Nicolson solution for a very coarsely-resolved
 Gaussian and a finely-resolved Gaussian.
@@ -509,21 +509,21 @@ \subsection{Convergence}
 \caption[Under-resolved Crank-Nicolson diffusion]
 {\label{fig:diff:cnunstable} Crank-Nicolson diffusion of a Gaussian
 on under-resolved initial data with a large timestep.  Here we use
-$64$ zones and $C= 10$.  \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diffusion_implicit.py}{diffusion\_implicit.py}}}
+$64$ zones and $C= 10$.  \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diffusion_implicit.py}{diffusion\_implicit.py}}}
 \end{figure}
 
 
 \begin{figure}
 \centering
 \includegraphics[width=0.75\linewidth]{diffimplicit-converge-0_8} \\
-\includegraphics[width=0.75\linewidth]{diffimplicit-converge-8_0} 
+\includegraphics[width=0.75\linewidth]{diffimplicit-converge-8_0}
 \caption[Convergence of diffusion methods]
 {\label{fig:diff:convergence} Convergence of the explicit, backward-difference,
 and Crank-Nicolson diffusion methods for $C = 0.8$ (top) and $C = 8.0$ (bottom).
 For the latter case, the explicit method is not valid and is not shown.
 We see that the Crank-Nicolson method has the lowest error, and when resolving
 the data, has second-order convergence.\\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/diffusion/diff_converge.py}{diff\_converge.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/diffusion/diff_converge.py}{diff\_converge.py}}}
 \end{figure}
 
 
@@ -551,7 +551,7 @@ \section{Non-constant Conductivity}
 For the case where $k = k(x)$,
   we discretize as:
   \begin{equation}
-  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
         \frac{ \{ k \nabla \phi \}_{i+\myhalf} -
                \{ k \nabla \phi \}_{i-\myhalf}}{\Delta x}
   \end{equation}
@@ -566,15 +566,15 @@ \section{Non-constant Conductivity}
  \frac{1}{k_{i+\myhalf}} = \frac{1}{2} \left (\frac{1}{k_i} + \frac{1}{k_{i+1}} \right )
  \end{equation}
 The actual form should be motivated by the physics
- 
+
 Slightly more complicated are state-dependent transport
 coefficients---the transport coefficients themselves depend on the
 quantity being diffused:
   \begin{equation}
-  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
         \frac{1}{2} \left \{
                \nabla \cdot [ k(\phi^n) \nabla \phi^n ]_i +
-               \nabla \cdot [ k(\phi^{n+1}) \nabla \phi^{n+1} ]_i 
+               \nabla \cdot [ k(\phi^{n+1}) \nabla \phi^{n+1} ]_i
                \right \}
   \end{equation}
   (for example, with thermal diffusion, the conductivity can
@@ -583,25 +583,25 @@ \section{Non-constant Conductivity}
   the transport coefficients evaluated at the old time, giving a provisional
   state, $\phi^\star$:
   \begin{equation}
-  \frac{\phi_i^\star - \phi_i^n}{\Delta t} = 
+  \frac{\phi_i^\star - \phi_i^n}{\Delta t} =
         \frac{1}{2} \left \{
                \nabla \cdot [ k(\phi^n) \nabla \phi^n ] +
-               \nabla \cdot [ k(\phi^n) \nabla \phi^\star ] 
+               \nabla \cdot [ k(\phi^n) \nabla \phi^\star ]
                \right \}
   \end{equation}
   Then we redo the diffusion, evaluating $k$ with $\phi^\star$ to
   center the righthand side in time, giving the new state, $\phi^{n+1}$:
   \begin{equation}
-  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} = 
+  \frac{\phi_i^{n+1} - \phi_i^n}{\Delta t} =
         \frac{1}{2} \left \{
                \nabla \cdot [ k(\phi^n) \nabla \phi^n ] +
-               \nabla \cdot [ k(\phi^\star) \nabla \phi^{n+1} ] 
+               \nabla \cdot [ k(\phi^\star) \nabla \phi^{n+1} ]
                \right \}
   \end{equation}
   This is the approach used, for example, in \cite{SNpaper}.
 
 
-\section{Diffusion in Hydrodynamics} 
+\section{Diffusion in Hydrodynamics}
 
 Often we find diffusion represented as one of many physical processes
 in a single equation.  For example, consider the internal energy
@@ -633,6 +633,5 @@ \section{Diffusion in Hydrodynamics}
  techniques described above.  Note: if the transport coefficients
  (including $c_v$, $e_\rho$) are dependent on $e$, then we still need
  to do a predictor-corrector method here.  This is discussed, for
- example, in \cite{SNpaper,malone:2011}.  A simple case for this type of 
+ example, in \cite{SNpaper,malone:2011}.  A simple case for this type of
  advection-diffusion is also shown in \S~\ref{sec:multiphys:diffreact}.
-
diff --git a/finite-volume/finite-volume.tex b/finite-volume/finite-volume.tex
index 77702b3..2a2c66a 100644
--- a/finite-volume/finite-volume.tex
+++ b/finite-volume/finite-volume.tex
@@ -118,7 +118,7 @@ \section{Grid basics}
 space, keep track of the total amount of $f$ in the zone (indicated as the shaded regions).  Since we generally don't know how $f$ varies in the zone,
 we will typically talk about the average of $f$, $\favg{f}_i$, over
 the zone, and represent this by a horizontal.  The total amount of $f$ in
-the zone is then simply $\Delta x \favg{f}_i$. 
+the zone is then simply $\Delta x \favg{f}_i$.
 We label the left and right edges of a zone with half-integer indices
 $i-\myhalf$ and $i+\myhalf$.  The physical coordinate of the center of the zone
 is the same as in the cell-centered finite-difference case.
@@ -235,7 +235,7 @@ \section{Finite-volume grids}
 \end{exercise}
 
 The {\sf Jupyter} notebook
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/finite-volume/conservative-interpolation.ipynb}{conservative-interpolation.ipynb}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/finite-volume/conservative-interpolation.ipynb}{conservative-interpolation.ipynb}}
 shows how to derive these interpolants using {\sf SymPy}, and gives
 higher-order interpolants.
 
@@ -372,7 +372,7 @@ \section{Numerical implementation details}
   ihi = nx + ng - 1
   state = np.zeros((nx + 2*ng, nvar))
 \end{lstlisting}
-in python (using NumPy), where {\tt nx} is the number of zones in the 
+in python (using NumPy), where {\tt nx} is the number of zones in the
 domain and {\tt ng} is the number of ghost cells.  For some problems,
 there might be more than one state variable, so {\tt nvar}, is the number
 of state variables that live in each zone.  We would then loop from {\tt ilo}
diff --git a/higher-order/higher-order-burgers.tex b/higher-order/higher-order-burgers.tex
index 4ee44ed..d78ad6a 100644
--- a/higher-order/higher-order-burgers.tex
+++ b/higher-order/higher-order-burgers.tex
@@ -58,7 +58,7 @@ \section{WENO methods, nonlinear equations, and flux-splitting}
 \includegraphics[width=0.8\linewidth]{weno-vs-plm-burger}
 \caption[Comparing PLM and WENO methods for Burgers' equation]
 {\label{fig:weno-vs-plm-burger} Numerical solutions to Burgers' equation using the piecewise linear method as in figure~\ref{fig:burgers-shock} and two flux-split WENO methods as outlined in section~\ref{sec:weno-burgers-fvs}. This uses the same sinusoidal data as figure~\ref{fig:burgers-shock} shortly after shock formation. Note that the WENO methods are not obviously any better at resolving the shock than the simpler piecewise linear approach.\\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/weno_burgers.py}{weno\_burgers.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/weno_burgers.py}{weno\_burgers.py}}}
 \end{figure}
 %
 
@@ -87,7 +87,7 @@ \section{WENO methods, nonlinear equations, and flux-splitting}
 For $r=3$ we see the expected fifth ($2 r - 1$) order convergence.
 For $r=5$ we see seventh order convergence briefly, before the error from the
 time integrator dominates. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/advection/weno_burgers.py}{weno\_burgers.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/advection/weno_burgers.py}{weno\_burgers.py}}}
 \end{figure}
 %
 
diff --git a/higher-order/higher-order-euler.tex b/higher-order/higher-order-euler.tex
index 54653ba..ece4e6f 100644
--- a/higher-order/higher-order-euler.tex
+++ b/higher-order/higher-order-euler.tex
@@ -46,7 +46,7 @@ \section{WENO methods for the Euler equations}
 \includegraphics[width=0.8\linewidth]{weno-euler-r3}
 \caption[WENO $r=3$ for the Sod test]
 {\label{fig:weno-euler-r3} The Sod test solved with WENO methods, $r=3$, using characteristic-wise and component-wise reconstruction, using 64 zones and a cfl of $0.5$. To compare with Godunov's method see figure~\ref{fig:Euler:sod:god} and with PPM see figure~\ref{fig:Euler:sod:ppm}. The WENO method does not capture discontinuities, especially the contact, as well as the PPM method which is specifically designed for the Euler equations. The small oscillations visible in the component-wise reconstruction are damped by using the more complex and expensive characteristic-wise approach. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
 \end{figure}
 %
 
@@ -56,7 +56,7 @@ \section{WENO methods for the Euler equations}
 \includegraphics[width=0.8\linewidth]{weno-euler-r5}
 \caption[WENO $r=5$ for the Sod test]
 {\label{fig:weno-euler-r5} The Sod test solved with WENO methods, $r=5$, using characteristic-wise and component-wise reconstruction, using 64 zones and a cfl of $0.5$. Compare with figure~\ref{fig:weno-euler-r3} to see the effect of the higher order of WENO scheme. The oscillations in the component-wise approach are much more pronounced as the order of the reconstruction is increased, and the characteristic-wise approach continues to help with this. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
 \end{figure}
 %
 
@@ -77,7 +77,7 @@ \section{WENO methods for the Euler equations}
 \includegraphics[width=0.8\linewidth]{weno-euler-rarefaction-r3}
 \caption[WENO $r=3$ for the double rarefaction test]
 {\label{fig:weno-euler-rarefaction-r3} The double rarefaction test solved with WENO methods, $r=3$, using characteristic-wise and component-wise reconstruction, using 64 zones and a cfl of $0.5$. To compare with Godunov's method see figure~\ref{fig:Euler:doublerare:god} and with PPM see figure~\ref{fig:Euler:doublerare:ppm}. The WENO method captures the edges of the rarefactions almost as well as the PPM method, but is as almost as poor as Godunov's method at the trivial contact at the center. A less diffusive flux-splitting may help here. Increasing the reconstruction order has little effect.\\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/compressible/weno_euler.py}{weno\_euler.py}}}
 \end{figure}
 %
 
diff --git a/hydro_examples/hydro_examples.tex b/hydro_examples/hydro_examples.tex
index 2a06be9..7577254 100644
--- a/hydro_examples/hydro_examples.tex
+++ b/hydro_examples/hydro_examples.tex
@@ -10,13 +10,13 @@ \section{Introduction}
 \section{Getting \hydroex}
 
 The \hydroex\ codes are hosted on github:\\
-\url{https://github.com/zingale/hydro_examples}
+\url{https://github.com/python-hydro/hydro_examples}
 
 There are a few ways to get them.  The simplest way is to just clone
 from github on the commandline:
 
 \begin{verbatim}
-git clone https://github.com/zingale/hydro_examples
+git clone https://github.com/python-hydro/hydro_examples
 \end{verbatim}
 
 This will create a local git repo on your machine with all of the
@@ -45,7 +45,7 @@ \section{\hydroex\ codes}
 \begin{itemize}
 
 \item {\tt advection/}
-  
+
   \begin{itemize}
   \item {\tt advection.py}: a 1-d second-order linear advection solver
     with a wide range of limiters.
@@ -72,7 +72,7 @@ \section{\hydroex\ codes}
 \item {\tt compressible/}
 
   \begin{itemize}
-  \item {\tt cfl.py}: a simple brute force script used to find the stability 
+  \item {\tt cfl.py}: a simple brute force script used to find the stability
     limit of a multidimensional advection discretization.
 
   \item {\tt euler.ipynb}: an {\sf Jupyter} notebook using {\sf SymPy}
@@ -83,9 +83,9 @@ \section{\hydroex\ codes}
     {\sf SymPy} that derives the eigensystem for the
     primitive-variable form of various forms of Euler equations with a
     general EOS (with either $(\rho e)$ or $\gamma_e$ augmenting the
-    system) and for gray radiation hydrodynamics.  Additionally, this 
+    system) and for gray radiation hydrodynamics.  Additionally, this
     notebook computes the full form of the interface states.
-      
+
   \item {\tt riemann.py}: the main Riemann class for gamma-law Euler
     equations.  The {\tt RiemannProblem} class has methods to find
     the star state, plot the Hugoniot curves, and sample the solution.
@@ -127,7 +127,7 @@ \section{\hydroex\ codes}
       is diffused.
 
   \end{itemize}
-  
+
 \item {\tt elliptic/}
 
   \begin{itemize}
@@ -145,7 +145,7 @@ \section{\hydroex\ codes}
 \item {\tt incompressible/}
 
   \begin{itemize}
-  \item {\tt project.py}: a simple example of using a projection 
+  \item {\tt project.py}: a simple example of using a projection
     to recover a divergence-free velocity field.
   \end{itemize}
 
@@ -158,7 +158,7 @@ \section{\hydroex\ codes}
     accuracy.
 
   \item {\tt mg\_test.py}: a simple driver for the multigrid solver.
-    This sets up and solves a Poisson problem and plots the behavior 
+    This sets up and solves a Poisson problem and plots the behavior
     of the solution as a function of V-cycle number.
 
   \item{\tt multigrid.py}: a multigrid class for cell-centered data.
@@ -180,7 +180,7 @@ \section{\hydroex\ codes}
   \begin{itemize}
   \item {\tt burgersvisc.py}: solve the viscous Burgers equation.
      The advective terms are treated explicitly with a second-order
-     accurate method.  The diffusive term is solved using an 
+     accurate method.  The diffusive term is solved using an
      implicit Crank-Nicolson discretization.  The overall coupling
      is second-order accurate.
 
@@ -191,10 +191,8 @@ \section{\hydroex\ codes}
     are evolved using the VODE ODE solver (via SciPy).  The two
     processes are coupled together using Strang-splitting to be
     second-order accurate in time.
-     
+
 
   \end{itemize}
 
 \end{itemize}
-
-
diff --git a/incompressible/incompressible.tex b/incompressible/incompressible.tex
index e496158..9a09186 100644
--- a/incompressible/incompressible.tex
+++ b/incompressible/incompressible.tex
@@ -20,9 +20,9 @@ \section{Incompressible flow}
 \begin{equation}
 \nabla \cdot \Ub = 0
 \end{equation}
-(since $D\rho / Dt = 0$).  The incompressible fluid approximation is 
+(since $D\rho / Dt = 0$).  The incompressible fluid approximation is
 a reasonable approximation when the Mach number of the fluid is small
-($\ll 1$).  To complete the system, we add the momentum equation.  If 
+($\ll 1$).  To complete the system, we add the momentum equation.  If
 we take the density to be constant everywhere in the domain (not just
 in our fluid element), then we have:
 \begin{equation}
@@ -32,7 +32,7 @@ \section{Incompressible flow}
 any equation of state.  This system is closed as written.  The value
 of $p$ is determined such that the velocity constraint is satisfied.
 
-We can gain more insight into the applicability of the incompressible 
+We can gain more insight into the applicability of the incompressible
 equations by doing an asymptotic expansion.  Starting with the momentum
 equation from the Euler equations:
 \begin{equation}
@@ -80,9 +80,9 @@ \section{Incompressible flow}
 \section{Projection methods}
 
 The basic idea behind a projection method is that any vector field
-can be decomposed into a divergence free part and the gradient of 
-a scalar (this is sometimes called a {\em Hodge decomposition}).  
-Given a velocity field $\Ub^\star$, we can express it in terms of the 
+can be decomposed into a divergence free part and the gradient of
+a scalar (this is sometimes called a {\em Hodge decomposition}).
+Given a velocity field $\Ub^\star$, we can express it in terms of the
 divergence free part $\Ub^d$ and a scalar, $\phi$ as:
 \begin{equation}
 \Ub^\star = \Ub^d + \nabla \phi
@@ -150,7 +150,7 @@ \section{Projection methods}
 by adding the gradient of a scalar, $\phi$, of the form:
 \begin{equation}
 \label{eq:incomp:projphi}
-\phi = \frac{1}{10} \cos(2\pi y) \cos(2\pi x) 
+\phi = \frac{1}{10} \cos(2\pi y) \cos(2\pi x)
 \end{equation}
 yielding a new velocity, $\Ub_p = \Ub + \nabla \phi$.
 The result is the middle panel in the figure, and is not divergence free.
@@ -163,7 +163,7 @@ \section{Projection methods}
 \begin{equation}
 \nabla \cdot \Ub = \frac{u_{i+1,j} - u_{i-1,j}}{2\dx} + \frac{v_{i,j+1} - v_{i,j-1}}{2\dy}
 \end{equation}
-Together, these represent an approximation projection. 
+Together, these represent an approximation projection.
 
 \begin{figure}[t]
 \centering
@@ -175,7 +175,7 @@ \section{Projection methods}
   scalar, Eq.~\ref{eq:incomp:projphi} (middle) and then a projection
   was done to recover the original velocity field (bottom).  In this
   case, the L$_2$ norm of the error on the velocity (comparing before
-  and after projection) is $0.000189$. \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/incompressible/project.py}{project.py}}}
+  and after projection) is $0.000189$. \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/incompressible/project.py}{project.py}}}
 \end{figure}
 
 \begin{exercise}[An approximate projection]
@@ -183,7 +183,7 @@ \section{Projection methods}
 the elliptic equation.  Start with the velocity field described
 by Eqs.~\ref{eq:incomp:proju} and \ref{eq:incomp:projv} and the scalar
 from Eq.~\ref{eq:incomp:projphi}.  Compute a poluted velocity field,
-$\Ub_p = \Ub + \nabla \phi$ and then project it to recover the 
+$\Ub_p = \Ub + \nabla \phi$ and then project it to recover the
 original velocity field.
 \end{exercise}
 
@@ -192,7 +192,7 @@ \section{Cell-centered approximate projection solver}
 Here we describe an incompressible algorithm that uses cell-centered
 data throughout---$\Ub$ and $p$ are both cell-centered.  The
 projection at the end is an approximate projection.  The basic
-algorithm flow is 
+algorithm flow is
 \begin{itemize}
 \item Create the time-centered advective velocities through the faces
   of the zones.
@@ -241,12 +241,12 @@ \subsection{Advective velocity}
 We follow ABS.  Our velocity evolution system (writing out the
 individual components of $\Ub$: $u$ and $v$) is
 \begin{eqnarray}
-\frac{\partial u}{\partial t} &=& -u \frac{\partial u}{\partial x} 
-                                  -v \frac{\partial u}{\partial y} 
+\frac{\partial u}{\partial t} &=& -u \frac{\partial u}{\partial x}
+                                  -v \frac{\partial u}{\partial y}
                                   -\frac{\partial p}{\partial x} = 0 \\
-\frac{\partial v}{\partial t} &=& -u \frac{\partial v}{\partial x} 
-                                  -v \frac{\partial v}{\partial y} 
-                                  -\frac{\partial p}{\partial y} = 0 
+\frac{\partial v}{\partial t} &=& -u \frac{\partial v}{\partial x}
+                                  -v \frac{\partial v}{\partial y}
+                                  -\frac{\partial p}{\partial y} = 0
 \end{eqnarray}
 Our goal in this step is to predict time-centered interface values of
 the normal velocity ($u$ on x-edges and $v$ on y-edges).  The
@@ -256,16 +256,16 @@ \subsection{Advective velocity}
 right states which we will resolve by solving a Riemann problem.  The
 left interface state of $u$ at $i+\myhalf,j$ is found as:
 \begin{eqnarray}
-u_{i+\myhalf,j,L}^{n+\myhalf} 
-  &=& u_{i,j} 
+u_{i+\myhalf,j,L}^{n+\myhalf}
+  &=& u_{i,j}
     + \frac{\Delta x}{2} \left . \frac{\partial u}{\partial x} \right |_{i,j}
     + \frac{\Delta t}{2} \left . \frac{\partial u}{\partial t} \right |_{i,j}\\
-  &=& u_{i,j} 
+  &=& u_{i,j}
     + \frac{\Delta x}{2} \left . \frac{\partial u}{\partial x} \right |_{i,j}
     + \frac{\Delta t}{2} \left . \left (-u \frac{\partial u}{\partial x}
                                 -v \frac{\partial u}{\partial y}
                                 -\frac{\partial p}{\partial x} \right ) \right |_{i,j}\\
-  &=& u_{i,j} 
+  &=& u_{i,j}
     + \frac{\Delta x}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u_{i,j} \right )
                          \left .  \frac{\partial u}{\partial x}\right |_{i,j}
     - \frac{\Delta t}{2} \left ( v \frac{\partial u}{\partial y}\right )_{i,j}
@@ -275,22 +275,22 @@ \subsection{Advective velocity}
 x$, where $\Dux_{i,j}$ is the limited slope of $u$ in the
 $x$-direction in zone $i,j$.  Our interface state is then:
 \begin{equation}
-u_{i+\myhalf,j,L}^{n+\myhalf} 
+u_{i+\myhalf,j,L}^{n+\myhalf}
     = \underbrace{u_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u \right ) \Dux_{i,j}}_{\equiv \hat{u}_{i+\myhalf,j,L}^{n+\myhalf}}
     - \underbrace{\frac{\Delta t}{2} \left ( v \frac{\partial u}{\partial y} \right )_{i,j}}_{\mathrm{transverse~term}}
     - \frac{\Delta t}{2} \left . \frac{\partial p}{\partial x} \right |_{i,j}
 \label{eq:uMAC_l}
-\end{equation} 
+\end{equation}
 
 
 Similarly, for $v$ through the $y$ faces, we find:
 \begin{equation}
-v_{i,j+\myhalf,L}^{n+\myhalf} 
+v_{i,j+\myhalf,L}^{n+\myhalf}
     = \underbrace{v_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta x} v \right ) \Dvy_{i,j}}_{\equiv \hat{v}_{i,j+\myhalf,L}^{n+\myhalf}}
     - \underbrace{\frac{\Delta t}{2} \left ( u \frac{\partial v}{\partial x} \right )_{i,j}}_{\mathrm{transverse~term}}
     - \frac{\Delta t}{2} \left . \frac{\partial p}{\partial y} \right |_{i,j}
 \label{eq:vMAC_l}
-\end{equation} 
+\end{equation}
 
 As indicated above (and following ABS and the similar notation used by
 Colella~\cite{colella:1990}), we denote the quantities that will be used to evaluate
@@ -306,51 +306,51 @@ \subsection{Advective velocity}
 \begin{eqnarray}
 \hat{u}_{i+\myhalf,j,L}^{n+\myhalf} &=&
    u_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u_{i,j} \right )
-       \Dux_{i,j} \\ 
+       \Dux_{i,j} \\
 %
 \hat{u}_{i+\myhalf,j,R}^{n+\myhalf} &=&
    u_{i+1,j} - \frac{1}{2} \left ( 1 + \frac{\Delta t}{\Delta x} u_{i+1,j} \right )
-       \Dux_{i+1,j} 
+       \Dux_{i+1,j}
 \end{eqnarray}
 
 \noindent $v$ on $x$-interfaces:
 \begin{eqnarray}
 \hat{v}_{i+\myhalf,j,L}^{n+\myhalf} &=&
    v_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta x} u_{i,j} \right )
-       \Dvx_{i,j} \\ 
+       \Dvx_{i,j} \\
 %
 \hat{v}_{i,j+\myhalf,R}^{n+\myhalf} &=&
    v_{i+1,j} - \frac{1}{2} \left ( 1 + \frac{\Delta t}{\Delta x} u_{i+1,j} \right )
-       \Dvx_{i+1,j} 
+       \Dvx_{i+1,j}
 \end{eqnarray}
 
 \noindent $u$ on $y$-interfaces:
 \begin{eqnarray}
 \hat{u}_{i,j+\myhalf,L}^{n+\myhalf} &=&
    u_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta y} v_{i,j} \right )
-       \Duy_{i,j} \\ 
+       \Duy_{i,j} \\
 %
 \hat{u}_{i,j+\myhalf,R}^{n+\myhalf} &=&
    u_{i,j+1} - \frac{1}{2} \left ( 1 + \frac{\Delta t}{\Delta y} v_{i,j+1} \right )
-       \Duy_{i,j+1} 
+       \Duy_{i,j+1}
 \end{eqnarray}
 
 \noindent $v$ on $y$-interfaces:
 \begin{eqnarray}
 \hat{v}_{i,j+\myhalf,L}^{n+\myhalf} &=&
    v_{i,j} + \frac{1}{2} \left ( 1 - \frac{\Delta t}{\Delta y} v_{i,j} \right )
-       \Dvy_{i,j} \\ 
+       \Dvy_{i,j} \\
 %
 \hat{v}_{i,j+\myhalf,R}^{n+\myhalf} &=&
    v_{i,j+1} - \frac{1}{2} \left ( 1 + \frac{\Delta t}{\Delta y} v_{i,j+1} \right )
-       \Dvy_{i,j+1} 
+       \Dvy_{i,j+1}
 \end{eqnarray}
 
 Note that the `right' state is constructed using the data to the right
 of the interface.  Also note that in constructing these transverse
 velocities, we do not include the $p$ term.
 
-Next we find the advective velocity through each interface.  
+Next we find the advective velocity through each interface.
 The incompressible
 velocity equation looks like the inviscid Burger's equation, and the
 Riemann solver follows that construction.  BCG provide the implementation
@@ -368,10 +368,10 @@ \subsection{Advective velocity}
 \end{equation}
 We solve this for each of the normal velocities, giving:
 \begin{eqnarray}
-\hat{u}^\mathrm{adv}_{i+\myhalf,j} &=& 
+\hat{u}^\mathrm{adv}_{i+\myhalf,j} &=&
     \mathcal{R}(\hat{u}_{i+\myhalf,j,L}^{n+\myhalf}, \hat{u}_{i+\myhalf,j,R}^{n+\myhalf}) \\
 %
-\hat{v}^\mathrm{adv}_{i,j+\myhalf} &=& 
+\hat{v}^\mathrm{adv}_{i,j+\myhalf} &=&
     \mathcal{R}(\hat{v}_{i,j+\myhalf,L}^{n+\myhalf}, \hat{v}_{i,j+\myhalf,R}^{n+\myhalf})
 \end{eqnarray}
 These advective velocities (sometimes called the {\em transverse
@@ -399,7 +399,7 @@ \subsection{Advective velocity}
 \hat{u}_{i,j+\myhalf} &=& \mathcal{U}[\hat{v}_{i,j+\myhalf}^\mathrm{adv}]
              (\hat{u}_{i,j+\myhalf,L}^{n+\myhalf}, \hat{u}_{i,j+\myhalf,R}^{n+\myhalf} ) \\
 \hat{v}_{i,j+\myhalf} &=& \mathcal{U}[\hat{v}_{i,j+\myhalf}^\mathrm{adv}]
-             (\hat{v}_{i,j+\myhalf,L}^{n+\myhalf}, \hat{v}_{i,j+\myhalf,R}^{n+\myhalf} ) 
+             (\hat{v}_{i,j+\myhalf,L}^{n+\myhalf}, \hat{v}_{i,j+\myhalf,R}^{n+\myhalf} )
 \end{eqnarray}
 
 Now we can construct the full left and right predictions for the
@@ -407,23 +407,23 @@ \subsection{Advective velocity}
 \ref{eq:vMAC_l}).  This involves simply adding the transverse term to
 the hat quantities and adding the pressure gradient.
 \begin{equation}
-u_{i+\myhalf,j,L}^{n+\myhalf} = \hat{u}_{i+\myhalf,j,L}^{n+\myhalf} 
-    - \frac{\Delta t}{2} 
+u_{i+\myhalf,j,L}^{n+\myhalf} = \hat{u}_{i+\myhalf,j,L}^{n+\myhalf}
+    - \frac{\Delta t}{2}
        \left [ \frac{1}{2} \left (\hat{v}^\mathrm{adv}_{i,j-\myhalf} +
                                   \hat{v}^\mathrm{adv}_{i,j+\myhalf} \right )
       \right ]
-      \left ( \frac{\hat{u}_{i,j+\myhalf}^{n+\myhalf} - 
+      \left ( \frac{\hat{u}_{i,j+\myhalf}^{n+\myhalf} -
                     \hat{u}_{i,j-\myhalf}^{n+\myhalf}}{\Delta y} \right )
     - \frac{\Delta t}{2} (Gp)^{(x),n-\myhalf}_{i,j}
 \end{equation}
 and
 \begin{equation}
-v_{i,j+\myhalf,L}^{n+\myhalf} = \hat{v}_{i,j+\myhalf,L}^{n+\myhalf} 
-    - \frac{\Delta t}{2} 
+v_{i,j+\myhalf,L}^{n+\myhalf} = \hat{v}_{i,j+\myhalf,L}^{n+\myhalf}
+    - \frac{\Delta t}{2}
        \left [ \frac{1}{2} \left (\hat{u}^\mathrm{adv}_{i-\myhalf,j} +
                                   \hat{u}^\mathrm{adv}_{i+\myhalf,j} \right )
       \right ]
-      \left ( \frac{\hat{v}_{i+\myhalf,j}^{n+\myhalf} - 
+      \left ( \frac{\hat{v}_{i+\myhalf,j}^{n+\myhalf} -
                     \hat{v}_{i-\myhalf,j}^{n+\myhalf}}{\Delta x} \right )
     - \frac{\Delta t}{2} (Gp)^{(y),n-\myhalf}_{i,j}
 \end{equation}
@@ -436,22 +436,22 @@ \subsection{Advective velocity}
 
 Finally, we do a Riemann solve (again, using the Burger's form of the
 Riemann problem) followed by upwinding to get the normal advective
-velocities.  This is basically the $\mathcal{R}(.,.)$ operation followed 
+velocities.  This is basically the $\mathcal{R}(.,.)$ operation followed
 by $\mathcal{U}(.,.)$.  Together, it is:
 \begin{equation}
 u^\mathrm{adv}_{i+\myhalf,j} = \left \{ \begin{array}{cl}
-   u_{i+\myhalf,j,L}^{n+\myhalf}  & 
-         \mathrm{~if~} u_{i+\myhalf,j,L}^{n+\myhalf} > 0, \qquad 
+   u_{i+\myhalf,j,L}^{n+\myhalf}  &
+         \mathrm{~if~} u_{i+\myhalf,j,L}^{n+\myhalf} > 0, \qquad
                        u_{i+\myhalf,j,L}^{n+\myhalf} + u_{i+\myhalf,j,R}^{n+\myhalf} > 0 \\
    \frac{1}{2} \left ( u_{i+\myhalf,j,L}^{n+\myhalf} +
                        u_{i+\myhalf,j,R}^{n+\myhalf} \right )  &
-          \mathrm{~if~} u_{i+\myhalf,j,L}^{n+\myhalf} \le 0, \qquad 
+          \mathrm{~if~} u_{i+\myhalf,j,L}^{n+\myhalf} \le 0, \qquad
                         u_{i+\myhalf,j,R}^{n+\myhalf} \ge 0 \\
    u_{i+\myhalf,j,R}^{n+\myhalf}  & \mathrm{~otherwise}
   \end{array}
   \right .
 \end{equation}
-and similar for $v^\mathrm{adv}_{i,j+\myhalf}$.  
+and similar for $v^\mathrm{adv}_{i,j+\myhalf}$.
 These velocities are sometimes referred to as the MAC velocities.
 
 
@@ -494,10 +494,10 @@ \subsection{MAC projection}
 \end{equation}
 using multigrid V-cycles and then update the MAC velocities as:
 \begin{eqnarray}
-u^\mathrm{adv}_{i+\myhalf,j} &=& u^\mathrm{adv}_{i+\myhalf,j} - 
+u^\mathrm{adv}_{i+\myhalf,j} &=& u^\mathrm{adv}_{i+\myhalf,j} -
     \frac{\phi_{i+1,j} - \phi_{i,j}}{\Delta x} \\
 %
-v^\mathrm{adv}_{i,j+\myhalf} &=& v^\mathrm{adv}_{i,j+\myhalf} - 
+v^\mathrm{adv}_{i,j+\myhalf} &=& v^\mathrm{adv}_{i,j+\myhalf} -
     \frac{\phi_{i,j+1} - \phi_{i,j}}{\Delta y}
 \end{eqnarray}
 
@@ -530,7 +530,7 @@ \subsection{Provisional update}
 We express the advective terms for $u$ as $A^{(u),n+\myhalf}_{i,j}$ and
 those for $v$ as $A^{(v),n+\myhalf}_{i,j}$.  These have the form:
 \begin{eqnarray}
-A^{(u),n+\myhalf}_{i,j} =   
+A^{(u),n+\myhalf}_{i,j} =
    && \frac{1}{2}
      \left ( u^\mathrm{adv}_{i-\myhalf,j} + u^\mathrm{adv}_{i+\myhalf,j} \right )
      \frac{u^{n+\myhalf}_{i+\myhalf,j} - u^{n+\myhalf}_{i-\myhalf,j}}{\Delta x} + \nonumber \\
@@ -538,7 +538,7 @@ \subsection{Provisional update}
      \left ( v^\mathrm{adv}_{i,j-\myhalf} + v^\mathrm{adv}_{i,j+\myhalf} \right )
      \frac{u^{n+\myhalf}_{i,j+\myhalf} - u^{n+\myhalf}_{i,j-\myhalf}}{\Delta y} \\
 %
-A^{(v),n+\myhalf}_{i,j} = 
+A^{(v),n+\myhalf}_{i,j} =
    && \frac{1}{2}
      \left ( u^\mathrm{adv}_{i-\myhalf,j} + u^\mathrm{adv}_{i+\myhalf,j} \right )
      \frac{v^{n+\myhalf}_{i+\myhalf,j} - v^{n+\myhalf}_{i-\myhalf,j}}{\Delta x} \nonumber +\\
@@ -551,7 +551,7 @@ \subsection{Provisional update}
 as:
 \begin{eqnarray}
 u^{\star}_{i,j} &=& u^n_{i,j} - \Delta t A^{(u),n+\myhalf}_{i,j} - \Delta t (G p)^{(x),n-\myhalf}_{i,j} \\
-v^{\star}_{i,j} &=& v^n_{i,j} - \Delta t A^{(v),n+\myhalf}_{i,j} - \Delta t (G p)^{(y),n-\myhalf}_{i,j} 
+v^{\star}_{i,j} &=& v^n_{i,j} - \Delta t A^{(v),n+\myhalf}_{i,j} - \Delta t (G p)^{(y),n-\myhalf}_{i,j}
 \end{eqnarray}
 Note that at this point, we don't have an updated $p$, so we use a lagged
 value from the previous step's projection.
@@ -561,7 +561,7 @@ \subsection{Provisional update}
 this update, giving an alternate provisional update:
 \begin{eqnarray}
 u^{\star\star}_{i,j} &=& u^n_{i,j} - \Delta t A^{(u),n+\myhalf}_{i,j}  \\
-v^{\star\star}_{i,j} &=& v^n_{i,j} - \Delta t A^{(v),n+\myhalf}_{i,j}  
+v^{\star\star}_{i,j} &=& v^n_{i,j} - \Delta t A^{(v),n+\myhalf}_{i,j}
 \end{eqnarray}
 
 
@@ -581,7 +581,7 @@ \subsection{Approximate projection}
 have the flexibility on whether to include the $G p^{n-\myhalf}$ term.  If
 we were doing an exact projection, then adding the gradient of a
 scalar would not change the divergence-free velocity field, so there
-would be no need to add it.  
+would be no need to add it.
 
 BCH did an exact projection on a cell-centered grid.  There, the
 divergence operator is:
@@ -590,7 +590,7 @@ \subsection{Approximate projection}
               \frac{v_{i,j+1} - v_{i,j-1}}{2\Delta y} \label{eq:ccdiv}
 \end{equation}
 This gives a cell-centered $D\Ub$.  If we want $\phi$ cell-centered, then
-the gradient, $G\phi$ must also be cell centered so $L\phi = DG\phi$ is 
+the gradient, $G\phi$ must also be cell centered so $L\phi = DG\phi$ is
 cell-centered.  This means that we must have
 \begin{eqnarray}
 (G\phi)^{(x)}_{i,j} &=& \frac{\phi_{i+1,j} - \phi_{i-1,j}}{2\Delta x} \\
@@ -602,7 +602,7 @@ \subsection{Approximate projection}
 (L\phi)_{i,j} = \frac{\phi_{i+2,j} -2\phi_{i,j} + \phi_{i-2,j}}{4 \Delta x^2} +
     \frac{\phi_{i,j+2} -2\phi_{i,j} + \phi_{i,j-2}}{4 \Delta y^2}
 \end{equation}
-This decomposes the domain into 4 distinct grids that are only linked 
+This decomposes the domain into 4 distinct grids that are only linked
 together at the boundaries.  While an exact projection, this decoupling
 can be undesirable.
 
@@ -610,7 +610,7 @@ \subsection{Approximate projection}
 projection, when you apply the projection operator, $P$, in
 succession, the result is unchanged ($P^2 = P$).  This is not the case
 for an approximate projection.  As a result, exactly what form you
-project matters.  For an approximate projection, we can use the 
+project matters.  For an approximate projection, we can use the
 standard 5-point stencil for the Laplacian,
 \begin{equation}
 (L\phi)_{i,j} = \frac{\phi_{i+1,j} -2\phi_{i,j} + \phi_{i-1,j}}{\Delta x^2} +
@@ -636,7 +636,7 @@ \subsection{Approximate projection}
 that has dimensions of pressure.  The procedure for the projection
 differs slightly depending on whether we project $\Ub^\star$ or $\Ub^{\star\star}$:
 \begin{itemize}
-\item case I: projecting $\Ub^{\star}/\Delta t$.  
+\item case I: projecting $\Ub^{\star}/\Delta t$.
 
 From the expression above,
 this looks like:
@@ -667,12 +667,12 @@ \subsection{Approximate projection}
 (see Martin 2.5 or ABC).  We store this for the next timestep.
 
 
-\item case II: projecting $\Ub^{\star\star}/\Delta t$.  
+\item case II: projecting $\Ub^{\star\star}/\Delta t$.
 
 From the expression above,
 this looks like:
 \begin{equation}
-\frac{\Ub^{\star\star}}{\Delta t} = \frac{\Ub^n}{\Delta t} - A^{(u),n+\myhalf} 
+\frac{\Ub^{\star\star}}{\Delta t} = \frac{\Ub^n}{\Delta t} - A^{(u),n+\myhalf}
 \end{equation}
 There is no explicit $Gp^{n-\myhalf}$ term.  We solve:
 \begin{equation}
@@ -685,7 +685,7 @@ \subsection{Approximate projection}
 \end{equation}
 
 Since there was no $Gp^{n-\myhalf}$ in what we projected, $p^{n+\myhalf} = \phi$,
-and 
+and
 \begin{equation}
 Gp^{n+\myhalf} = G\phi
 \end{equation}
@@ -701,7 +701,7 @@ \section{Boundary conditions}
 
 For the advection portion of the algorithm, the boundary conditions
 on $u$ and $v$ are implemented in the usual way, using ghost cells.
-For the projection, 
+For the projection,
 
 For a periodic domain, the boundary conditions on $\phi$ are likewise
 periodic.  At a solid wall or inflow boundary, we already predicted
@@ -711,7 +711,7 @@ \section{Boundary conditions}
 \begin{equation}
 \Ub^{n+1} = \Ub^\star - \nabla \phi
 \end{equation}
-we want $\nabla \phi \cdot n = 0$.  
+we want $\nabla \phi \cdot n = 0$.
 
 At outflow boundaries, we do not want to introduce any shear as we go
 through the boundary.  This means that we do not want any tangential
@@ -753,7 +753,7 @@ \subsection{Convergence test}
 solution.  This is run on a doubly-periodic unit square domain.  The
 main utility of this set of initial conditions is that we can use the
 analytic solution to measure the convergence behavior of the
-algorithm. 
+algorithm.
 
 
 \section{Extensions}
@@ -767,7 +767,7 @@ \section{Extensions}
   This can arise, for instance, in modeling the Rayleigh-Taylor
   instability.
 
-  The basic idea follows the method described above.  Now the 
+  The basic idea follows the method described above.  Now the
   mass continuity equation is also evolved:
   \begin{equation}
   \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \Ub) = 0
@@ -799,23 +799,23 @@ \section{Extensions}
   for $u$ and one for $v$).  These are differenced using
   Crank-Nicolson centering:
   \begin{eqnarray}
-  \frac{u^\star - u^n}{\Delta t} &=& - A^{(u),n+\myhalf} - \nabla p^{(x),n-\myhalf} 
+  \frac{u^\star - u^n}{\Delta t} &=& - A^{(u),n+\myhalf} - \nabla p^{(x),n-\myhalf}
        + \frac{\epsilon }{2} \nabla^2 (u^n + u^\star) \\
-  \frac{v^\star - v^n}{\Delta t} &=& - A^{(v),n+\myhalf} - \nabla p^{(y),n-\myhalf} 
-       + \frac{\epsilon }{2} \nabla^2 (v^n + v^\star) 
+  \frac{v^\star - v^n}{\Delta t} &=& - A^{(v),n+\myhalf} - \nabla p^{(y),n-\myhalf}
+       + \frac{\epsilon }{2} \nabla^2 (v^n + v^\star)
   \end{eqnarray}
   This involves two separate multigrid solves.  Once $\Ub^\star$ is found,
   the final projection is done as usual.
 
 \item {\em Low Mach number combustion}: In low-Mach number combustion
-  flows, the fluid is allowed to respond to local heat release 
+  flows, the fluid is allowed to respond to local heat release
   by expanding.  The velocity constraint is derived by differentiating
   the equation of state along particle paths, leading to the appearance
   of a source term:
   \begin{equation}
   \nabla \cdot \Ub = S
   \end{equation}
-  Here, $S$, incorporates the compressibility effects due to the 
+  Here, $S$, incorporates the compressibility effects due to the
   heat release and diffusion.  This system is used when modeling
   burning fronts (flames).  This type of flow can be thought
   of as linking two incompressible states (the fuel and the ash)
@@ -852,5 +852,3 @@ \section{Extensions}
   projections stumble on.
 
 \end{itemize}
-
-
diff --git a/intro/intro.tex b/intro/intro.tex
index d7ebd42..e66aa39 100644
--- a/intro/intro.tex
+++ b/intro/intro.tex
@@ -49,7 +49,7 @@ \section{What is simulation?}
 simulations, you should do parameter and convergence studies.
 
 There is no ``\"uber-code''.  Every algorithm begins with
-approximations and has limitations.  
+approximations and has limitations.
 Comparisons between different codes are important and common in our
 field (see, for example,
 \cite{frenk:1999,dimonte:2004,devalborro:2006}), and build confidence
@@ -99,7 +99,7 @@ \subsection{Sources of error}
 In your choice of programming language, create a floating point
 variable and initialize it to 0.1.  Now, print it out in full
 precision (you may need to use format specifiers in your language to
-get all significant digits the computer tracks).  
+get all significant digits the computer tracks).
 
 You should see that it is not exactly 0.1 to the computer---this is
 the floating point error.  The number 0.1 is not exactly representable
@@ -166,7 +166,7 @@ \subsection{Differentiation and integration}
   number of points, $x_0$, $x_1$, $\ldots$, and we want to compute the
   derivative at a point or integration over a range of points.
 \item We know a function analytically and we want to construct a
-  derivative or integral of this function numerically. 
+  derivative or integral of this function numerically.
 \end{enumerate}
 
  In these notes, we will mainly be concerned with the first case.
@@ -193,8 +193,8 @@ \subsubsection{Differentiation of discretely-sampled function}
 \end{equation}
 Solving for ${\partial a}/{\partial x} |_i$, we see
 \begin{align}
-\left . \frac{\partial a}{\partial x} \right |_i &= 
-    \frac{a_i - a_{i-1}}{\Delta x} - \frac{1}{2}\Delta x \left . \frac{\partial^2 a}{\partial x^2} \right |_i + \ldots \\ 
+\left . \frac{\partial a}{\partial x} \right |_i &=
+    \frac{a_i - a_{i-1}}{\Delta x} - \frac{1}{2}\Delta x \left . \frac{\partial^2 a}{\partial x^2} \right |_i + \ldots \\
 %
  &= \frac{a_i - a_{i-1}}{\Delta x} + \mathcal{O}(\Delta x)
 \end{align}
@@ -225,7 +225,7 @@ \subsubsection{Differentiation of discretely-sampled function}
 \caption[Difference approximations to the derivative of $\sin(x)$]
         {\label{fig:derivs} A comparison of one-sided and centered
           difference approximations to the derivative of $\sin(x)$.
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/basic_numerics/derivatives/derivatives.py}{derivatives.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/basic_numerics/derivatives/derivatives.py}{derivatives.py}}}
 \end{figure}
 
 Figure~\ref{fig:derivs} shows the left- and right-sided first-order
@@ -272,7 +272,7 @@ \subsubsection{Differentiation of an analytic function}
     \sin(x)$ at $x_0 = 1$ as a function of the spacing $\delta x$.  For
     small $\delta x$, roundoff error dominates the error in the
     approximation.  For large $\delta x$, truncation error dominates. \\
-    \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/basic_numerics/derivatives/deriv_error.py}{deriv\_error.py}}}
+    \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/basic_numerics/derivatives/deriv_error.py}{deriv\_error.py}}}
 \end{figure}
 
 
@@ -410,7 +410,7 @@ \subsection{Root finding}
   defines the next approximation to the root.  The vertical dotted
   line to the function shows the new slope that will be used for
   extrapolation in the next iteration. \\
-\hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/basic_numerics/roots/roots_plot.py}{roots\_plot.py}}}
+\hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/basic_numerics/roots/roots_plot.py}{roots\_plot.py}}}
 \end{figure}
 
 
@@ -447,15 +447,15 @@ \subsection{Norms}
 $i = 0, \ldots, N-1$.  The error at each point $i$ is $\epsilon_i =
 |f_i - f(x_i)|$.  But this is $N$ separate errors---we want a single
 number that represents the error of our approximation.  This is the
-job of a vector norm.  
+job of a vector norm.
 
 There are many different norms we can define.  For a vector ${\bf q}$,
 we write the norm as $\|q\|$.  Often, we'll put a subscript after the
 norm brackets to indicate which norm was used.  Some popular norms
 are:
 \begin{itemize}
-\item {\em inf norm}: 
-  \begin{equation} 
+\item {\em inf norm}:
+  \begin{equation}
     \| q \|_\infty = \max_i |q_i|
   \end{equation}
 
@@ -540,9 +540,9 @@ \subsection{ODEs}
 % rk4_plot.py creates rk4_k[1-4].pdf, rk4_final.pdf
 \begin{figure}[p]
 \centering
-\includegraphics[width=0.48\linewidth]{rk4_k1} 
+\includegraphics[width=0.48\linewidth]{rk4_k1}
 \includegraphics[width=0.48\linewidth]{rk4_k2} \\
-\includegraphics[width=0.48\linewidth]{rk4_k3} 
+\includegraphics[width=0.48\linewidth]{rk4_k3}
 \includegraphics[width=0.48\linewidth]{rk4_k4} \\
 \includegraphics[width=0.6\linewidth]{rk4_final}
 %
@@ -759,15 +759,15 @@ \subsection{FFTs}
   following transforms:
   \begin{itemize}
   \item $\sin(2\pi k_0 x)$ with $k_0 = 0.2$.  The transform should
-    have all of the power in the imaginary component only at the 
+    have all of the power in the imaginary component only at the
     frequency $0.2$.
 
   \item $\cos(2\pi k_0 x)$ with $k_0 = 0.2$.  The transform should
-    have all of the power in the real component only at the 
+    have all of the power in the real component only at the
     frequency $0.2$.
 
   \item $\sin(2\pi k_0 x + \pi/4)$ with $k_0 = 0.2$.  The transform should
-    have equal power in the real and imaginary components, only at the 
+    have equal power in the real and imaginary components, only at the
     frequency $0.2$.  Since the power is 1, the amplitude of the real and
     imaginary parts will be $1/\sqrt{2}$.
   \end{itemize}
@@ -802,7 +802,7 @@ \subsection{FFTs}
   imaginary parts.  The third pane shows the power (the absolute value
   of the transform).  Finally, the bottom panel shows the inverse
   transform of our transform, giving us back our original function. \\
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/basic_numerics/FFT/fft_simple_examples.py}{fft\_simple\_examples.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/basic_numerics/FFT/fft_simple_examples.py}{fft\_simple\_examples.py}}}
 \end{figure}
 
 
diff --git a/multigrid/multigrid.tex b/multigrid/multigrid.tex
index e75d48f..2ecc866 100644
--- a/multigrid/multigrid.tex
+++ b/multigrid/multigrid.tex
@@ -41,7 +41,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
 \nabla^2 \phi = f
 \end{equation}
 We will difference this in a second-order accurate fashion---see
-Figure~\ref{fig:mg:laplacian}.  
+Figure~\ref{fig:mg:laplacian}.
 \begin{figure}[t]
 \centering
 \includegraphics[width=\linewidth]{laplacian}
@@ -69,7 +69,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
 $\phi$ on edges:
 \begin{align}
 [\nabla \phi \cdot \hat{x}]_{i+\myhalf,j} &= \frac{\phi_{i+1,j} - \phi_{i,j}}{\dx} \\
-[\nabla \phi \cdot \hat{y}]_{i,j+\myhalf} &= \frac{\phi_{i,j+1} - \phi_{i,j}}{\dy} 
+[\nabla \phi \cdot \hat{y}]_{i,j+\myhalf} &= \frac{\phi_{i,j+1} - \phi_{i,j}}{\dy}
 \end{align}
 Again, since this is defined on edges, this represents a centered
 difference, and is therefore second-order accurate.  We then
@@ -84,7 +84,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
 %
   &= \frac{\phi_{i+1,j} - 2\phi_{i,j} + \phi_{i-1,j}}{\dx^2} +
      \frac{\phi_{i,j+1} - 2\phi_{i,j} + \phi_{i,j-1}}{\dy^2} = f_{i,j}
-\end{align}        
+\end{align}
 Again, since we used a centered-difference of the edge values, this
 expression is second-order accurate.  This is the standard {\em
   5-point stencil} for the 2-d Laplacian\footnote{There are other
@@ -112,7 +112,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
 \phi_{i,j} &= \frac{1}{MN} \sum_{k_x = 0}^{M-1} \sum_{k_y = 0}^{N-1}
   \Phi_{k_x,k_y} e^{2\pi\imag i k_x/M} e^{2\pi\imag j k_y/N} \\
 f_{i,j} &= \frac{1}{MN} \sum_{k_x = 0}^{M-1} \sum_{k_y = 0}^{N-1}
-  F_{k_x,k_y} e^{2\pi\imag i k_x/M} e^{2\pi\imag j k_y/N} 
+  F_{k_x,k_y} e^{2\pi\imag i k_x/M} e^{2\pi\imag j k_y/N}
 \end{align}
 
 Inserting these into the differenced equation, we have:
@@ -137,7 +137,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
   &\frac{\Phi_{k_x,k_x}}{\dx^2}
   \left [e^{2\pi\imag k_x/M} + e^{-2\pi\imag k_x/M} - 2 \right ] + \nonumber \\
   &\frac{\Phi_{k_x,k_x}}{\dy^2}
-  \left [e^{2\pi\imag k_y/N} + e^{-2\pi\imag k_y/N} - 2 \right ] 
+  \left [e^{2\pi\imag k_y/N} + e^{-2\pi\imag k_y/N} - 2 \right ]
   - F_{k_x,k_y}
   \biggr \} = 0
 \end{align}
@@ -183,7 +183,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
   Eq.~\ref{eq:mg:fftsource}. (right) Error vs.\ the true solution as a
   function of resolution for the Fourier method, showing second-order
   convergence.
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/elliptic/poisson_fft.py}{poisson\_fft.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/elliptic/poisson_fft.py}{poisson\_fft.py}}}
 \end{figure}
 
 
@@ -205,7 +205,7 @@ \section{\label{elliptic:sec:fft} Fourier Method}
 \section{Relaxation}
 
 Relaxation is an iterative technique, and as we will see shortly, it
-provides the basis for the multigrid technique.  
+provides the basis for the multigrid technique.
 
 Consider Poisson's equation, again differenced as:
 \begin{equation}
@@ -216,7 +216,7 @@ \section{Relaxation}
 1$ in $y$.  For the moment, consider the case where $\dx = \dy$.  If
 we solve this discretized equation for $\phi_{i,j}$, then we have:
 \begin{equation}
-\phi_{i,j} = \frac{1}{4} (\phi_{i+1,j} + \phi_{i-1,j} + 
+\phi_{i,j} = \frac{1}{4} (\phi_{i+1,j} + \phi_{i-1,j} +
                           \phi_{i,j+1} + \phi_{i,j-1} - \dx^2 f_{i,j} )
 \end{equation}
 A similar expression exists for every zone in our domain, coupling all
@@ -225,20 +225,20 @@ \section{Relaxation}
 technique called {\em relaxation} (also sometimes called {\em
   smoothing} because generally speaking the solution to elliptic
 equations is a smooth function) to find the solution for $\phi$
-everywhere.  
+everywhere.
 
 Imagine an initial guess to $\phi$: $\phi_{i,j}^{(0)}$.
 We can improve that guess by using our difference equation to define a
 new value of $\phi$, $\phi_{i,j}^{(1)}$:
 \begin{equation}
-\phi_{i,j}^{(1)} = \frac{1}{4} (\phi_{i+1,j}^{(0)} + \phi_{i-1,j}^{(0)} + 
-                                \phi_{i,j+1}^{(0)} + \phi_{i,j-1}^{(0)} - 
+\phi_{i,j}^{(1)} = \frac{1}{4} (\phi_{i+1,j}^{(0)} + \phi_{i-1,j}^{(0)} +
+                                \phi_{i,j+1}^{(0)} + \phi_{i,j-1}^{(0)} -
                                  \dx^2 f_{i,j} )
 \end{equation}
 or generally, the $k+1$ iteration will see:
 \begin{equation}
-\phi_{i,j}^{(k+1)} = \frac{1}{4} (\phi_{i+1,j}^{(k)} + \phi_{i-1,j}^{(k)} + 
-                                  \phi_{i,j+1}^{(k)} + \phi_{i,j-1}^{(k)} - 
+\phi_{i,j}^{(k+1)} = \frac{1}{4} (\phi_{i+1,j}^{(k)} + \phi_{i-1,j}^{(k)} +
+                                  \phi_{i,j+1}^{(k)} + \phi_{i,j-1}^{(k)} -
                                    \dx^2 f_{i,j} )
 \end{equation}
 This will (slowly) converge to the true solution\footnote{Formally,
@@ -256,8 +256,8 @@ \section{Relaxation}
 will see a mix of the old and new solutions.  We can express this in-place
 updating as:
 \begin{equation}
-\phi_{i,j} \leftarrow \frac{1}{4} (\phi_{i+1,j} + \phi_{i-1,j} + 
-                                   \phi_{i,j+1} + \phi_{i,j-1} - 
+\phi_{i,j} \leftarrow \frac{1}{4} (\phi_{i+1,j} + \phi_{i-1,j} +
+                                   \phi_{i,j+1} + \phi_{i,j-1} -
                                    \dx^2 f_{i,j} )
 \end{equation}
 This only requires a single copy of $\phi$ to be stored.  This
@@ -272,7 +272,7 @@ \section{Relaxation}
 \end{equation}
 We can discretize this as:
 \begin{equation}
-\alpha \phi_{i,j} - \beta \left ( 
+\alpha \phi_{i,j} - \beta \left (
     \frac{\phi_{i+1,j} - 2 \phi_{i,j} + \phi_{i-1,j}}{\dx^2} +
     \frac{\phi_{i,j+1} - 2 \phi_{i,j} + \phi_{i,j-1}}{\dy^2} \right )
 = f_{i,j}
@@ -283,10 +283,10 @@ \section{Relaxation}
      \left . \left ( f_{i,j} + \frac{\beta}{\dx^2} \phi_{i+1,j}
                              + \frac{\beta}{\dx^2} \phi_{i-1,j}
                              + \frac{\beta}{\dy^2} \phi_{i,j+1}
-                             + \frac{\beta}{\dy^2} \phi_{i,j-1} \right ) 
+                             + \frac{\beta}{\dy^2} \phi_{i,j-1} \right )
 \middle / \left ( \alpha + \frac{2 \beta}{\dx^2} + \frac{2 \beta}{\dy^2} \right ) \right .
 \end{equation}
-Notice that if $\alpha = 0$, $\beta = -1$, and $\dx = \dy$, we 
+Notice that if $\alpha = 0$, $\beta = -1$, and $\dx = \dy$, we
 recover the relaxation expression for Poisson's equation from above.
 
 
@@ -305,9 +305,9 @@ \subsection{Boundary conditions}
 \begin{figure}[h]
 \centering
 \includegraphics[width=\linewidth]{fv-fd_grid_bc}
-\caption[Node-centered vs.\ cell-centered data at boundaries]{\label{mg:fig:bcs} 
+\caption[Node-centered vs.\ cell-centered data at boundaries]{\label{mg:fig:bcs}
   A node-centered (top) and cell-centered (bottom) finite difference
-  grid showing the data and domain boundaries.  Notice that for the 
+  grid showing the data and domain boundaries.  Notice that for the
   cell-centered grid, there is no data point precisely at the boundary.}
 \end{figure}
 
@@ -318,7 +318,7 @@ \subsection{Boundary conditions}
     boundary condition}} We'll label the first zone inside the domain,
 at the left boundary, $\mathrm{lo}$, and the last zone inside the
 domain, at the right boundary, $\mathrm{hi}$---see
-Figure~\ref{mg:fig:bcs}.  To second order, we can average the zone values on 
+Figure~\ref{mg:fig:bcs}.  To second order, we can average the zone values on
 either side of the interface to get the boundary condition:
 \begin{eqnarray}
 \phi_l &=& \frac{1}{2} ( \phi_\mathrm{lo} + \phi_\mathrm{lo-1} ) \\
@@ -354,7 +354,7 @@ \subsection{Residual and true error}
 satisfies the discretized equation.  For the Poisson equation, we
 can the residual as:
 \begin{equation}
-r_{i,j} = f_{i,j} - (L \phi)_{i,j} 
+r_{i,j} = f_{i,j} - (L \phi)_{i,j}
 \end{equation}
 and the residual error as:
 \begin{equation}
@@ -366,7 +366,7 @@ \subsection{Residual and true error}
 the true solution.  If $\phi^\mathrm{true}$ satisfies $\nabla^2
 \phi^\mathrm{true} = f$, then the true error in each zone is
 \begin{equation}
-e_{i,j} = \phi^\mathrm{true}(x_i,y_j) - \phi_{i,j} 
+e_{i,j} = \phi^\mathrm{true}(x_i,y_j) - \phi_{i,j}
 \end{equation}
 and
 \begin{equation}
@@ -396,9 +396,9 @@ \subsection{Norms}
 \begin{equation}
 \|e\|_\infty = \max_{i,j} \{ |e_{i,j}| \}
 \end{equation}
-This will pick up on local errors.  
+This will pick up on local errors.
 
-The $L_1$ norm and $L_2$ norms are more global. 
+The $L_1$ norm and $L_2$ norms are more global.
 \begin{eqnarray}
 \|e\|_1 &=& \frac{1}{N} \sum_{i,j} |e_{i,j} | \\
 \|e\|_2 &=& \left ( \frac{1}{N} \sum_{i,j} |e_{i,j} |^2 \right )^{1/2}
@@ -423,7 +423,7 @@ \subsection{Performance}
   L2 norm of the error compared with the true solution (solid lines)
   and the L2 norm of the residual (dotted lines) for 3 different
   resolutions (16, 32, and 64 zones). \\
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/smooth.py}{smooth.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/smooth.py}{smooth.py}}}
 \end{figure}
 
 Consider the simple Poisson problem\footnote{This is the test problem
@@ -433,7 +433,7 @@ \subsection{Performance}
 \begin{equation}
 \phi_{xx} = \sin(x), \qquad \phi(0) = \phi(1) = 0
 \end{equation}
-The analytic solution to this is simply 
+The analytic solution to this is simply
 \begin{equation}
 \phi^a(x) = -\sin(x) + x \sin(1)
 \end{equation}
@@ -465,11 +465,11 @@ \subsection{Performance}
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{smooth-error-norms}
-\caption[Convergence of smoothing in different norms]{\label{fig:smoothnorms} 
+\caption[Convergence of smoothing in different norms]{\label{fig:smoothnorms}
   Gauss-Seidel relaxation applied to $\phi_{xx} = \sin(x)$ with
   $\phi(0) = \phi(1) = 0$.  This is like figure \ref{fig:smootherror}, but
   now we show the error in 3 different norms: $L_1$, $L_2$, and $L_\infty$.
-  \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/smooth-norms.py}{smooth-norms.py}}}
+  \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/smooth-norms.py}{smooth-norms.py}}}
 \end{figure}
 
 To demonstrate the influence of the boundary conditions,
@@ -479,7 +479,7 @@ \subsection{Performance}
 boundary value, instead of averaging to the interface:
 \begin{align}
 \phi_\mathrm{lo-1} &= \phi_\mathrm{lo} \\
-\phi_\mathrm{hi+1} &= \phi_\mathrm{hi} 
+\phi_\mathrm{hi+1} &= \phi_\mathrm{hi}
 \end{align}
 We see that with this mistake at the boundaries, the error of the
 entire solution is affected, and we get only first-order convergence
@@ -493,13 +493,13 @@ \subsection{Performance}
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{smooth-badBCs}
-\caption[Convergence of smoothing in first-order BCs]{\label{fig:smooth-badbcs} 
+\caption[Convergence of smoothing in first-order BCs]{\label{fig:smooth-badbcs}
   The same problem as in figure \ref{fig:smootherror}, but
   now we done with a naive first-order boundary conditions---just initializing
-  the ghost cell to the boundary value (solid lines).  We see that this 
+  the ghost cell to the boundary value (solid lines).  We see that this
   achieves only first-order convergence in the true error.  The correct
   second-order implmentation is shown as the dotted lines.
-  \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/smooth-badbcs.py}{smooth-badbcs.py}}}
+  \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/smooth-badbcs.py}{smooth-badbcs.py}}}
 \end{figure}
 
 
@@ -510,7 +510,7 @@ \subsection{Performance}
   Error in the solution to $\phi'' = 0$ given an initial guess with 3
   different wavenumbers of noise.  A 128 zone grid was used.  The
   different curves are different numbers of smoothing
-  iterations. \\ \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/smooth-modes.py}{smooth-modes.py}}}
+  iterations. \\ \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/smooth-modes.py}{smooth-modes.py}}}
 \end{figure}
 
 \subsection{Frequency/wavelength-dependent error}
@@ -555,11 +555,11 @@ \subsection{Frequency/wavelength-dependent error}
 the Laplace equation, $\phi^{\prime\prime} = 0$ in 1-d on $[0,1]$ with
 homogeneous Dirichlet boundary conditions.  The solution is simply
 $\phi = 0$.  We use Eq.~\ref{eq:mg:phimodes} as an initial guess---this is a superposition
-of 3 different modes.  We see that the after just a few iterations, the 
+of 3 different modes.  We see that the after just a few iterations, the
 shortest wavelength mode, $\sin(2\pi 16x)$ is no longer apparent in the error,
 and just the two longer wavelength modes dominate.  By 100 iterations,
 the error appears to be only due to the longest wavelength mode, $\sin(2\pi x)$.
-Even after 1000 iterations, we still see this mode in the error.  This 
+Even after 1000 iterations, we still see this mode in the error.  This
 demonstrates that the longest wavelength modes in the error take the
 longest to smooth away.
 
@@ -579,7 +579,7 @@ \section{Multigrid}
 \subsection{Prolongation and restriction on cell-centered grids}
 
 Multigrid relies on transferring the problem up and down a hierarchy of
-grids. 
+grids.
 The process of moving the data from the fine grid to the coarse grid
 is called {\em restriction}.  The reverse process, moving data from
 the coarse grid to the fine grid is called {\em prolongation}.  If the
@@ -607,21 +607,21 @@ \subsubsection{1-d}
 \begin{equation}
 \phi^c_j = \frac{1}{\dx^c} (\dx^f \phi^f_i + \dx^f \phi^f_{i+1} )
 \end{equation}
-or, for a jump in 2 in resolution ($\dx^c = 2 \dx^f$), 
+or, for a jump in 2 in resolution ($\dx^c = 2 \dx^f$),
 \begin{equation}
 \phi^c_j = \frac{1}{2} (\phi^f_i + \phi^f_{i+1} )
 \end{equation}
 This latter form appears as a simple average to the interface of the
 two fine cells / center of the corresponding coarse cell.
 
-The simplest type of 
+The simplest type of
 prolongation is simply {\em direct injection}:
 \begin{align}
 \phi^f_i &= \phi^c_j \\
-\phi^f_{i+1} &= \phi^c_j 
+\phi^f_{i+1} &= \phi^c_j
 \end{align}
 A higher-order method is to do linear reconstruction of the coarse
-data and then average under the profile, e.g., 
+data and then average under the profile, e.g.,
 \begin{equation}
 \phi(x) = \frac{\phi^c_{i+1} - \phi^c_{i-1}}{2\dx^c} (x - x_i^c) + \phi_i^c
 \end{equation}
@@ -636,7 +636,7 @@ \subsubsection{1-d}
 giving
 \begin{align}
 \phi^f_i     &= \phi^c_i - \frac{1}{8} (\phi^c_{i+1} - \phi^c_{i-1}) \\
-\phi^f_{i+1} &= \phi^c_i + \frac{1}{8} (\phi^c_{i+1} - \phi^c_{i-1}) 
+\phi^f_{i+1} &= \phi^c_i + \frac{1}{8} (\phi^c_{i+1} - \phi^c_{i-1})
 \end{align}
 Notice that this is conservative, since $\dx^f (\phi^f_i + \phi^f_{i+1}) = \dx^c \phi^c_j$.
 
@@ -667,7 +667,7 @@ \subsubsection{2-d}
 this reconstruction to determine what the fine cell values are.  For
 instance, a linear reconstruction of the coarse data in $x$ and $y$ is:
 \begin{equation}
-\phi(x,y) = \frac{m_x}{\dx^c} (x - x_i^c) + 
+\phi(x,y) = \frac{m_x}{\dx^c} (x - x_i^c) +
             \frac{m_y}{\dy^c} (y - y_j^c) + \phi_{i,j}^c
 \end{equation}
 with slopes:
@@ -719,9 +719,9 @@ \subsection{Multigrid cycles}
 error equations are now homogeneous.  This must be understood by
 any ghost cell filling routines.
 
-There are many different forms of the multigrid process.  The simplest 
+There are many different forms of the multigrid process.  The simplest
 is called the {\em V-cycle}.  Here you start of the fine grid, restrict
-down to the coarsest, solve, and then prolong back up to the finest. 
+down to the coarsest, solve, and then prolong back up to the finest.
 The flow looks like a `V'.  You continue with additional V-cycles
 until the residual error is smaller than your tolerance.
 
@@ -797,7 +797,7 @@ \subsection{Boundary conditions throughout the hierarchy}
 Now the left side looks precisely like the differenced Poisson equation
 with homogeneous Dirichlet BCs.  The RHS has an additional `charge' that
 captures the boundary value.  By modifying the source term, $f$, in the
-multigrid solver to include this charge, we can use the homogeneous 
+multigrid solver to include this charge, we can use the homogeneous
 ghost cell filling routines throughout the multigrid algorithm.
 This technique is discussed a bit in~\cite{colellanotes}.
 
@@ -818,7 +818,7 @@ \subsection{Stopping criteria}
 $\|f\|$ is called the {\em source norm}.  If $\|f\| = 0$, then we stop
 when
 \begin{equation}
-\| r \| < \epsilon 
+\| r \| < \epsilon
 \end{equation}
 Picking the tolerance $\epsilon$ is sometimes problem-dependent, and
 generally speaking, a problem with a large number of zones will require
@@ -841,7 +841,7 @@ \subsection{Stopping criteria}
   see that the true error, $\|e\|$ stalls at truncation error while
   the residual error, $\|r\|$ reaches roundoff error, the same
   behavior as seen with smoothing alone (as expected). \\
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/mg_test.py}{mg\_test.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/mg_test.py}{mg\_test.py}}}
 \end{figure}
 
 The overall convergence of the multigrid algorithm is limited by the
@@ -855,8 +855,8 @@ \subsection{Stopping criteria}
 \includegraphics[width=0.85\linewidth]{mg-converge}
 \caption[Convergence of the multigrid
   algorithm]{\label{fig:mg:convergence} Convergence of the multigrid
-  algorithm. \\ 
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multigrid/mg_converge.py}{mg\_converge.py}}}
+  algorithm. \\
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multigrid/mg_converge.py}{mg\_converge.py}}}
 \end{figure}
 
 
@@ -882,7 +882,7 @@ \section{Solvability}
 \begin{equation}
 \phi^{\prime\prime} = 0
 \end{equation}
-on $[a, b]$ with 
+on $[a, b]$ with
 \begin{align}
 \phi^\prime(a) &= A \\
 \phi^\prime(b) &= B
diff --git a/multiphysics/multiphysics.tex b/multiphysics/multiphysics.tex
index 4149c02..fa0371c 100644
--- a/multiphysics/multiphysics.tex
+++ b/multiphysics/multiphysics.tex
@@ -62,7 +62,7 @@ \section{Ex: diffusion-reaction}
 where $R_{\Delta t/2}$ represents reacting for a step of $\Delta t/2$
 and $D_{\Delta t}$ represents diffusing for a step of $\Delta t$.  In
 each case, these operators act as if the other were not present, but
-they see the effect of the previous operation on the input $\phi$.  
+they see the effect of the previous operation on the input $\phi$.
 Note that no explicit source terms describing one process appear in the other
 process's update.  The procedure for updating appears as:
 \begin{enumerate}
@@ -70,7 +70,7 @@ \section{Ex: diffusion-reaction}
 
   Define $\phi^\star$ as the the solution to the ODE:
    \begin{equation}
-     \frac{d\phi}{dt} = \frac{1}{\tau} R(\phi), \quad 
+     \frac{d\phi}{dt} = \frac{1}{\tau} R(\phi), \quad
    \phi(t^n) = \phi^n, \quad t \in [t^n, t^{n+\myhalf}]
    \end{equation}
 
@@ -88,7 +88,7 @@ \section{Ex: diffusion-reaction}
 
   Define $\phi^{n+1}$ as the the solution to the ODE:
    \begin{equation}
-     \text{define } \phi^{n+1}: \, \frac{d\phi}{dt} = \frac{1}{\tau} R(\phi), \quad 
+     \text{define } \phi^{n+1}: \, \frac{d\phi}{dt} = \frac{1}{\tau} R(\phi), \quad
    \phi(t^{n+\myhalf}) = \phi^{\star\star}, \quad t \in [t^{n+\myhalf}, t^{n+1}]
    \end{equation}
 
@@ -113,7 +113,7 @@ \section{Ex: diffusion-reaction}
    with 256 zones, and $\kappa = 0.1$, $\tau = 1.0$.  The lines shown
    are spaced 8.0 time-units apart.  We see the initial smoothed tophat
    profile giving rise to a traveling front. \\
-   \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multiphysics/diffusion-reaction.py}{diffusion-reaction.py}}
+   \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multiphysics/diffusion-reaction.py}{diffusion-reaction.py}}
    }
 \end{figure}
 
@@ -122,10 +122,10 @@ \section{Ex: diffusion-reaction}
  some details of this system).
 
 \begin{exercise}[Diffusion-reaction system]
- Solve the the difusion reaction equation with the source given in 
+ Solve the the difusion reaction equation with the source given in
  Eq.~\ref{eq:multiphyiscs:react}.  Note that you should begin with
- some smooth initial conditions---if they are too sharp than the 
- C-N discretization will cause jagged features to appear.  
+ some smooth initial conditions---if they are too sharp than the
+ C-N discretization will cause jagged features to appear.
  Compare your results to Figure~\ref{fig:diffreact}.
 \end{exercise}
 
@@ -144,36 +144,36 @@ \section{Ex: advection-diffusion}
 a physical width.
 
 As we saw earlier, there are efficient, accurate methods for handling
-the advective parts explicitly, but for diffusion, we often want to 
-solve it implicitly.  We can split the solution up, but couple the 
+the advective parts explicitly, but for diffusion, we often want to
+solve it implicitly.  We can split the solution up, but couple the
 two processes together to make a method that is overall second-order
 accurate in time.  We write our equation as:
 \begin{equation}
 u_t + A(u) = D(u)
 \end{equation}
-with $A(u) = [\frac{1}{2} u^2]_x$ and $D(u) = \epsilon u_{xx}$.  Then our update 
+with $A(u) = [\frac{1}{2} u^2]_x$ and $D(u) = \epsilon u_{xx}$.  Then our update
 appears in two steps.
 \begin{enumerate}
 \item {\em Find the advective update over the timestep}:
    We use an approximation of the diffusion term at time-level $n$, $D(u^n)$
-   as a source in the construction of the interface states for the 
+   as a source in the construction of the interface states for the
    advective part.  Once the interface states, $u_{i+\myhalf}^{n+\myhalf}$ are
    known, we construct the advective update term as:
    \begin{equation}
-   A_i^{n+\myhalf} = 
+   A_i^{n+\myhalf} =
      \frac{\left [ \frac{1}{2} \left (u_{i+\myhalf}^{n+\myhalf}\right)^2\right ] -
            \left [ \frac{1}{2} \left (u_{i-\myhalf}^{n+\myhalf}\right)^2\right ]}
           {\Delta x}
     \end{equation}
 
 \item {\em Solve the diffusion equation with the advective source}:
-    We use a Crank-Nicolson discretization of the diffusion part of 
+    We use a Crank-Nicolson discretization of the diffusion part of
     our equation, with the advective update term appearing as a source.
     \begin{equation}
-    \frac{u^{n+1} - u^n}{\Delta t} = 
+    \frac{u^{n+1} - u^n}{\Delta t} =
         \frac{1}{2}D(u^n) + \frac{1}{2}D(u^{n+1}) - A^{n+\myhalf}
     \end{equation}
-    
+
     This is a linear system that can be solved as a tridiagonal matrix
     or with multigrid.  The result of this step is that $u^{n+1}$ is
     updated with both the advection and diffusion terms.
@@ -191,7 +191,7 @@ \section{Ex: advection-diffusion}
   with a variety of different viscosities.  The initial conditions was
   a single wavelength of a sine wave for $x \in [1/3,2/3]$, and $u = 1$
   otherwise. \\
-  \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multiphysics/burgersvisc.py}{burgersvisc.py}}}
+  \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multiphysics/burgersvisc.py}{burgersvisc.py}}}
 \end{figure}
 
 For step 1, the addition of the explicit diffusion source requires
@@ -227,7 +227,7 @@ \subsection{Convergence without an analytic solution}
 compare coarser simulations to a coarsened version of the benchmark solution.
 This is done in Figure~\ref{fig:multiphysics:burgersconverge}.  The
 fine zones of the benchmark are averaged into a corresponding coarse grid zone
-to produce a coarsened representation of the benchmark, and the norm 
+to produce a coarsened representation of the benchmark, and the norm
 of the error with the coarse simulation is computed.
 
 \begin{figure}[t]
@@ -238,6 +238,6 @@ \subsection{Convergence without an analytic solution}
     solution with $\epsilon = 0.005$.  To test convergence a run with 4096
     zones was done as a benchmark solution.  We see second-order
     convergence for the higher-resolution runs in this sample.\\
-   \hydroexdoit{\href{https://github.com/zingale/hydro_examples/blob/master/multiphysics/burgersvisc_converge.py}{burgersvisc\_converge.py}}
+   \hydroexdoit{\href{https://github.com/python-hydro/hydro_examples/blob/master/multiphysics/burgersvisc_converge.py}{burgersvisc\_converge.py}}
    }
 \end{figure}
diff --git a/preface/preface.tex b/preface/preface.tex
index 889566c..bbbb6da 100644
--- a/preface/preface.tex
+++ b/preface/preface.tex
@@ -49,7 +49,7 @@
   script will be noted in the figure caption.
 
   You can get this set of scripts from github at:\\
-  \url{https://github.com/zingale/hydro_examples/}
+  \url{https://github.com/python-hydro/hydro_examples/}
 
   References to the scripts in \hydroex\ are shown throughout
   the text as: \\[0.5em]
@@ -79,7 +79,7 @@
   hydrodynamics code that implements the piecewise parabolic method
   from Chapter~\ref{ch:compressible}.  It can be obtained from
   github at:\\
-  \url{https://github.com/zingale/hydro1d/}
+  \url{https://github.com/python-hydro/hydro1d/}
 
   Details on it are given in
   Appendix~\ref{app:hydro1d}.