# Ray tracing through a plane where the refractive index is a function of the distance from the origin

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?

• Υour link The Refractive Invariant refers to this one Deriving the Differential Equation for Refraction. – Frobenius Jan 14 '17 at 6:09
• Yes, but I'm not entirely sure how finding the refraction differential in terms of $\theta$ (and $dn$) leads to a straightforward way of finding the path. I'd prefer some sort of equation in the plane in terms of $y$ and $z$, and was wondering if there was a different way of deriving the invariance relation to achieve this. – user38762 Jan 14 '17 at 6:21
• Perhaps there is a simpler formula for the zenith distance? – user38762 Jan 14 '17 at 7:47

Inspired by the solution given to the brachistochrone problem by the calculus of variations I found a differential equation but in cylindrical plane coordinates $\:(r,\phi)\:$. Note that the brachistochrone problem is identical to find the light path in a medium with variable refractive index of the form $\:n(y)=\textrm{constant}\cdot y^{-1/2}\:$ where $\:y\:$ a rectilinear coordinate, for example on the vertical line.

In our case the time to travel from point 1 to point 2 is $$t_{12}=\int\limits_{1}^{2}\mathrm{d}t=\int\limits_{1}^{2} \dfrac{\mathrm{d}s}{v} \tag{01}$$ where $\:s\:$ the arc-length parameter $$\mathrm{d}s=\sqrt{\mathrm{d}r^{2}+\left(r\mathrm{d}\phi\right)^{2}}=\sqrt{r^{2}+r'^{\,2}}\:\mathrm{d}\phi, \quad r'\equiv \dfrac{\mathrm{d}r}{\mathrm{d}\phi} \tag{02}$$ and $\:v\:$ the speed of light, a function of the radius $\:r\:$ $$v(r)=\dfrac{c_{0}}{n(r)}, \quad c_{0}=\textrm{speed of light in empty space} \tag{03}$$ where $\:n(r)\:$ the variable refractive index.

Least time means $$c_{0}\cdot t_{12}=\int\limits_{1}^{2}L\left(r,r',\phi\right)\mathrm{d}\phi=\text{extremum} \tag{04}$$ where $$L\left(r,r',\phi\right) \equiv n(r)\sqrt{r^{2}+r'^{\,2}} \tag{05}$$ is the Lagrangian of the problem. Note that the Lagrangian does not depend explicitly on the independent variable $\:\phi\:$ so by the Beltrami Identity (1) the Euler-Lagrange equation is equivalent to this one $$r'\dfrac{\partial L}{\partial r'}-L = A = \text{constant of the motion} \tag{06}$$ Inserting (05) in (06) we have the differential equation in cylindrical plane coordinates $$n(r)\dfrac{r^{2}}{\sqrt{r^{2}+r'^{\,2}}}= A = \text{constant}=\textrm{The Refractive Invariant} \tag{07}$$ I believe that comparing this with the Refractive Invariant of the link given we must have (2) $$\sin \mathsf{z} = \dfrac{r}{\sqrt{r^{2}+r'^{\,2}}}= \dfrac{r}{\sqrt{r^{2}+\left(\dfrac{\mathrm{d}r}{\mathrm{d}\phi}\right)^{2}}} \tag{08}$$ The difficulty (or even the possibility) to have analytical solutions to the differential equation (07) depends upon the function $\:n(r)\:$.

Note that for $\:n(r)=\textrm{constant}\:$, the differential equation (07) becomes $$\dfrac{r^{2}}{\sqrt{r^{2}+r'^{\,2}}}= A' = \textrm{constant} \tag{09}$$ with a general solution $$r=\dfrac{A'}{\cos\left(\phi+\phi_{0} \right)} \tag{10}$$ that is the equation of a straight line as expected (no refraction).

(1)See here: Beltrami identity

(2) equation (08) is written as $$\sin \mathsf{z}=\dfrac{1}{\sqrt {1+\left(\dfrac{\mathrm{d}r}{r\mathrm{d}\phi}\right)^{2}}} \tag{08a}$$ that is $$\cot \mathsf{z}=\dfrac{\mathrm{d}r}{r\mathrm{d}\phi} \tag{08b}$$ identical to the 3rd equation in Deriving the Differential Equation for Refraction with $\:\phi \equiv \theta, r \equiv R\:$. We can express the right hand side of equation (08b) in terms of the $\:\left(x,y\right)\:$ cartesian coordinates from equations \begin{align} r & =\sqrt{x^{2}+y^{2}} \tag{09a}\\ \tan \phi & =\dfrac{y}{x} \tag{09b} \end{align} From equations (09) \begin{align} \mathrm{d}r & =\dfrac{x\mathrm{d}x+y\mathrm{d}y}{\sqrt{x^{2}+y^{2}}} \tag{10a}\\ \mathrm{d}\phi & =\dfrac{-y\mathrm{d}x+x\mathrm{d}y}{x^{2}+y^{2}} \tag{10b} \end{align} Inserting (09a),(10a),10b) in (08b) we have $$\cot \mathsf{z} =\dfrac{1+\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)\biggl(\dfrac{y}{x}\biggr)}{\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right) - \biggl(\dfrac{y}{x}\biggr)}=\cot \left[\tan^{-1}\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)-\tan^{-1}\biggl(\dfrac{y}{x}\biggr) \right] \tag{11}$$ or $$\mathsf{z} =\tan^{-1}\left(\dfrac{\mathrm{d}y}{\mathrm{d}x}\right)-\tan^{-1}\biggl(\dfrac{y}{x}\biggr) \tag{12}$$ Replacing $\:\mathsf{z}\rightarrow\theta'\:,\: x\rightarrow y\:,\:y\rightarrow z\:$ we have the OP's equation but with wrong sign $$\theta' =\tan^{-1}\left(\dfrac{\mathrm{d}z}{\mathrm{d}y}\right)-\tan^{-1}\biggl(\dfrac{z}{y}\biggr)=-\theta \tag{12}$$

(3) Two Examples : In the first example, with $\:n(r)=\mathrm{constant}\cdot r^{-2}\:$, the solution of the differential equation (07) is $\: r(\phi)=c_{1}\,\sin(\phi+c_{2})\:$ where $\:c_{1},c_{2}\:$ are constants. This is the parametric equation of a circle. So the light path is a circular arc. In the second example, with $\:n(r)=\mathrm{constant}\cdot r^{-1}\:$, the solution of the differential equation (07) is $\: r(\phi)=c_{1}\,e^{c_{2}\,\phi} \:$ where $\:c_{1},c_{2}\:$ are positive constants. This is the parametric equation of a spiral. So the light path is a spiral arc.