I'm trying to use the formalism from 'Generalized beam matrices. III. Application to diffraction analysis' and 'Generalized beam matrices: Gaussian beam propagation in misaligned complex optical systems' by Anthony A. Tovar and Lee W. Casperson to calculate how misaligned optical elements reshape Gaussian beam parameters.
The formalism uses a 3x3 matrix approach instead of the conventional 2x2 matrix. Additionally to the complex beam parameter q, there is a complex displacement parameter S.
$$\frac{1}{q_x}=\frac{1}{R_x}-i\frac{2}{\beta_0w_x^2}$$ $$S_x=\beta_0\left(-\frac{d_{xa}}{q_x}+d_{xa}'\right)$$
with radius of curvature R, beam size w, and $\beta_0=\frac{2\pi n_0}{\lambda}$ with refractive index $n_0$ and wavelength $\lambda$. The index x denotes that in this case it's in x-direction only, because beams and elements are not necesarily spherically symmetric. The displacement parameter S contains the parameters $d_{xa}$ and $d_{xa}'$. The paper describes this as "position" and "slope of the beam". So my interpretation was that $d_{xa}$ is a displacement function with
$$d(z)=(z-z_0)\tan(\epsilon)+d_0$$ that describes that the beam's optical axis is not necesarily aligned with the "lab axis" z. $z_0$ is the position of the beam waist, $d_0$ the offset at the waist position, and $\epsilon$ the angle between the z-axis and the optical axis. Performing a derivation leads to $$d'=\frac{dd(z)}{dz}=\tan(\epsilon)$$ which would be the slope of the beam. As an optical axis of the beam I chose: $$z_\text{beam}=\frac{z-z_0}{\cos(\epsilon)}$$ and then with this there's of course beam size $w(z_\text{beam})$ and radius of curvature $R(z_\text{beam})$: $$w(z_\text{beam})=w_0\sqrt{1+\left(\frac{z_\text{beam}}{z_R}\right)^2}=w_0\sqrt{1+\left(\frac{z-z_0}{z_R\cos(\epsilon)}\right)^2}$$
$$R(z_\text{beam})=z_\text{beam}\cdot\left(1+\left(\frac{z_R}{z_\text{beam}}\right)^2\right)=\frac{z-z_0}{\cos(\epsilon)}\left(1+\left(\frac{z_R\cos(\epsilon)}{z-z_0}\right)^2\right)$$
with the Rayleigh length $z_R=\frac{\pi w_0^2n}{\lambda}$. As a simple case I use the matrix for a tilted linear surface from the paper:
$$ \begin{bmatrix} A_x & B_x & 0 \\ C_x & D_x & 0 \\ G_x & H_x & 1 \\ \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \frac{k_{01}}{k_{02}} & 0 \\ (k_{01}-k_{02})\tan\theta & 0 & 1 \\ \end{bmatrix} $$ with $k_0\approx\beta_0$. $\theta$ is the tilt of the surface. q and S then transform like this: $$\frac{1}{q_2}=\frac{C+\frac{D}{q_1}}{A+\frac{B}{q_1}}$$ $$S_2=\frac{S_1+G+\frac{H}{q1}}{A+\frac{B}{q_1}}$$
Both q and S are complex, so splitting both into real and imaginary parts allows to extract the values $w_2$, $R_2$ from $q_2$, and $d_2$ and $d'_2$ from $S_2$ of the transformed beam. I'll skip the extraction of $w_2$ and $R_2$ because that works as intended. Even $d_2$ works. The only thing that does not work for some reason is the $d'$ transformation.
If I seperate $S_2$ inte real and imaginary parts, then the solution I get is $$\Re(S_2)=\beta_{0,2}\left(-d_2\Re\left(\frac{1}{q_2}\right)+d'\right)$$ $$\Im(S_2)=-\beta_{0,2}d_2\Im\left(\frac{1}{q_2}\right)$$
which then yields
$$d_2=\frac{\Im(S_2)}{\beta_{0,2}\Im\left(\frac{1}{q_2}\right)}$$ $$d'_2=\frac{1}{\beta_{0,2}}\left(\Re(S_2)-\Im(S_2)\frac{\Re\left(\frac{1}{q_2}\right)}{\Im\left(\frac{1}{q_2}\right)}\right)$$
and then since $d'_2=\tan(\epsilon_2)$ I can obtain the new tilt as $\epsilon_2=\arctan(d'_2)$. But it does not work. If I use a beam let's say with these parameters: $w_0=258\cdot10^{-6}\,$m, $z_0=0.646\,$m, $\lambda=1550\cdot10^{-9}\,$m, without any misaligment, meaning $d_0=0$, $\epsilon=0$ with a boundary at $z=0.273\,$m that has an angle of $\theta=\frac{\pi}{2}$ with $n_1=1.0003$ (air) and $n_2=1.444$, I get an angle of $\epsilon_2=-17.08^\circ$. According to Snells law however, this angle should be $15.67^\circ$.
At this point I've checked every single step many times and I absolutely cannot find any errors. The only possible issue that's really left in my mind is that I might not correctly interpret what $d$ and $d'$ are. Edit: To maybe make this more digestible: If the initial beam has no misalignment, then $S_1=0$ and subsequently $S_2=G=(k_{01}-k_{02})\tan\theta$, and $d'_2=\frac{(k_{01}-k_{02})\tan\theta}{\beta_{0,2}}$.
