I asked a similar question on the Physics site but not get the answer I was looking for. Every derivation I have read skips over the details when it comes to deriving Snell's law and $\omega_i = \omega_r = \omega_t$ for the electric fields on both sides of the boundary. So, I have attempted to prove it (just the latter one since it is more simple and the former proof should be similar).
Let the boundary be at $z=0$. Starting from the boundary conditions we have
$ \mathbf{E}_{Ti} \exp i( \omega_i t - \mathbf{k}_i \cdot \mathbf{r} ) + \mathbf{E}_{Tr} \exp i( \omega_r t - \mathbf{k}_r \cdot \mathbf{r} ) = \mathbf{E}_{Tt} \exp i( \omega_t t - \mathbf{k}_t \cdot \mathbf{r} ) $
where $T$ stands for the transverse components. Since this needs to hold for all $x$ and $y$ then it needs to hold for $x=y=0$. So we have
$ \mathbf{E}_{Ti} \exp i( \omega_i t ) + \mathbf{E}_{Tr} \exp i( \omega_r t ) = \mathbf{E}_{Tt} \exp i( \omega_t t ) $
Now, lets assume that at least one of these is different from the others. WLOG let $\omega_i \neq \omega_r = \omega_t$. Then we can write
$ \frac{\mathbf{E}_{Ti}}{\mathbf{E}_{Tt} - \mathbf{E}_{Tr}} = \exp i(\omega_r - \omega_i)$t.
But since the left side is constant, the only way for the right side to be constant for all $t$ is for $\omega_r = \omega_i $. Therefore the assumption that at least one frequency is different is wrong and so all of them are equal.
Is this proof correct? Is there a way to prove it directly?
I think your proof is in the right direction, but you don't explain why $\omega_r=\omega_t$ can be assumed WLOG. The fact that $\vec{E}$ is a vector also means you can't divide by it.
Another way of looking at it
Another proof would be to proceed by showing linear dependence. Every value of $t$ gives you a linear equation. I can't count the number of linear degrees of freedom here off the top of my head (Maxwell's equations and boundary conditions put other conditions on $E_{i,r,t}$), but once you plug in more values of $t$ than there are degrees of freedom, you'll have an overdetermined problem. Either some of the equations are linearly dependent, or there is no nonzero solution. The linear independence of the equations (for all $t$) gives you $\omega_i=\omega_r=\omega_t$.
Yet another way of looking at it
This is essentially just the Fourier transform in time. A case simpler than Maxwell's equations but which has all the same features is the simple wave equation. So, first solve $\partial_t^2 \psi-c(\vec{r})^{2}\nabla^2\psi=0$, where $c(\vec{r})$ is the speed of light at position $r$. You can immediately Fourier transform in time to get $(-i\omega)^2\phi-c(\vec{r})^{-2}\nabla^2\phi=0$. That is to say: we can always find a representation of $\psi$ as a sum (integral) of different single frequency components, each of which also satisfies the wave equation.
We can go a bit further too. If $c(x,y,z)$ has y and z symmetry (so, maybe x<0 is shallow water and x>0 is deep water), then you can do the same thing with the wavenumbers $k_y$ and $k_z$.
A more physicsy way to say this is that "the frequency, the $k_x$ component, and the $k_z$ component of an incoming plane wave don't change when reflecting on a boundary".