A complete treatment of your question (more than you ever wanted) is given at http://homepage.tudelft.nl/q1d90/FBweb/diss.pdf, especially section 2.1.1 "Differential equation of light rays in inhomogeneous media". Trying to extract the most useful expression from that dissertation, I believe that the equation you are looking for is:
$$\nabla \Phi a = \frac{2\pi}{\lambda}n(R)$$
where
$\Phi$ = phase
$R$ = position vector
$a$ = unit vector pointing along ray
With a bit of manipulation, that turns into
$$\frac{d}{ds}\left(n\frac{dR}{ds}\right) = \nabla n$$ (equation 2.1.8 in the above reference).
The factor $ds$ can be a bit tricky since it is pointing along the ray - if you want things in X,Y coordinates then you need to worry about the length of $ds$ when it is no longer at a small angle to the X axis - it becomes $\sqrt{dx^2+dy^2}$
Snell's law in its sine-over-sine form is not the best tool to handle situations that involve incidence angles above the critical angle, and the complex angles that you're getting will only lead you astray.
Instead, the way to see things is through an explicit evanescent wave, as in this recent answer of mine. When light refracts, you have an incoming wave of the form
$$
\mathbf E_\mathrm{in}(\mathbf r,t) = \mathrm{Re}\mathopen{}\left( \mathbf E_{0,\mathrm{in}} e^{i(\mathbf k_\mathrm{in}\cdot \mathbf r- \omega t)}\right)\mathclose{},
$$
and on the other side you get a transmitted beam of the form
$$
\mathbf E_\mathrm{tr}(\mathbf r,t) = \mathrm{Re}\mathopen{}\left( \mathbf E_{0,\mathrm{tr}} e^{i(\mathbf k_\mathrm{tr}\cdot \mathbf r- \omega t)}\right)\mathclose{},
$$
where Snell's law comes in through the matching of the two wavevectors, $\mathbf k_\mathrm{in}$ and $\mathbf k_\mathrm{tr}$. These are subject to two conditions:
- The $x$ components need to match, $k_{x,\mathrm{in}} = k_{x,\mathrm{tr}} = k_x$, in order for both waves to match continuously at the boundary.
- The squares of the wavevectors need to be related to the frequency through $k^2 = n^2\omega^2$, to reflect the change in the phase velocity, and since $\omega$ is the same on both sides, they need to be related to each other through
$$k_{z,\mathrm{tr}}^2 +k_x^2 = \frac{n_2^2}{n_1^2}\left(k_x^2 +k_{z,\mathrm{in}}^2 \right).$$
The critical angle comes in when the right-hand side is too small for $k_{z,\mathrm{tr}}^2$ to be positive, i.e. when $n_1$ is so much bigger than $n_2$, and $k_x$ is so big, that the square root in
$$k_{z,\mathrm{tr}} = \pm \frac{1}{n_1}\sqrt{ n_2^2k_{z,\mathrm{in}}^2 +(n_2^2-n_1^2)k_x^2 }$$
is no longer real.
If this is the case, it's not actually a problem: you just write $k_{z,\mathrm{tr}} = i\kappa_z$, and it gives you a perfectly workable 'transmitted' field, of the form
$$
\mathbf E_\mathrm{tr}(\mathbf r,t) = \mathrm{Re}\mathopen{}\left( \mathbf E_{0,\mathrm{tr}} e^{-\kappa_{z} z}e^{i(k_x x- \omega t)}\right)\mathclose{},
$$
where what used to be an oscillatory $e^{ik_z z}$ has turned into the decaying exponential $e^{-\kappa_{z} z}$: this is the evanescent wave, which is the decaying trace of the light that's been unable to make it past the boundary, but which cannot drop discontinuously to zero.
Similarly, if you have a complex index of refraction, it means that your medium absorbs light, and there you have no alternative but to use explicit exponential decays in your description of the beams: if $n$ is complex, then the wavevector still needs to satisfy
$$k^2 = n^2\omega^2,$$
which means that it too will be complex, and you'll have a decaying field in that medium.
In this situation, I would advise you to stay well away from a complex-angles formalism, and to keep things grounded in the cartesian components of the wavevector - or to be very, very careful with those complex angles if you insist on using them. The reason for this is that the cartesian components are related to the angle and the magnitude via
\begin{align}
k_{x,\mathrm{tr}} & = k_\mathrm{tr} \sin(\theta_\mathrm{tr}) \\
k_{z,\mathrm{tr}} & = k_\mathrm{tr} \cos(\theta_\mathrm{tr}),
\end{align}
and while it is tempting to make $\theta_\mathrm{tr}$ complex and leave it at that, it's important to note that $k_\mathrm{tr}$ also needs to be complex (because $k^2 = n^2\omega^2$), but that the two imaginary parts need to cancel out when you form $k_{x,\mathrm{tr}} = k_\mathrm{tr} \sin(\theta_\mathrm{tr})$, and that's not easy to visualize at all.
That said, you do have a good point that when the second medium has a nonzero absorbance, it is no longer easy to distinguish between the 'sub-critical' incidence angles and over-critical cases with total internal reflection - and indeed there is no longer any case with full reflection of the incident power, since the transmitted wave always has a nonzero $\mathrm{Re}(k_z)$ and therefore always carries up some amount of power away from the boundary.
The way to analyze this in a way that will bring out the relevant behaviour is to plot the wavevector component $k_z$ in the complex plane as the angle of incidence is shifted from $0$ to $\pi/2$. If $n_2$ is real, then $k_z$ will be real at normal incidence, it will approach the origin as $\theta_\mathrm{in}$ increases, and once it reaches the origin it will advance back up along the imaginary axis.
If $n_2$ has a nonzero imaginary part, though, then $k_z$ will follow a similar path but it will be detached: off with a nonzero imaginary part before reaching the origin, and then up the imaginary axis but with a positive real part. The image below shows this behaviour, tracking $k_z$ for $\theta_\mathrm{in}$ between $0$ and $\pi/2$, and with $n_2$ between $2$ (dark blue) to $2+i$ (light blue).
Mathematica source through Import["http://halirutan.github.io/Mathematica-SE-Tools/decode.m"]["http://i.stack.imgur.com/f5Zz6.png"]
As you can see, if $\mathrm{Im}(n_2)$ is very small (which it is in real materials), then you can define the critical angle as the 'kink' when $k_z$ is close to the origin. However, as that imaginary part increases, the concept of critical angle simply makes less and less sense: you just have a decaying wave away from the boundary, and the incidence angle has very little effect on how fast that wave decays. It mostly stops being a useful concept to begin with, so it doesn't make a lot of sense to dig deep into the mathematics to look for a precise definition of it in that regime.
Best Answer
Let a unit vector $\;\mathbf{n}=(\rm n_1,n_2,n_3)\,, \Vert\mathbf{n}\Vert=1$. Any vector $\;\mathbf{r}\;$ could be decomposed in two components with respect to $\;\mathbf{n}\;$, see Figure-01 in the bottom \begin{equation} \mathbf{r}=\mathbf{r}_\|+\mathbf{r}_\bot \tag{01}\label{01} \end{equation} one parallel and the other normal to axis $\mathbf{n}$ respectively \begin{align} \mathbf{r}_\| &=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{r}\right)\mathbf{n} \tag{02.1}\label{02.1}\\ \mathbf{r}_\bot &=\left(\mathbf{n}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\times}\mathbf{n}= \mathbf{r}-(\mathbf{n}\boldsymbol{\cdot}\mathbf{r})\mathbf{n} \tag{02.2}\label{02.2} \end{align} that is \begin{equation} \mathbf{r}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{r}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\times} \mathbf{n} \tag{03}\label{03} \end{equation} The vectors $\;\mathbf{t},\mathbf{i}\;$ are decomposed as follows \begin{equation} \mathbf{t}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{t}\right)\boldsymbol{\times} \mathbf{n} \tag{04}\label{04} \end{equation} \begin{equation} \mathbf{i}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{i}\right)\boldsymbol{\times}\mathbf{n} \tag{05}\label{05} \end{equation} Now, Snell's Law is expressed as \begin{equation} \left(\mathbf{n}\boldsymbol{\times}\mathbf{t}\right)=\mu\left(\mathbf{n}\boldsymbol{\times}\mathbf{i}\right) \tag{06}\label{06} \end{equation} see Figure-02 in the bottom.
Equation \eqref{04} combined with \eqref{05} and \eqref{06} yields \begin{equation} \mathbf{t}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}+\mu\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right] \tag{07}\label{07} \end{equation} Taking norms in \eqref{07} and since the vector $\;\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}\;$ is normal to the vector $\;\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right]\;$ \begin{equation} \Vert\mathbf{t}\Vert^2=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)^2+\mu^2\Vert\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\Vert^2 \tag{08}\label{08} \end{equation} or \begin{equation} 1=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)^2+\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right] \tag{09}\label{09} \end{equation} so \begin{equation} \left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)=\pm\sqrt{1-\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right]} \tag{10}\label{10} \end{equation} Since the angle between $\;\mathbf{n},\mathbf{t}\;$ is less than $\;\pi/2\;$ we keep the plus sign in \eqref{10} and \eqref{07} yields finally \begin{equation} \mathbf{t}=\sqrt{1-\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right]}\mathbf{n}+\mu\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right] \tag{11}\label{11} \end{equation}