diff --git a/.gitignore b/.gitignore
index 163d09fd929ab0a4e414915592c68e0293c11eb8..f843002a1c03af25e89a601c75ac14f1203cca84 100644
--- a/.gitignore
+++ b/.gitignore
@@ -185,6 +185,7 @@ theory/Multipoles/potential_short.pdf
theory/Multipoles/force_short.pdf
theory/Cosmology/cosmology.pdf
theory/Cooling/eagle_cooling.pdf
+theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.pdf
m4/libtool.m4
m4/ltoptions.m4
diff --git a/theory/Gizmo/gizmo-implementation-details/bert.tex b/theory/Gizmo/gizmo-implementation-details/bert.tex
new file mode 100644
index 0000000000000000000000000000000000000000..99975349d90ef99f736c37dd3a4ec8d2092189cb
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/bert.tex
@@ -0,0 +1,104 @@
+
+\section{GIZMO~: variables and basic equations}
+
+This section is left unchanged and as the first in this document such that the comments in the code still correspond to the correct equations.
+
+
+\subsection{Variables}
+6 standard variables~: $\vec{x}$, $\vec{v}$\\ \\
+5 conserved variables ($C$)~: $m$, $\vec{p}$, $E$\\ \\
+5 primitive variables ($X$)~: $\rho$, $\vec{w}$, $P$\\
+It could be that $\vec{w} = \vec{v}$, but this is not generally true for all methods...\\ \\
+6 matrix elements + the volume (can be computed from $m$ and $\rho$ if this is
+accurate enough)\\ \\
+15 gradients (3 for every primitive variable)\\
+
+In total: 37 quantities per particle (without the volume), 34 if the velocity is equal to the fluid velocity
+
+\subsection{Equations by Neighbour Loops}
+
+\subsubsection{Loop 1: volumes (=densities) + matrix elements}
+Volumes are given by (combination of eqns. (7) and (27) from Hopkins)~:
+\begin{equation}
+V_i = \frac{1}{\sum_j W(|\vec{x_i}-\vec{x_j}|, h_i)}
+\end{equation}
+
+The 6 elements of the (symmetric) matrix $E_i$ by (eqn. (14) in Hopkins without the normalization, because $E_i$ is only used in combination with $\psi_i$ (eqn. (12)), which contains the same normalization and cancels it out again)~:
+\begin{equation}
+E_i^{\alpha \beta} = \sum_j (\vec{x_j}-\vec{x_i})^\alpha (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i)
+\end{equation}
+
+The kernel functions in SWIFT do not include the normalization factor $\frac{1}{h_i^3}$, so this is added after the loop in the ghost.
+
+\subsubsection{Ghost 1: primitive variables + matrix inversion}
+Use the volume to convert conserved variables to primitive variables (eqn. (5), (7) and (9) in the AREPO paper)~:
+\begin{align}
+\rho &= \frac{m}{V}\\
+\vec{v} &= \frac{\vec{p}}{m}\\
+P &= \frac{\gamma - 1}{V}\left(E - \frac{1}{2} \frac{|\vec{p}|^2}{m}\right)
+\end{align}
+
+Invert the matrix and call the inverse matrix $B_i$ (eqn. (13) in Hopkins).
+
+\subsubsection{Loop 2: gradients}
+Calculate gradients for the primitive variables (eqn. (12) in Hopkins with eqn. (6) inserted and the normalization constant eliminated)~:
+\begin{equation}
+ \left(\vec{\nabla} X_i\right)^\alpha = \sum_j \left(X_j - X_i\right) \sum_\beta B_i^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i)
+\end{equation}
+If you want to use a slope limiter of some sorts, then this would also be the time to collect the necessary data.
+
+\subsubsection{Ghost 2: nothing?}
+Perform the slope limiting if wanted.
+
+\subsubsection{Loop 3: hydro}
+For every neighbour $j$, calculate an interface area (combination of eqn. (12) in Hopkins and the definition of $\vec{A_{ij}}$, given in between eqns. (18) and (19))
+\begin{equation}
+ \left(\vec{A_{ij}}\right)^\alpha = V_i \sum_\beta B_i^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_i) + V_j \sum_\beta B_j^{\alpha \beta} (\vec{x_j}-\vec{x_i})^\beta W(|\vec{x_i}-\vec{x_j}|, h_j)
+\end{equation}
+
+Calculate a position for the interface (eqn. (20) in Hopkins)~:
+\begin{equation}
+ \vec{x_{ij}} = \vec{x_i} + \frac{h_i}{h_i+h_j} (\vec{x_j} - \vec{x_i})
+\end{equation}
+and a velocity (eqn. (21) in Hopkins)~:
+\begin{equation}
+ \vec{v_{ij}} = \vec{v_i} + (\vec{v_j}-\vec{v_i})\left[ \frac{(\vec{x_{ij}}-\vec{x_i}).(\vec{x_j}-\vec{x_i})}{|\vec{x_i}-\vec{x_j}|^2} \right]
+\end{equation}
+
+Use the velocity to boost the primitive variables $X_k$ ($k=i,j$) to the rest frame of the interface ($X'_k$). This comes down to applying a correction to the fluid velocities.
+
+Use the gradients to interpolate the primitive variables from positions $\vec{x_i}$ and $\vec{x_j}$ to $\vec{x_{ij}}$ and predict them forward in time
+by half a time step to obtain second order accuracy in time (eqn. (A3) in Hopkins)~:
+\begin{equation}
+ X''_k = X'_k + \left(\vec{\nabla}X\right)_k.(\vec{x_{ij}}-\vec{x_k}) + \frac{\Delta t}{2} \frac{\partial X'_k}{\partial t}.
+\end{equation}
+The time derivatives are given by (eqn. (A4) in Hopkins)
+\begin{equation}
+ \frac{\partial X'_k}{\partial t} = -\begin{pmatrix}
+ \vec{w'} & \rho' & 0\\
+ 0 & \vec{w'} & 1/\rho'\\
+ 0 & \gamma P' & \vec{w'}
+ \end{pmatrix} \vec{\nabla}X_k,
+\end{equation}
+or, for example~:
+\begin{equation}
+ \frac{\rho_k^{n+1/2}-\rho_k^n}{\Delta t/2} = -\left( \vec{w_k}-\vec{v_{ij}}^n \right).\vec{\nabla}\rho_k^n - \rho_k^n \vec{\nabla}.\vec{w_k}^n
+\end{equation}
+
+We then feed the $X''_k$ to a 1D Riemann solver, that gives us a vector of fluxes $\vec{F''_{ij}}$, which we deboost to a static frame of reference ($\vec{F_{ij}}$). To account for the arbitrary orientation of the interface, we also feed the normal vector to the interface to the Riemann solver. The Riemann solver internally uses the fluid velocity along the interface normal to solve an effective 1D Riemann problem. The velocity solution along the interface normal is also internally added to the velocities (see GIZMO source code). Mathematically, this is the same as first rotating the velocities to a frame aligned with the interface and then rotating the solution back to the original frame (e.g. AREPO paper), but it is computationally cheaper and causes less round off error in the solution.
+
+Given the solution vector $X_\text{half}$ returned by the Riemann problem, the fluxes are given by (eqns. (1)-(3) in Hopkins)
+\begin{align}
+ \vec{F}_\rho &=
+ \rho_\text{half} \left( \vec{w}_\text{half} - \vec{v}_{ij}\right)\\
+ \vec{F}_{w_k} &=
+ \rho_\text{half} w_{k, \text{half}} \left( \vec{w}_\text{half} - \vec{v}_{ij} \right) + P_\text{half} \vec{n}_k\\
+ \vec{F}_P &=
+ \rho_\text{half} e_\text{half} \left( \vec{w}_\text{half} - \vec{v}_{ij}\right) + P_\text{half} \vec{w}_\text{half},
+\end{align}
+with $\vec{n}_k$ the unit vector along the coordinate axes ($k=x,y,z$) and $e = \frac{P_\text{half}}{\left(\gamma -1\right)\rho_\text{half}} + \frac{1}{2}\vec{w}_\text{half}^2$.
+
+Finally, we use the fluxes to update the conserved variables (eqn. (23) in Hopkins)~:
+\begin{equation}
+ \Delta C_k = -\Delta t \sum_j \vec{F_{ij}}.\vec{A_{ij}}.
+\end{equation}
\ No newline at end of file
diff --git a/theory/Gizmo/gizmo-implementation-details/computations.tex b/theory/Gizmo/gizmo-implementation-details/computations.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3d1ba66f6ace87a8f2615c632294aec24d44b85a
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/computations.tex
@@ -0,0 +1,282 @@
+%=========================================
+\section{Explicit Computations}
+%=========================================
+
+Main purpose of this section is to have the explicit computations written down clearly somewhere.
+Implicitly swallowing indices in the formulae can make programming life a nightmare, so here we go.
+
+
+
+%=========================================
+\subsection{Normalization}
+%=========================================
+
+
+To compute the normalisations \ref{omega}, we need to sum over all neighbouring particles and sum the kernels correctly.
+To evaluate the kernels, we use the \verb|kernel_deval(xij, wij, wij_dx)| function in SWIFT.
+
+If a kernel is defined as
+\begin{align*}
+ W_i(\x) = W(\x - \x_i, h(\x)) = \frac{1}{h(\x)^\nu} w\left(\frac{| \x - \x_i |}{h(\x)} \right)
+\end{align*}
+
+then \verb|kernel_deval| computes
+
+\begin{align*}
+ \mathtt{wij} &= w(\mathtt{xij}) \\
+ \text{and \quad} \mathtt{wij\_dx} &= \frac{\del w(r)}{\del r} \big{|} _{r = \mathtt{xij}} \\
+ \text{with \quad} \mathtt{xij} &= \frac{| \x_i - \x_j |}{h(\x_i)}
+\end{align*}
+
+
+
+So for a specific particle position $i$, we need to compute
+\begin{align*}
+ \omega(\x_i) &= \sum_j W(\x_i - \x_j, h(\x_i)) \\
+ &= \sum_j \frac{1}{h(\x_i)^\nu} w\left( \frac{| \x_i - \x_j |}{h(\x_i)} \right) \\
+ &= \sum_j \frac{1}{h_i^\nu} \ \mathtt{ wij }
+\end{align*}
+
+
+with $h_i = h(\x_i)$.
+
+
+
+
+
+
+
+
+%================================================
+\subsection{Analytical gradients of $\psi(\x)$}
+%================================================
+
+For the \cite{ivanovaCommonEnvelopeEvolution2013} expression of the effective surfaces \Aij, we need analytical gradients in Cartesian coordinates of $\psi_i(\x_j)$.
+
+From eq. \ref{psi} we have that
+\begin{align*}
+ \psi_j(\x_i) &= \frac{ W(\x_i - \x_j, h(\x_i))}{\omega(\x_i)}\\
+\end{align*}
+
+Let $r_{ij} \equiv |\x_i - \x_j|$ and $q_{ij} \equiv \frac{r_{ij}}{h_i}$. Then
+
+\begin{align}
+ \deldx \psi_j(\x_i) &= \deldx \frac{ W(\x_i - \x_j, h(\x_i))}{\omega(\x_i)}
+ = \deldx \frac{ W(r_{ij}, h_i)}{\omega(\x_i)} \nonumber \\
+%
+ &= \frac{
+ \DELDX{W}(r_{ij}, h_i) \ \omega(\x_i) -
+ W(r_{ij}, h_i) \DELDX{\omega}(\x_i)
+ }{\omega(\x_i)^2} \nonumber \\
+%
+ &= \frac{1}{\omega(\x_i)} \DELDX{W}(r_{ij}, h_i) -
+ \frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \DELDX{}\sum_k W(r_{ik}, h_i) \nonumber \\
+%
+ &= \frac{1}{\omega(\x_i)} \DELDX{W}(r_{ij}, h_i) -
+ \frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \sum_k \DELDX{W}(r_{ik}, h_i) \label{grad_psi}
+\end{align}
+
+
+
+
+
+
+If a kernel\footnote{
+Helpful way to think of the indices: $W_j(\x_i) = W(|\x_j - \x_i|, h_i)$ is the kernel value of particle $i$ at position $\x_i$ evaluated at the position $\x_j$.
+Personally I would've used the indices the other way around for simplified thinking, but I'm going to stick to this notation because Hopkins also uses it.
+}
+is defined as
+
+\begin{align*}
+ W_j(\x_i) = W(\x_i - \x_j, h_i))
+ = \frac{1}{h_i^\nu} w \left(\frac{| \x_i - \x_j |}{h_i} \right)
+ = \frac{1}{h_i^\nu} w ( q_{ij} )
+\end{align*}
+
+
+
+
+
+and we assume that the smoothing length $h_i$ is treated as constant at this point, then the gradient of the kernel is given by
+
+\begin{align}
+ \deldx W_j (\x_i) &= \deldx \left( \frac{1}{h_i^\nu} w ( q_{ij} ) \right) \nonumber \\
+ &= \frac{1}{h_i^\nu} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\del q_{ij}(r_{ij})}{\del r_{ij}} \frac{\del r_{ij}}{\del x} \label{grad_kernel}
+\end{align}
+
+
+
+
+We now use
+\begin{align}
+ \frac{\del q_{ij}(r_{ij})}{\del r_{ij}} &= \frac{\del}{\del r_{ij}} \frac{r_{ij}}{h_i} = \frac{1}{h_i} \label{dqdr} \\
+ \DELDX{r_{ij}} &= \deldx \sqrt{(\x_i - \x_j)^2}
+ = \frac{1}{2} \frac{1}{\sqrt{(\x_i - \x_j)^2}} \cdot 2 (\x_i - \x_j) \nonumber \\
+ &= \frac{\x_i - \x_j}{r_{ij}} \label{drdx}
+\end{align}
+
+
+
+
+
+
+
+
+
+
+
+To be perfectly clear, we should in fact write
+\begin{align*}
+ r_j(\x) & \equiv | \x - \x_j | \\
+\end{align*}
+
+which again leads to
+
+\begin{align*}
+ \DELDX{r_{ij}} &= \DELDX{r_j(\x_i)} = \DELDX{r_j(\x)} \bigg{|}_{\x = \x_i} \\
+ &= \deldx \sqrt{(\x - \x_j)^2} \big{|}_{\x = \x_i}
+ = \frac{1}{2} \frac{1}{\sqrt{(\x - \x_j)^2}} \cdot 2 (\x - \x_j) \big{|}_{\x = \x_i} \\
+ &= \frac{\x - \x_j}{r_{j}(\x)} \big{|}_{\x = \x_i}
+ = \frac{\x_i - \x_j}{r_{ij}}
+\end{align*}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Inserting expressions \ref{dqdr} and \ref{drdx} in \ref{grad_kernel}, we obtain
+\begin{equation}
+ \deldx W_j (\x_i) = \frac{1}{h_i^{\nu + 1}} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\x_i - \x_j}{r_{ij}} \label{grad_kernel_final}
+\end{equation}
+
+
+
+$\frac{\del w(q_{ij})}{\del q_{ij}}$ is given by \verb|wij_dx| of \verb|kernel_deval|.
+
+
+
+
+
+
+
+Finally, inserting \ref{grad_kernel_final} in \ref{grad_psi} we get
+\begin{align}
+ \deldx \psi_j (\x_i) &=
+ \frac{1}{\omega(\x_i)} \frac{1}{h_i^{\nu + 1}} \deldr{w}(r_{ij}, h_i) \frac{\x_i - \x_j}{r_{ij}} -
+ \frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \sum_k \frac{1}{h_i^{\nu + 1}} \deldr{w}(r_{ik}, h_i) \frac{\x_i - \x_k}{r_{ik}} \nonumber \\
+ &=
+ \frac{1}{\omega(\x_i)} \frac{1}{h_i^{\nu + 1}} \mathtt{wij\_dx} \frac{\x_i - \x_j}{r_{ij}} -
+ \frac{1}{\omega(\x_i)^2} W(r_{ij}, h_i) \frac{1}{h_i^{\nu + 1}} \sum_k \mathtt{wik\_dx} \frac{\x_i - \x_k}{r_{ik}}
+\end{align}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+The definition of $r_{ij}$ requires a bit more discussion.
+Since kernels used in hydrodynamics (at least in those methods currently implemented in SWIFT) are usually taken to be spherically symmetric, we might as well have defined
+
+\begin{align*}
+r'_{ij} = |\x_j - \x_i |
+\end{align*}
+
+which would leave the evaluation of the kernels invariant [$r'_{ij} = r_{ij}$], but the gradients would have the opposite direction:
+
+\begin{align*}
+\DELDX{r'_{ij}}
+= \frac{\x_j - \x_i}{r'_{ij}}
+= - \frac{\x_i - \x_j}{r_{ij}}
+= - \DELDX{r_{ij}}
+\end{align*}
+
+So which definition should we take? \\
+
+
+Consider a one-dimensional case where we choose two particles $i$ and $j$ such that $x_j > x_i$ and $q_{ij} = | x_j - x_i | / h_i < H$, where $H$ is the compact support radius of the kernel of choice.
+Because we're considering a one-dimensional case with $x_j > x_i$, we can now perform a simple translation such that particle $i$ is at the origin, i.e. $x'_i = 0$ and $x'_j = x_j - x_i = | x_j - x_i | = r_{ij}$.
+In this scenario, the gradient in Cartesian coordinates and in spherical coordinates must be the same:
+
+
+\begin{align}
+ \frac{\del}{\del x'} W (|x_i' - x'|, h_i) \big{|}_{x' = x_j' } &= \frac{\del}{\del r_{ij}} W (r_{ij}, h_i) \nonumber\\
+%
+ \Rightarrow \quad \frac{1}{h_i^\nu} \frac{\del w(q'_{ij})}{\del q'_{ij}} \frac{\del q'_{ij}(r'_{ij})}{\del r'_{ij}} \frac{\del r'_{i}(x')}{\del x'} \big{|}_{x' = x_j' } &=
+ \frac{1}{h_i^\nu} \frac{\del w(q_{ij})}{\del q_{ij}} \frac{\del q_{ij}(r_{ij})}{\del r_{ij}} \label{spher_cart}
+\end{align}
+
+
+
+
+We have the trivial case where
+\begin{align*}
+ r'_{ij} &= |x_i' - x_j'| = | x_i - x_j | = r_{ij} \\
+ q'_{ij} &= r'_{ij}/h_i = q_{ij} \\
+ \Rightarrow \quad \frac{\del w(q'_{ij})}{\del q'_{ij}} &= \frac{\del w(q_{ij})}{\del q_{ij}}, \quad \frac{\del q'_{ij}(r'_{ij})}{\del r'_{ij}} = \frac{\del q_{ij}(r_{ij})}{\del r_{ij}}
+\end{align*}
+
+
+
+giving us the condition from \ref{spher_cart}:
+
+
+
+
+
+
+\begin{equation}
+ \frac{\del r'_{i}(x')}{\del x'} \big{|}_{x' = x_j' } = 1 = \frac{\del r_{i}(x)}{\del x}
+\end{equation}
+
+
+this is satisfied for
+\begin{align*}
+ r_j(\x) &= | \x - \x_j | \\
+ \Rightarrow r_{ij} &= | \x_i - \x_j |,\ \text{ not } \ r_{ij} = | \x_j - \x_i |
+\end{align*}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/theory/Gizmo/gizmo-implementation-details/equations.tex b/theory/Gizmo/gizmo-implementation-details/equations.tex
new file mode 100644
index 0000000000000000000000000000000000000000..66e51fccfbdc807ee9aaec4493a7700ef6f54c1c
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/equations.tex
@@ -0,0 +1,108 @@
+
+
+\section{Equations}
+
+
+
+\subsection{Partition of Unity and Related Quantities}
+
+The partition of unity is defined as:
+
+\begin{align}
+ \psi_i(\x) &= \frac{1}{\omega(\x)} W(\x - \x_i, h(\x)) \label{psi} \\
+ \omega(\x) &= \sum_j W(\x - \x_j, h(\x)) \label{omega}
+\end{align}
+
+where $h(\x)$ is some ``kernel size'' and $\omega(\x)$ is used to normalise the volume partition at any point $\x$.
+We use the compact support radius $H$ of the kernels as $h$.
+
+
+
+
+
+It can be shown that (provided $W(\x)$ is normalized, i.e. $\int_V W(\x) \D V = 1$):
+
+\begin{align}
+ V_i &= \int_V \psi_i(\x) \de V = \frac{1}{\omega(\x_i)} \\ \label{psi_volume_integral}
+ V &= \sum_i V_i \\
+ \int_V f(\x) \D V &= \sum_i f(\x_i) V_i + \mathcal{O}(h^2) \label{vol_integral_fx}
+\end{align}
+
+Eq. (\ref{vol_integral_fx}) can be derived by Taylor-expanding $f(\x)$ around $\x_i$ and integrating the first order term by parts to show that it is zero.
+
+
+
+
+
+
+
+\subsection{Meshless Hydrodynamics \`a la Hopkins}
+
+
+
+Following \cite{hopkinsGIZMONewClass2015}, we arrive at the equation
+
+\begin{equation}
+ \frac{\D}{\D t} (V_i \U_{k,i}) + \sum_j \F_{k,ij} \cdot \Aijm = 0
+\end{equation}
+
+
+with
+
+\begin{equation}
+ \Aijm^\alpha = V_i \psitilde_j^\alpha (\x_i) - V_j\psitilde_i^\alpha (\x_j) \label{Hopkins}
+\end{equation}
+
+for every component $k$ of the Euler equations and every gradient component $\alpha$
+
+The $\psitilde(\x)$ come from the $\mathcal{O}(h^2)$ accurate discrete gradient expression from \cite{lansonRenormalizedMeshfreeSchemes2008}:
+
+\begin{align}
+ \frac{\del}{\del x_{\alpha}} f(\x) \big{|}_{\x_i} &=
+ \sum_j \left( f(\x_j) - f(\x_i) \right) \psitilde_j^\alpha (\x_i) \\ \label{gradient}
+ \psitilde_j^\alpha (\x_i) &= \sum_{\beta = 1}^{\beta=\nu} \mathbf{B}_i^{\alpha \beta}
+ (\x_j - \x_i)^\beta \psi_j(\x_i) \\
+ \mathbf{B}_i &= \mathbf{E_i} ^ {-1} \\
+ \mathbf{E}_i^{\alpha \beta} &= \sum_j (\x_j - \x_i)^\alpha (\x_j - \x_i)^\beta \psi_j(\x_i)
+\end{align}
+
+
+
+where $\alpha$ and $\beta$ again represent the coordinate components for $\nu$ dimensions.
+
+
+
+
+
+
+
+
+
+
+
+
+\subsection{Meshless Hydrodynamics \`a la Ivanova}
+
+
+
+Following \cite{ivanovaCommonEnvelopeEvolution2013}, we arrive at the equation
+
+\begin{equation}
+ \frac{\D}{\D t} (V_i \U_{k,i}) + \sum_j \F_{k,ij} \cdot \Aijm = 0
+\end{equation}
+
+
+In the paper, no discrete expression for \Aij is given, instead:
+\begin{equation}
+ \Aijm = \int_V \left[ \psi_j(\x) \nabla \psi_i(\x) - \psi_i(\x) \nabla \psi_j(\x) \right] \D V
+\end{equation}
+
+
+By Taylor-expanding the gradients of $\psi$ the volume integral can be discretised as:
+
+\begin{equation}
+ \Aijm^\alpha = V_i \nabla^\alpha \psi_j (\x_i) - V_j \nabla^\alpha \psi_i (\x_j) + \mathcal{O}(h^2) \label{Ivanova}
+\end{equation}
+
+This time we get analytical gradients of $\psi$ instead of $\psitilde$.
+Note that the indices $i$ and $j$ in these equations are switched compared to the one given in \cite{ivanovaCommonEnvelopeEvolution2013} so that the evolution equations for Hopkins and Ivanova versions are identical.
diff --git a/theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.tex b/theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.tex
new file mode 100644
index 0000000000000000000000000000000000000000..237ac1ee405701eb25a2e4dadab0049fde5dd175
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/gizmo-implementation-details.tex
@@ -0,0 +1,81 @@
+\include{./header}
+
+
+
+%------------------------------------------
+%: Metadata
+%------------------------------------------
+
+\title{GIZMO Implementation Details}
+%\subtitle{Subtitle}
+\author{Bert Vandenbroucke (bert.vandenbroucke@ugent.be), Mladen Ivkovic (mladen.ivkovic@epfl.ch)}
+\date{}
+
+%------------------------------------------
+
+
+
+
+
+
+\newcommand{\Aij}{$\mathbf{A}_{ij}$}
+\newcommand{\Aijm}{\mathbf{A}_{ij}} % A_ij math
+\newcommand{\U}{\mathbf{U}}
+\newcommand{\F}{\mathbf{F}}
+\newcommand{\psitilde}{\tilde{\boldsymbol{\psi}}}
+
+
+
+
+
+\begin{document}
+
+%--------------------------------------------
+% Stuff that needs to be done before all else
+%--------------------------------------------
+%\pagestyle{plain}
+%\nocite{*} % show all entries of bibliography, even if they are not cited.
+
+
+
+%Titlepage
+\maketitle
+
+
+
+
+\input{bert}
+\input{equations}
+\input{computations}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%------------
+%Bibliography
+%------------
+\bibliography{references}
+
+
+
+
+
+
+
+\end{document}
diff --git a/theory/Gizmo/gizmo-implementation-details/header.tex b/theory/Gizmo/gizmo-implementation-details/header.tex
new file mode 100644
index 0000000000000000000000000000000000000000..a09935b761813e87ad4c469042cdf8e20fd34308
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/header.tex
@@ -0,0 +1,115 @@
+\documentclass[12pt, a4paper, english, singlespacing, parskip]{scrartcl}
+
+%\documentclass[
+%11pt, % The default document font size, options: 10pt, 11pt, 12pt
+%oneside, % Two side (alternating margins) for binding by default, uncomment to switch to one side
+%chapterinoneline, % Have the chapter title next to the number in one single line
+%english, % ngerman for German
+%singlespacing, % Single line spacing, alternatives: onehalfspacing or doublespacing
+%draft, % Uncomment to enable draft mode (no pictures, no links, overfull hboxes indicated)
+%nolistspacing, % If the document is onehalfspacing or doublespacing, uncomment this to set spacing in lists to single
+%liststotoc, % Uncomment to add the list of figures/tables/etc to the table of contents
+%toctotoc, % Uncomment to add the main table of contents to the table of contents
+%parskip, % Uncomment to add space between paragraphs
+%nohyperref, % Uncomment to not load the hyperref package
+%headsepline, % Uncomment to get a line under the header
+%]{scrartcl or scrreprt or scrbook} % The class file specifying the document structure
+
+\usepackage{lmodern} % Diese beiden packages sorgen für echte
+\usepackage[T1]{fontenc} % Umlaute.
+
+\usepackage{amssymb, amsmath, color, graphicx, float, setspace, tipa}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage[pdfpagelabels,
+ pdfstartview = FitH,
+ bookmarksopen = true,
+ bookmarksnumbered = true,
+ linkcolor = black,
+ plainpages = false,
+ hypertexnames = false,
+ citecolor = black,
+ breaklinks]{hyperref}
+\usepackage{url}
+
+\usepackage{authblk} % titlepage stuff
+\usepackage[titletoc, title]{appendix}
+
+%===================
+% BIBLIOGRAPHY
+%===================
+\usepackage[]{natbib}
+\bibliographystyle{unsrtnat}
+%DONT FORGET TO COMPILE THE BIBLIOGRAPHY WITH BIBTEX WHEN CHANGES ARE MADE.
+
+
+
+
+%--------------------------------------------
+% NEW COMMANDS
+%--------------------------------------------
+
+\newcommand{\corresponds}{\mathrel{\widehat{=}}} % equals with hat
+
+\newcommand {\arctanh}{\mathrm{arctanh}} % Atanh
+\newcommand{\arccot}{\mathrm{arccot }} % Acotanh
+
+\newcommand{\limz}[1]{\lim\limits_{#1 \rightarrow 0}} % Limes of something towars zero
+
+\newcommand{\bm}{\boldmath} % Bold font in math
+\newcommand{\dps}{\displaystyle}
+
+\newcommand{\e}{\mbox{e}} % e noncursive in math mode
+
+\newcommand{\del}{\partial} % partial diff operator
+\newcommand{\de}{\mathrm{d}} % differential d
+\newcommand{\D}{\mathrm{d}} % differential d
+\newcommand{\GRAD}{\mathrm{grad}\ } % gradient
+\newcommand{\DIV}{\mathrm{div}\ } % divergence
+\newcommand{\ROT}{\mathrm{rot}\ } % rotation
+
+\newcommand{\CONST}{\mathrm{const.\ }} % constant
+\newcommand{\var}{\mathrm{var}} % variance
+
+\newcommand{\g}{^\circ} % degrees
+\newcommand{\degr}{^\circ} % degrees
+
+\newcommand{\msol}{M_\odot} % solar mass
+
+
+\newcommand{\x}{\mathbf{x}} % x vector
+\newcommand{\xdot}{\dot{\mathbf{x}}} % x dot vector
+\newcommand{\xddot}{\ddot{\mathbf{x}}} % x doubledot vector
+\newcommand{\R}{\mathbf{r}} % r vector
+\newcommand{\rdot}{\dot{\mathbf{r}}} % r dot vector
+\newcommand{\rddot}{\ddot{\mathbf{r}}} % r doubledot vector
+\newcommand{\vel}{\mathbf{v}} % v vector
+\newcommand{\V}{\mathbf{v}} % v vector
+\newcommand{\vdot}{\dot{\mathbf{v}}} % v dot vector
+\newcommand{\vddot}{\ddot{\mathbf{v}}} % v doubledot vector
+
+\newcommand{\dete}{\mathrm{d}t} % dt
+\newcommand{\delte}{\del t} % partial t
+\newcommand{\dex}{\mathrm{d}x} % dx
+\newcommand{\delx}{\del x} % partial x
+\newcommand{\der}{\mathrm{d}r} % dr
+\newcommand{\delr}{\del r} % partial r
+
+\newcommand{\deldt}{\frac{\del}{\del t}} % shortcut partial derivative, in line
+\newcommand{\ddt}{\frac{\de}{\de t}} % shortcut total derivative, in line
+\newcommand{\DELDT}[1]{\frac{\del #1}{\del t}} % shortcut partial derivative, on fraction
+\newcommand{\DDT}[1]{\frac{\de #1}{\de t}} % shortcut total derivative, on fraction
+
+\newcommand{\deldx}{\frac{\del}{\del x}} % shortcut partial derivative, in line
+\newcommand{\ddx}{\frac{\de}{\de x}} % shortcut total derivative, in line
+\newcommand{\DELDX}[1]{\frac{\del #1}{\del x}} % shortcut partial derivative, on fraction
+\newcommand{\DDX}[1]{\frac{\de #1}{\de x}} % shortcut total derivative, on fraction
+
+\newcommand{\deldr}{\frac{\del}{\del r}} % shortcut partial derivative, in line
+\newcommand{\ddr}{\frac{\de}{\de r}} % shortcut total derivative, in line
+\newcommand{\DELDR}[1]{\frac{\del #1}{\del r}} % shortcut partial derivative, on fraction
+\newcommand{\DDR}[1]{\frac{\de #1}{\de r}} % shortcut total derivative, on fraction
+
+% replace \sum with \sum\limits
+\let\oldsum\sum
+\renewcommand{\sum}{\oldsum\limits}
diff --git a/theory/Gizmo/gizmo-implementation-details/references.bib b/theory/Gizmo/gizmo-implementation-details/references.bib
new file mode 100644
index 0000000000000000000000000000000000000000..ca2a3f57ff8f7f9e374dae7cf246a347500e9da9
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/references.bib
@@ -0,0 +1,66 @@
+
+@article{hopkinsGIZMONewClass2015,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1409.7395},
+ title = {{{GIZMO}}: {{A New Class}} of {{Accurate}}, {{Mesh}}-{{Free Hydrodynamic Simulation Methods}}},
+ volume = {450},
+ issn = {0035-8711, 1365-2966},
+ shorttitle = {{{GIZMO}}},
+ abstract = {We present two new Lagrangian methods for hydrodynamics, in a systematic comparison with moving-mesh, SPH, and stationary (non-moving) grid methods. The new methods are designed to simultaneously capture advantages of both smoothed-particle hydrodynamics (SPH) and grid-based/adaptive mesh refinement (AMR) schemes. They are based on a kernel discretization of the volume coupled to a high-order matrix gradient estimator and a Riemann solver acting over the volume 'overlap.' We implement and test a parallel, second-order version of the method with self-gravity \& cosmological integration, in the code GIZMO: this maintains exact mass, energy and momentum conservation; exhibits superior angular momentum conservation compared to all other methods we study; does not require 'artificial diffusion' terms; and allows the fluid elements to move with the flow so resolution is automatically adaptive. We consider a large suite of test problems, and find that on all problems the new methods appear competitive with moving-mesh schemes, with some advantages (particularly in angular momentum conservation), at the cost of enhanced noise. The new methods have many advantages vs. SPH: proper convergence, good capturing of fluid-mixing instabilities, dramatically reduced 'particle noise' \& numerical viscosity, more accurate sub-sonic flow evolution, \& sharp shock-capturing. Advantages vs. non-moving meshes include: automatic adaptivity, dramatically reduced advection errors \& numerical overmixing, velocity-independent errors, accurate coupling to gravity, good angular momentum conservation and elimination of 'grid alignment' effects. We can, for example, follow hundreds of orbits of gaseous disks, while AMR and SPH methods break down in a few orbits. However, fixed meshes minimize 'grid noise.' These differences are important for a range of astrophysical problems.},
+ number = {1},
+ journal = {Monthly Notices of the Royal Astronomical Society},
+ doi = {10.1093/mnras/stv195},
+ author = {Hopkins, Philip F.},
+ month = jun,
+ year = {2015},
+ keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics,Astrophysics - Instrumentation and Methods for Astrophysics,Physics - Computational Physics,Astrophysics - Astrophysics of Galaxies,Physics - Fluid Dynamics,_tablet},
+ pages = {53-110},
+ file = {/home/mivkov/Zotero/storage/23ALA3M9/Hopkins_2015_GIZMO.pdf;/home/mivkov/Zotero/storage/SHV6QDMB/1409.html}
+}
+
+@article{lansonRenormalizedMeshfreeSchemes2008,
+ title = {Renormalized {{Meshfree Schemes I}}: {{Consistency}}, {{Stability}}, and {{Hybrid Methods}} for {{Conservation Laws}}},
+ volume = {46},
+ issn = {0036-1429},
+ shorttitle = {Renormalized {{Meshfree Schemes I}}},
+ abstract = {This paper is devoted to the study of a new kind of meshfree scheme based on a new class of meshfree derivatives: the renormalized meshfree derivatives, which improve the consistency of the original weighted particle methods. The weak renormalized meshfree scheme, built from the weak formulation of general conservation laws, turns out to be \$L\^2\$ stable under some geometrical conditions on the distribution of particles and some regularity conditions of the transport field. A time discretization is then performed by analogy with finite volume methods, and the \$L\^1\$, \$L\^\textbackslash{}infty\$, and \$BV\$ stabilities of the obtained time discretized scheme are studied. From the same analogy with finite volume methods, a hybrid particle scheme is built using the Godunov method and is numerically compared to the weak renormalized scheme.},
+ number = {4},
+ journal = {SIAM J. Numer. Anal.},
+ doi = {10.1137/S0036142903427718},
+ author = {Lanson, Nathalie and Vila, Jean-Paul},
+ month = apr,
+ year = {2008},
+ keywords = {finite volume,hybrid particle schemes,meshfree methods,nonlinear conservation law,renormalization,SPH,weak scheme,_tablet},
+ pages = {1912--1934},
+ file = {/home/mivkov/Zotero/storage/QPBH9IJ6/Lanson_Vila_2008_Renormalized Meshfree Schemes I.pdf}
+}
+
+@article{ivanovaCommonEnvelopeEvolution2013,
+ archivePrefix = {arXiv},
+ eprinttype = {arxiv},
+ eprint = {1209.4302},
+ title = {Common {{Envelope Evolution}}: {{Where}} We Stand and How We Can Move Forward},
+ volume = {21},
+ issn = {0935-4956, 1432-0754},
+ shorttitle = {Common {{Envelope Evolution}}},
+ abstract = {This work aims to present our current best physical understanding of common-envelope evolution (CEE). We highlight areas of consensus and disagreement, and stress ideas which should point the way forward for progress in this important but long-standing and largely unconquered problem. Unusually for CEE-related work, we mostly try to avoid relying on results from population synthesis or observations, in order to avoid potentially being misled by previous misunderstandings. As far as possible we debate all the relevant issues starting from physics alone, all the way from the evolution of the binary system immediately before CEE begins to the processes which might occur just after the ejection of the envelope. In particular, we include extensive discussion about the energy sources and sinks operating in CEE, and hence examine the foundations of the standard energy formalism. Special attention is also given to comparing the results of hydrodynamic simulations from different groups and to discussing the potential effect of initial conditions on the differences in the outcomes. We compare current numerical techniques for the problem of CEE and also whether more appropriate tools could and should be produced (including new formulations of computational hydrodynamics, and attempts to include 3D processes within 1D codes). Finally we explore new ways to link CEE with observations. We compare previous simulations of CEE to the recent outburst from V1309 Sco, and discuss to what extent post-common-envelope binaries and nebulae can provide information, e.g. from binary eccentricities, which is not currently being fully exploited.},
+ number = {1},
+ journal = {The Astronomy and Astrophysics Review},
+ doi = {10.1007/s00159-013-0059-2},
+ author = {Ivanova, N. and Justham, S. and Chen, X. and De Marco, O. and Fryer, C. L. and Gaburov, E. and Ge, H. and Glebbeek, E. and Han, Z. and Li, X.-D. and Lu, G. and Marsh, T. and Podsiadlowski, Ph and Potter, A. and Soker, N. and Taam, R. and Tauris, T. M. and van den Heuvel, E. P. J. and Webbink, R. F.},
+ month = nov,
+ year = {2013},
+ keywords = {Astrophysics - High Energy Astrophysical Phenomena,Astrophysics - Solar and Stellar Astrophysics,_tablet},
+ file = {/home/mivkov/Zotero/storage/V6IL3NUY/Ivanova et al_2013_Common Envelope Evolution.pdf;/home/mivkov/Zotero/storage/MMFFSMV7/1209.html}
+}
+
+@phdthesis{vandenbrouckeAdvancedModelsSimulating,
+ title = {{Advanced models for simulating dwarf galaxy formation and evolution}},
+ language = {nl},
+ author = {Vandenbroucke, Bert},
+ year = {2016},
+ file = {/home/mivkov/Zotero/storage/T3UZJAWQ/Vandenbroucke - Advanced models for simulating dwarf galaxy format.pdf}
+}
+
+
diff --git a/theory/Gizmo/gizmo-implementation-details/run.sh b/theory/Gizmo/gizmo-implementation-details/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..befd7f3a73a4760703953948f7043781cbb47e07
--- /dev/null
+++ b/theory/Gizmo/gizmo-implementation-details/run.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+echo "Generating PDF..."
+pdflatex -jobname=gizmo-implementation-details gizmo-implementation-details.tex
+bibtex gizmo-implementation-details
+pdflatex -jobname=gizmo-implementation-details gizmo-implementation-details.tex
+pdflatex -jobname=gizmo-implementation-details gizmo-implementation-details.tex