The index of refraction, $n$, for some medium, is the speed of light in a vacuum divided by the speed of light in the medium ($n = \frac{c}{v}$).
Consider some light source, $p_s$ (at coordinates $[x_p,y_p,z_p]$), and some light receiver, $p_r$ (at coordinates $[x_r,y_r,z_r]$). Neglecting the effects of refraction, we might calculate the time of flight for a light signal, moving in a straight line from $p_s$ to $p_r$, as:
$T = \frac{1}{c} \cdot \sqrt{(x_p - x_r)^2 + (y_p - y_r)^2 + (z_p - z_r)^2 }$
For a constant index of refraction, we simply replace $c$ by $\frac{c}{n}$. How, however, might we calculate the time of flight for a light signal traveling in a "straight line" when the index of refraction varies as a function of some coordinate (say, depth ($z$))? In this case, of course, light doesn't actually travel in a straight line, in the common sense.
What I have so far...
Note: This could be very wrong. It's just where I've gotten so far.
Using the principle of least time (see), we have that:
$T = \int dt$
We have that $v = \frac{ds}{dt}$. Thus, $T = \int \frac{ds}{v}$
Suppose $n$ is a function of $z$. Then, $n = n(z)$. Using $v = \frac{c}{n(z)}$, we have that:
$T = \int \frac{n(z)}{c} ds$
Now, by Pythagoras $ds^s = dx^2 + dy^2 + dz^2$. So, $ds = \sqrt{dx + dy +dz}$. Then,
$T = \int \frac{n(z)}{c} \sqrt{dx + dy + dz}$
What I want, in the end, is to have a formula such that I can plug in $[x_p,y_p,z_p]$, $[x_r,y_r,z_r]$, $n(z)$, $c$, and get $T$.
EDIT
I've found a reference, here. Unfortunately, I don't quite follow everything. I recreate the argument here. I use the above example, where $n$ is a function of depth, $z$. We simplify this to a problem in 2-dimensions.
We know that the index of refraction, $n$, at a point $z,x$ is equal to $c/v$, where $v$ is the velocity of light at that point. We parameterize the path by equations $x = x(u)$ and $z = z(u)$. The "optical path length" from a point $A$ to a point $B$ is then given by:
$$L = \int_A^B n\;ds = \int_A^B n\;\sqrt{z^2 + x^2}\;du$$
I don't quite follow what is being done at this point. I know that we want to minimize $L$... I'm not sure, however, how this step achieves this.
To make this integral an extremum, let $f$ denote the integrand function:
$f(z,x,\dot{z},\dot{x}) = n(z,x)\sqrt{z^2+x^2}$
Then the Euler equations are:
$$\frac{\partial n}{\partial z} = \frac{d}{du}\Bigg(\frac{\partial f}{\partial \dot{z}}\Bigg)$$
$$\frac{\partial n}{\partial x} = \frac{d}{du}\Bigg(\frac{\partial f}{\partial \dot{x}}\Bigg)$$
Which gives:
$\frac{\partial n}{\partial z} \sqrt{\dot{z}^2 + \dot{x}^2} = \frac{d}{du} \Bigg[\ \frac{n\dot{z}}{\sqrt{\dot{z}^2 + \dot{x}^2}}\Bigg]$
$\frac{\partial n}{\partial x} \sqrt{\dot{z}^2 + \dot{x}^2} = \frac{d}{du} \Bigg[\ \frac{n\dot{x}}{\sqrt{\dot{z}^2 + \dot{x}^2}}\Bigg]$
If we define parameter $u$ as spatial path length $s$, then $z^2 + x^2 = 1$, and the above equations reduce to:
$$\frac{\partial n}{\partial z} = \frac{d}{ds}\Bigg(n \frac{dz}{ds}\Bigg)$$ $$\frac{\partial n}{\partial x} = \frac{d}{ds}\Bigg(n \frac{dx}{ds}\Bigg)$$
I could be looking at this wrong, however this is my interperetation...
So, since (in my example) $n$ is a function of $z$, it occurs to me that we can drop the partial derivatives, and have:
$$\frac{dn}{dz} = \frac{d}{ds}\Bigg(n \frac{dz}{ds}\Bigg)$$
Then:
$$L = \frac{d}{ds} = \frac{\frac{dn}{dz}}{\Bigg(n \frac{dz}{ds}\Bigg)}$$
Suppose $n(z) = e^z$, for example...
$$L = \frac{e^z}{\Bigg(e^z \frac{dz}{ds}\Bigg)}$$
$\frac{dz}{ds}$ should be a differential equation, yes? Perhaps I'm being slow about this, but it isn't immediately clear to me how I'd then use starting and ending values for $x$,$y$, and $z$ to obtain an $L$. As I recall, the derivative can't simply be "flipped" (i.e. this isn't $= \frac{ds}{dz}$)... but maybe I'm wrong.
Short answer:
There is no formula for arbitrary $n(z)$. However, for many specific functional forms of $n(z)$ (e.g. $n(z) = e^z$) you can derive a formula for $T$ as a function of $[x_p,y_p,z_p]$ and $[x_r,y_r,z_r]$ (and $c$ of course). To derive such a formula you need to solve a differential equation involving $n(z)$.
Background:
Finding paths of least time is one of the classic problems of a subject called calculus of variations. Basically what we try to do in the calculus of variations is find a function which minimizes some integral. This is similar to what we do in ordinary calculus when we try to find the maximum or minimum value of some function, but slightly more complicated. (The idea is either case is to look for where the derivative is zero.) You've already written down the integral that we want to minimize in this case, namely $$T=\int_{\text{path}} \frac{n(z)}{c} ds$$ Once we know which path minimizes this integral, we can just plug it in and evaluate the integral to find the minimum time.
How to solve it:
There is a set of differential equations for the path of least time, which you've also written down: $$\begin{align}\frac{\partial n}{\partial z} = \frac{d}{ds}\left(n \frac{dz}{ds}\right) \\ \frac{\partial n}{\partial x} = \frac{d}{ds}\left(n \frac{dx}{ds}\right) \end{align}$$ Here's how to interpret these formulas: $x=x(s)$ and $z=z(s)$ are coordinates of the path we want to find, parametrized by the path length $s$. Then $n=n(z)=n(z(s))$ is also implicitly a function of $s$. If we know the form of $n(z)$ then we can try to solve these differential equations. E.g. let's take $n(z) = e^z$. Then the two equations above become $$\begin{align} e^z = \frac{d}{ds}\left(e^z \frac{dz}{ds}\right) \\ 0 = \frac{d}{ds}\left(e^z \frac{dx}{ds}\right) \end{align}$$ The first of these equations can be written $e^z = \frac{d^2}{ds^2} e^z$, which implies $e^z = c_1 e^s + c_2 e^{-s}$, or $z(s) = \log\left(c_1 e^s + c_2 e^{-s}\right)$ (for some integration constants $c_1, c_2$). Then the second equation above can be integrated once to give $c_3 = e^z \frac{dx}{ds}$, which means $\frac{dx}{ds} = c_3 e^{-z} = c_3/\left(c_1 e^s + c_2 e^{-s}\right)$. This can be integrated to give a formula for $x(s)$. In the end we find $$\begin{align} z(s) = \log\left(c_1 e^s + c_2 e^{-s}\right) \\ x(s) = \tan^{-1}\left(\sqrt{\frac{c_1}{c_2}} e^s\right)\frac{c_3}{\sqrt{c_1 c_2}} + c_4 \end{align}$$ You can see how this is a bit messy already, but doable. It will get messier as we proceed to the next steps, and in the interest of clarity I will only say what you need to do, rather than writing out all the math.
The next step is to invert these formulas for $x(s)$ and $z(s)$ to find the constants of integration $c_i$ in terms of the starting and ending points $[x_p,z_p]$ and $[x_r,z_r]$. Two of your equations for this step will be $x(0)=x_p$, $z(0)=z_p$. There is also a consistency condition $x'(0)^2 + z'(0)^2 =1$ which comes from the fact that we parametrized by arc length $s$. Then there will be the equations for the final point $x(s_f) = x_r$ and $z(s_f) = z_r$ which must be solved for both the remaining unknown integration constants and the final path length $s_f$.
Once you have found all of the $c_i$'s and the final path distance $s_f$, you can plug all of this into the formula for $T$, $T=\int \frac{n(z)}{c} ds$, to get your final formula.
Summary:
As you can see, the hard part of this comes from solving the differential equation for the path of least time. It can get algebraically heavy, and for some forms of $n(z)$ it will have no analytical solution. That's why there will be no general formula for time $T$. However, for many forms of $n(z)$ it is solvable, and in general you can solve the differential equations numerically if you really need that information.
If you want to learn more about this, I would recommend looking into any book on classical mechanics. Lagrangian mechanics, in particular, proceeds very much along the same lines.