Differential equation for path in radial field Suppose that in a plane parallel to the yz-plane, the index of refraction $n$ is a function of the distance from the origin, $R$, i.e., $n = n(R)$.
We know that (e.g., http://aty.sdsu.edu/explain/atmos_refr/invariant.html), that $n(R)R\sin\theta = \mathrm{constant}$ at every point, where $\theta$ is analogous to the "$z$" angles shown in the diagram:
With the relation above I have the information to calculate $\theta$ at every point, but I'd like to recast it into a differential equation that gives me the path of a ray in terms of $y$ and $z$. Is it possible to find a form for $\frac{dz}{dy}$? I could write, for instance, $\theta = \tan^{-1}\left(\frac{z}{y}\right) - \tan^{-1}\left(\frac{dz}{dy}\right)$. Or would I have to parametrize $y$ and $z$? How would I do this? What is the easiest way of numerically finding the path?

I think it is better to work in polar coordinates and then convert back to Cartesian if necessary.
Specifically, consider the $\{\log w\}$ plane (or rather, strip) for $w = y + iz$. The circles will map to horizontal lines $\mathop{\mathrm{Im}} \log w = \mathrm{const}$, the ray you need to find will map to... well, something, but because the logarithm is conformal, the angles between the two at each point will stay the same. Therefore, for $\log(y + iz) = \log R + i\phi$, you will get (schematically)
$$n(R)R\sin(\arctan d\phi/d\log R) = C,$$
$$d\phi/d\log R = \tan\arcsin(C / n(R)R) = \sqrt{[n(R)R/C]^2 - 1}\,,$$
$$d\phi/dR = \sqrt{[n(R)/C]^2 - R^2}\,,$$
which looks a lot more tractable to me.
