Maxwell's equations in the dielectric material are
\[
\DeclareMathOperator{\Rot}{curl}
\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\Grad}{grad}
\begin{split}
&\Div \vec{D} = 0 \quad \Rot \vec{E} = -\frac{\partial \vec{B}}{\partial t}\\
&\Div \vec{B} = 0 \quad \Rot \vec{B} = \mu_0 \frac{\partial \vec{D}}{\partial t}
\end{split},
\]
where $\vec{D}$ is the electric displacement field related to the electric field $E$ by the $\vec{D} = \varepsilon_0 \varepsilon \vec{E}$, where $\varepsilon$ is a relative permittivity of the material. Keep in mind that $c^2 \varepsilon_0 \mu_0=1$, where $c$ is the speed of light.
Maxwell's equation could be solved by plane waves
\[
\begin{split}
\vec{E}(\vec{r}, t) = \mathfrak{Re} \vec{E}_0 e^{i\vec{k}\vec{r} - i \omega t} \\
\vec{B}(\vec{r}, t) = \mathfrak{Re} \vec{B}_0 e^{i\vec{k}\vec{r} - i \omega t}
\end{split},\]
where $\mathfrak{Re}$ is a real part operator, $\vec{k}$ is a wavevector, $\omega$ is a wave frequency, $\vec{E}_0$ and $\vec{B}_0$ are amplitudes of oscillations of the electric and magnetic fields, respectively.
Remark: The divergence $\Div \vec{A}$ of a vector field $\vec{A}(\vec{r})$ is a scalar value. It has the meaning of the ratio of the flux through the surface of a tiny volume to its volume and is given by
\[ \Div \vec{A} = \frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y} + \frac{\partial A_z}{\partial z}.\]
The $\Rot \vec{A}$ of a vector field $\vec{A}(\vec{r})$ is a vector variable and each of its components means the ratio of the circulation around some shape to its area. It's given by
\[ \Rot \vec{A} = \begin{vmatrix}
\vec{x} & \vec{y} & \vec{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\
A_x & A_y & A_z
\end{vmatrix}
=
\vec{x}
\begin{vmatrix}
\frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\
A_y & A_z
\end{vmatrix}
-
\vec{y}
\begin{vmatrix}
\frac{\partial}{\partial x} & \frac{\partial}{\partial z} \\
A_x & A_z
\end{vmatrix}
+
\vec{z}
\begin{vmatrix}
\frac{\partial}{\partial x} & \frac{\partial}{\partial y} \\
A_x & A_y
\end{vmatrix}
=
\begin{pmatrix}
\frac{\partial A_z}{\partial y} - \frac{\partial A_y}{\partial z} \\
-\left[ \frac{\partial A_z}{\partial x} - \frac{\partial A_x}{\partial z} \right] \\
\frac{\partial A_y}{\partial x} - \frac{\partial A_x}{\partial y} \\
\end{pmatrix}.\]
The Lapplacian $\Delta \vec{A}$ of the vector field $\vec{A}$ is a vector which is given by
\[\Delta \vec{A} = \frac{\partial^2 \vec{A}}{\partial x^2} +
\frac{\partial^2 \vec{A}}{\partial y^2} +
\frac{\partial^2 \vec{A}}{\partial z^2} \]
The operators $\mathfrak{Re}$ and $\mathfrak{Im}$ return only the real and imaginary parts of a complex number:\[ \mathfrak{Re} \lbrace a+bi \rbrace= a, \quad \mathfrak{Im} \lbrace a + bi \rbrace= b\]
These operations are linear, so you can interchange it with a differentiation: $ \left( \mathfrak{Re} \lbrace f(x) \rbrace \right)' = \mathfrak{Re} \lbrace f'(x) \rbrace$.
\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\Grad}{grad} \Rot \Rot \vec{A} = \Grad \Div \vec{A} - \Delta \vec{A}$ find a wave equation for the vector $\vec{E}$:
\[ \Delta \vec{E} - \frac{1}{v^2}\frac{\partial^2 E}{\partial t^2} = 0.\]
What is the phase velocity $v$ of the electromagnetic wave in the medium with relative permittivity $\varepsilon$? What is the refraction index $n$ of the medium with relative permittivity $\varepsilon$?
Remark: partial derivatives with respect to different variables could be interchanged:
\[ \frac{\partial^2 f(x,y)}{\partial x \partial y} = \frac{\partial^2 f(x,y)}{\partial y \partial x} \quad \Rightarrow \quad \Rot \frac{\partial}{\partial t} \vec{A} = \frac{\partial}{\partial t} \Rot \vec{A}\]
To describe the electromagnetic wave, it's usually convenient to separate the time and coordinate dependence of the electric field by introducing a vector $\tilde{E}(\vec{r})$ of the complex amplitude:
\[ \vec{E}(\vec{r}, t) = \mathfrak{Re} \left\lbrace \tilde{E}(\vec{r}) e^{- i \omega t} \right\rbrace.\]
You can make sure that the Helmoltz equation holds for the complex amplitude $\tilde{E}$:
\[ \Delta \tilde{E} + \frac{\omega^2}{v^2} \tilde{E} = 0.\]
Its solution is a coordinate dependent part of the plane wave $\vec{A} e^{i\vec{k}\cdot \vec{r}}$.
Consider a plane wave with the angle of incidence $\theta_1$ at the plane interface between two dielectrics. Then the complex amplitude is
\[
\tilde{E}(\vec{r}) = \begin{cases}
\vec{p}_0 E_0 e^{i k_{1y}y + i k_{1x}x} + r \vec{p}_r E_0 e^{i k_{1y}y - i k_{1x}x},& \text{ for } x < 0\\
t \vec{p}_t E_0 e^{i k_{2y}y + i k_{2x}x},& \text{ for } x > 0
\end{cases},\]
where $\vec{p}_0$, $\vec{p}_r$, $\vec{p}_t$ are unit vectors that determine the direction of the electric field oscillations, $E_0$ is the amplitude of the electric field oscillations in the incident wave, $r$ and $t$ are the amplitude coefficients of reflection and transmission, respectively.
There are two different cases of polarization: $s$- and $p$-polarized waves. The electric field $\vec{E}$ is parallel to the surface for $s$-polarization and $\vec{p}_0 = \vec{p}_r = \vec{p}_t$. The electric field $\vec{E}$ is in the plane of incidence for $p$-polarization and the vectors $\vec{p}_0$, $\vec{p}_r$, $\vec{p}_t$ are chosen so that their cross product with the wave vectors has the same direction.
The next consequence of the Maxwell's equations are boundary conditions at the plane interface between two dielectrics with relative permittivities $\varepsilon_1$ and $\varepsilon_2$. The $1,2$ indices for each variable indicate the medium we are describing. From the Maxwell's equations with divergence it follows that the normal to the surface components are equal in pairs $\vec{D}_1$, $\vec{D}_2$ and $\vec{B}_1$, $\vec{B}_2$:
\[ D_{1\perp}=D_{2\perp}, \quad B_{1\perp}=B_{2\perp}.\]
From the Maxwell''s equations with curl it follows that the planar components are equal in pairs $\vec{E}_{1,2}$ and $\vec{B}_{1,2}$:
\[E_{1\parallel} = E_{2\parallel}, \quad B_{1\parallel} = B_{2\parallel}.\]
Vectors $\vec{E}$, $\vec{D}$ and $\vec{B}$ are strictly related to the complex amplitude $\tilde{E}$, so the discussed boundary conditions could be rewritten in terms of the complex amplitude. Note, that the frequency can't change during reflection or transmission because it's impossible to satisfy boundary conditions in such a case.
Let $\theta_2$ the angle of refraction.
\[ r_s = \frac{n_1 \cos \theta_1 - n_2 \cos \theta_2}{n_1 \cos \theta_1+ n_2 \cos \theta_2}\]
\[ r_p = \frac{n_2 \cos \theta_1 - n_1\cos \theta_2}{n_2 \cos \theta_1+ n_1 \cos \theta_2}\]
The light intensity $I$ is the average energy flux brought by the radiation. It can be found as the absolute value of the average over the period of the oscillations of the electric field of the Poynting vector:
\[ \vec{S} = \frac{1}{\mu_0} \vec{E} \times \vec{B}\]
\[ I = \frac{c \varepsilon_0 n}{2}|\tilde{E}|^2,\]
where $n$ is the refractive index of the medium.
The reflection coefficient $R$ is the ratio of the intensity of the reflected waves to the intensity of the incident wave.
With fresnel.py you can plot $R_p/R_s$ as a function of incident angle $\theta_1$ for different $\varepsilon_1$ and $\varepsilon_2$.
The figure below shows the graph of $R_p/R_s$ for the data obtained in air ($n_1=1$) for two unknown materials $A$ and $B$ with unknown refractive indices $n_A$ and $n_B$, respectively. Note, that the refractive index of metamaterials can be almost any.
