The boundary conditions for the electric field are
\[ \begin{cases}
1+r = A+B \\
A e^{in_1kd} + B e^{-in_1kd} = t
\end{cases}.\]
At the same time, as the magnetic field $\vec{B} = \vec{k} \times \vec{E}/\omega$ we have boundary conditions for the magnetic field
\[\begin{cases}
1 - r = n_1 (A - B) \\
n_1 (A e^{in_1kd} - B e^{-in_1kd} ) = n_2 t
\end{cases}.\]
First, let's exclude $t$:
\[
A e^{in_1kd} + Br^{-in_1kd} = \frac{n_1}{n_2} \left( A e^{in_1kd} - B e^{-in_1kd} \right) \quad \Rightarrow \quad \frac{B}{A} = -\frac{n_2-n_1}{n_2 + n_1} e^{2in_1kd} .\]
Also
\[\begin{cases}
2 = A (n_1 + 1) + B (1-n_1) \\
2r = A(1 - n_1) + B(1 + n_1)
\end{cases} \quad \Rightarrow \quad r = \frac{-\frac{n_1-1}{n_1+1}+ \frac{B}{A}}{ 1-\frac{B}{A}\frac{n_1 - 1}{n_1+1}}.\]
Substituting the eqution which connects $A$ and $B$ gives us
\[r=\frac{-\frac{n_1-1}{n_1+ 1} - \frac{n_2-n_1}{n_2+n_1} e^{2in_1 kd} }{1+\frac{n_2 - n_1}{n_2 + n_1} \frac{n_1 - 1}{n_1 + 1} e^{2 in_1 kd}},\]
that could be easly transform to the required form for instance with substituting
\[ a= \frac{n_1-n_2}{n_1+n_2}, \quad b = \frac{1-n_1}{1+n_1}, \quad \phi = -2n_1kd.\]
Interesting fact about the complex number
\[z=\frac{a+be^{i\phi}}{ab+ e^{i\phi}}\]
is that if we try to represent $z(\phi)$ as a curve in the coordinates $\operatorname{\mathfrak{Im}} z$ vs. $\operatorname{\mathfrak{Re}} z$, we get a circle. It could be easly prooved with algebra of complex numbers. Let's assume that
\[z_0 = \frac{1}{2} \left( \frac{a+b}{ab+1} + \frac{a-b}{ab-1} \right) = \frac{a^2b - b}{a^2b^2-1}\]
is a center of this circle. Thus
\[z-z_0 = \frac{a+be^{i\phi}}{ab+e^{i\phi}} - \frac{a^2b - b}{a^2b^2-1} = \frac{a(b^2-1)(1+abe^{i\phi})}{(a^2b^2-1)(ab+e^{i\phi})}\]
and
\[|z-z_0| = \left| \frac{a(b^2-1)}{(a^2b^2-1)} \right| \sqrt{\frac{1+abe^{i\phi}+a^*b^*e^{-i\phi}+aa^*bb^*}{aa^*bb^*+abe^{-i\phi}+a^*b^*e^{i\phi}+1}}= \left| \frac{a(b^2-1)}{(a^2b^2-1)} \right|\]
is a constant, so our assumption is correct. To find the refractive index let's find a ratio $\Gamma$ between maximal and minimal value of the reflection index $R$:
\[\Gamma = \frac{(|z_0|+|z-z_0|)^2}{(|z_0|-|z-z_0|)^2}\]
For the given data $I_\text{max}=391.0~\text{a.u.}$ and $I_\text{min}=54.5~\text{a.u.}$, so $\Gamma=7.17$. From the graph above we can obtain that $n_1=1.64$. You can sertianly obtain the same result manually without the following analysis of $R$ with the program thin_film.py.
For $n_1=1.64$ a graph of $R$ vs. $d$ could be plotted with thin_film.py. The minimum reflection is reached for $d_\text{min}=75.8~\text{nm}$ and for $N_\text{min}=711$, so we can estimate $d_0 = d_\text{min}/N_\text{min}=0.107~\text{nm}$.
To check the value of $d_0$ we can plot a graph $I(N)$ vs $R(Nd_0)$ and it should be linear. According to the error in $N_\text{min}$ that we just got from the data, the best fit corresponds to $d_0=0.110~\text{nm}$