Both results are correct. As Luboš Motl pointed out in his comment, we can get the geometrical optics approach answer from the wave method answer by averaging it over 1 full period.
Perhaps you made a mistake somewhere when averaging.
If we calculate the average carefully:
$\bar{T}=\frac{1}{2\pi}\int _{-\pi}^{\pi}T d\varphi$
$\bar{T}=\frac{1}{2\pi} \times T_{01}T_{12} \int_{-\pi}^{\pi}\frac{1}{(r_{10}r_{12})^2-2r_{10}r_{12}cos\varphi+1}d\varphi$
we'll get this
$\bar{T}=\frac{1}{2\pi} \times T_{01}T_{12}\left|\frac{2tan^{-1}\left(\frac{1+r_{10}r_{12}}{1-r_{10}r_{12}}tan{\frac{\varphi}{2}}\right)}{1-(r_{10}r_{12})^2}\right|_{-\pi}^{\pi} $
you can check here
now it's not hard to see that
$\bar{T}=\frac{1}{2\pi} \times T_{01}T_{12}\frac{(2\pi)}{1-R_{10}R_{12}}=\frac{T_{01}T_{12}}{1-R_{10}R_{12}}$
same as the one obtained from geometric optics method
EDIT(to respond to the EDIT part of the question):
First of all, let's make everything clear so that the equation that you quoted from Wikipedia make sense. The amplitude here means the magnitude of electric field. The $t$ simply means the ratio of transmitted and incoming electric field (We were using different definition of $t$ in the calculations above, the $t$ that we use above includes other factors beside the ratio of electric fields). And the transmittance $T$ here represents the fraction of power transmitted to the medium 2. Here the power $P$ is proportional to
$P\propto IA$
where the intensity is proportional to $I\propto nE^2$. And the beam area is proportional to $A\propto cos\theta$, it's because the beam cross section gets smaller as it bends toward the boundary plane (see this image). Putting them together we have
$P\propto nE^2cos\theta$
so from the new definition of transmittance $T$ we get
$T=\frac{n_2 cos\theta_t}{n_1 cos\theta_i}t^2$
and since the incoming and reflected rays have the same $cos\theta$ and $n$, the reflectance is simple
$R=r^2$
we have calculated that
$\bar{t^2}=\frac{t_{01}^2t_{12}^2 }{1-R_{10}R_{12}}$
multiply both sides with $\frac{n_2 cos\theta_t}{n_0 cos\theta_i}$
$\bar{T}=\frac{t_{01}^2t_{12}^2 }{1-R_{10}R_{12}}\frac{n_2 cos\theta_t}{n_0 cos\theta_i}$
We can check
$T_{01}T_{12}=\left(\frac{n_1 cos\theta_m}{n_0 cos\theta_i}t_{01}^2\right)\left(\frac{n_2 cos\theta_t}{n_1 cos\theta_m}t_{12}^2\right)=\frac{n_2 cos\theta_t}{n_0cos\theta_i}t_{01}^2t_{12}^2$
Thus again we have
$\bar{T}=\frac{T_{01}T_{12}}{1-R_{10}R_{12}}$
Suppose we define the following $\zeta = \ln{\varepsilon}$ and $\xi = \ln{\mu}$, where $\varepsilon$ and $\mu$ are the permittivity and permeability, respectively. In a system with no sources (i.e., $\mathbf{j} = 0$ and $\rho_{c} = 0$), then we know that $\nabla \cdot \mathbf{D} = 0$, where $\mathbf{D} = \varepsilon \ \mathbf{E}$ and $\mathbf{B} = \mu \ \mathbf{H}$. After a little vector calculus we can show that:
$$
\nabla \cdot \mathbf{E} = - \nabla \zeta \cdot \mathbf{E} \tag{0}
$$
Using this and some manipulation of Faraday's law and Ampêre's law, we can show that the general differential equation in terms of electric fields only is given by:
$$
\left( \mu \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \nabla^{2} \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \nabla \zeta + \left( \nabla \zeta \cdot \nabla \right) \mathbf{E} + \nabla \left( \zeta + \xi \right) \times \left( \nabla \times \mathbf{E} \right) \tag{1}
$$
We can get a tiny amount of reprieve from this by assuming that the permeability is that of free space, i.e., $\nabla \xi = 0$. If we further argue that the only direction in which gradients matter is along $\hat{x}$ and that the incident wave vector, $\mathbf{k}$, is parallel to this, then we can further simplify Equation 1 to:
$$
\left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \zeta' \hat{x} + \left( \zeta' \frac{ \partial }{ \partial x } \right) \mathbf{E} + \zeta' \hat{x} \times \left( \nabla \times \mathbf{E} \right) \tag{2}
$$
where $\zeta' = \tfrac{ \partial \zeta }{ \partial x }$.
After some more manipulation, we can break this up into components to show that:
$$
\begin{align}
\text{x : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{x} & = E_{x} \zeta'' + \zeta' \frac{ \partial E_{x} }{ \partial x } \tag{2a} \\
\text{y : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{y} & = 0 \tag{2b} \\
\text{z : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{z} & = 0 \tag{2c}
\end{align}
$$
In the limit where the incident wave is entirely transverse, then $E_{x} = 0$ and the x-component (Equation 2a) is entirely zero.
Next you assume that $\mathbf{E} = \mathbf{E}_{o}\left( x \right) e^{i \omega t}$, where $\omega$ is the frequency of the incident wave. Then there will be incident, reflected, and transmitted contributions to the total field at any given point (well the transmitted is always zero in the first medium, of course). Any incident and transmitted contributions with have $\mathbf{k} \cdot \mathbf{x} > 0$ while reflected waves will satisfy $\mathbf{k} \cdot \mathbf{x} < 0$. You define the ratio of the reflected-to-incident fields (well impedances would be more appropriate) to get the coefficient of reflection.
Simpler Approach
A much simpler approach is to know where to look for the answer to these types of questions. As I mentioned in the comments, there has been a ton of work on this very topic (i.e., spatially dependent index of refraction) done for the ionosphere. If we look at, for instance, Roettger [1980] we find a nice, convenient equation for the reflection coefficient, $R$, as a function of the index of refraction, given by:
$$
R = \int \ dx \ \frac{ 1 }{ 2 \ n\left( x \right) } \frac{ \partial n\left( x \right) }{ \partial x } \ e^{-i \ k \ x} \tag{3}
$$
There is no analytical expression for $R$ for your specific index of refraction. However, numerical integration is not difficult if one knows the values for $d$ and $\epsilon$. Note that if we do a Taylor expansion for small $\epsilon$, then the integrand (not including the exponential) is proportional to cosine, to first order in $\epsilon$ (cosine times sine if we go to second order).
References
- Gossard, E.E. "Refractive index variance and its height distribution in different air masses," Radio Sci. 12(1), pp. 89-105, doi:10.1029/RS012i001p00089, 1977.
- Roettger, J. "Reflection and scattering of VHF radar signals from atmospheric refractivity structures," Radio Sci. 15(2), pp. 259-276, doi:10.1029/RS015i002p00259, 1980.
- Roettger, J. and C.H. Liu "Partial reflection and scattering of VHF radar signals from the clear atmosphere," Geophys. Res. Lett. 5(5), pp. 357-360, doi:10.1029/GL005i005p00357, 1978.
Best Answer
The phase-shift occurs due to the complex nature of the reflection coefficient. The so-called critical angle is given by:
$sin(\theta_M) = \frac{n_2}{n_1}$
As long as $\theta<\theta_M$ we have only partial reflection and a real valued reflection coefficient $R$.
As soon as the critical angle is exceeded $(\theta>\theta_M)$, we have $\mid R \mid=1$ and total reflection of the light occurs. $R$ is now complex and a phase shift is imposed on the reflected light, we write:
$R=e^{2j\phi}$
wehre $\phi$ is the phase shift between the incident and reflected wave.
Now, we know what the value of $R$ is, both for TE and for TM modes, i.e. the Fresnel coefficients for the two different states of polarization of the electric and magnetic fields. To get the formula for the phase-shift you only need to equal this complex expression for $R$ to the Fresnel formulas and solve for $\phi$. The formula you wrote is the one corresponding to TE modes, so it's the TE phase-shift (with electric fields perpendicular to the plane of incidence spanned by the wave normal and the normal to the interface). For TM modes you would get exactly the same result but multiplied by a factor of $n_{1}^{2}/n_{2}^{2}$.