Arthur Lupton’s ray of sunshine.

As a ray passes through a material with non-uniform refractive index, its angle α relative to the normal to the surface of constant refractive index μ obeys the relation:
      μ(r) sin α(r,φ)   =   constant

We are trying to find the equation of the sunray
      r = r(φ)
     Along this ray we can express α = α(φ) as r is dependent on φ.

the refractive index is a known function of the radius from the centre of the earth, and independent of the angle φ:
      μ = μ(r)

We take:
      μ(r) sin α(φ)   =   constant
         and differentiate with respect to φ:
           μ '(r) dr
dφ 		  sin α    +    μ(r) cos α   dα
dφ 		  =   0

However, this is misleading, because the axis relevant to which α is measured does not have a fixed direction in space.
          The term   dα
dφ 		  should be replaced by   dα
dφ 		+   1
We can see the correctness of this by considering the case where the refractive index is constant so that
           μ '(r)   =   0
The direction of the ray is constant and any increase in φ results in a corresponding decrease in α.

Our equation should be:
           μ '(r) dr
dφ 		  sin α    +    μ(r) cos α   dα
dφ 		   +    μ(r) cos α   =   0
Dividing through by sin α (which is only zero for a straight up ray in which we are not interested):
           μ '(r) dr
dφ 		   +    μ(r) cot α   dα
dφ 		   +    μ(r) cot α   =   0

The above equation is obtained from the consideration of the refraction. The goemetry of the situation relates α to r and φ by
           cot α   =     1
r 		  dr
dφ

We substitue for dr/dφ  in the previous equation:
           r μ '(r)   cot α    +    μ(r) cot α   dα
dφ 		   +    μ(r) cot α   =   0
           r μ '(r)    +    μ(r) dα
dφ 		   +    μ(r)   =   0

           –   dα
dφ 		  =   1   +   r μ '(r)
μ(r)
which can be integrated using Runge-Kutta when combined with:
           dr
dφ 		  =   r cot α

If we substitute θ + φ   =   ½ π – α, we get:
           dθ
dφ 		  +   1   =   1   +   r μ '(r)
μ(r)
           dθ
dφ 		  =   r μ '(r)
μ(r)
        and
           dr
dφ 		  =   r tan ( θ + φ )
which is Arthur’s pair of equations but for a minus sign.

I've not been able to understand Arthur’s figure 4 well enough to see his reasoning. I reckon that dα/dφ  is certainly –1 when the refractive index is constant, and α must tend to zero as one approaches the sun. In short, I reckon that dα/dφ  is always negative. By the same argument d(θ+φ)/dφ  must be positive, but if I understand Arthur’s angles correctly θ  is the angle between the ray and a fixed direction which is tangential to the earth’s surface at the observer, and I believe that this decreases as one approaches the sun. Note that μ '(r) is negative.