python pandas matplotlib reflection fresnel

Why is the plot of the refractive index wavelength-dependent Fresnel's equation not showing as expected?

I want to reproduce the reflectance spectrum of a thin film whose complex refractive index is wavelength-dependent (the complex refractive index data, named N2 in code, can be obtained from here).

Using fresnel's equations for medium 1: air, medium 2: thin film, and medium 3: air : ;

The reflection coefficient

where \delta is defined as

And finally the reflectance equation

Since the refractive indices have a dependence on the wavelength, then the refraction angle have the same dependence, using the Snell equation:

I wrote the following code:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

## Function definition
def Ang_Refrac(na,nb,angle_a):
ang_refrac = np.arcsin( (na/nb).real * np.sin(angle_a) )
return ang_refrac

# Data
Thin_film = pd.read_csv("Si_refractive index.txt", delimiter='\t')

## Variable definition

# thin film thickness, [d]: nm
d = 200

# Wave-length
lamb = Thin_film[Thin_film.columns[0]]

# Complex refractive index
N1 = np.ones(len(lamb))
N2 = Thin_film[Thin_film.columns[1]] + Thin_film[Thin_film.columns[1]]*1j
N3 = np.ones(len(lamb))

# Angle:
ang_1_s = 0  #sexagesimal

ang_1 = ang_1_s*np.pi/180  #radians
ang_2 = [] #radians
ang_3 = [] #radians
for i in range(len(N2)):
    ang_refrac_12 = Ang_Refrac(N1[i],N2[i],ang_1)
    ang_refrac_23 = Ang_Refrac(N2[i],N3[i],ang_refrac_12)
    ang_2.append(ang_refrac_12)
    ang_3.append(ang_refrac_23)

## Reflectance
R_s = np.zeros(len(lamb))
R_p = np.zeros(len(lamb))

for i in range(len(lamb)):
    # S-type polarization
    r12_s = (N1[i]*np.cos(ang_1) - N2[i]*np.cos(ang_2[i]))/(N1[i]*np.cos(ang_1) + N2[i]*np.cos(ang_2[i]))
    r23_s = (N2[i]*np.cos(ang_2[i]) - N3[i]*np.cos(ang_3[i]))/(N2[i]*np.cos(ang_2[i]) + N3[i]*np.cos(ang_3[i]))

    # P-type polarization
    r12_p = (N2[i]*np.cos(ang_1) - N1[i]*np.cos(ang_2[i])) / (N2[i]*np.cos(ang_1) + N1[i]*np.cos(ang_2[i]))
    r23_p = (N3[i]*np.cos(ang_2[i]) - N2[i]*np.cos(ang_3[i])) / (N3[i]*np.cos(ang_2[i]) + N2[i]*np.cos(ang_3[i]))

    # Phase shift
    delta = 2 * np.pi * (1/lamb[i]) * d * np.sqrt(abs(N2[i])**2 - (np.sin(ang_1)**2))

    # Reflection coefficient
    r123_s = (r12_s + r23_s*np.exp(2j*delta)) / (1 - r12_s*r23_s*np.exp(2j*delta))
    r123_p = (r12_p + r23_p*np.exp(2j*delta)) / (1 - r12_p*r23_p*np.exp(2j*delta))

    # Reflectance
    R_s[i] = abs(r123_s)**2
    R_p[i] = abs(r123_p)**2


# Reflectance normalization
R_s_normalized = (R_s - min(R_s)) / (max(R_s) - min(R_s))
R_p_normalized = (R_p - min(R_p)) / (max(R_p) - min(R_p))

# Plotting
plt.title("Refletância $\phi=" + str(ang_1_s) + "$°")
plt.xlim(300, 1000)  # Definindo os limites do eixo x
plt.plot(lamb, R_s_normalized, 'go', label="Polarização S", markersize=4)
plt.plot(lamb, R_p_normalized, 'b-', label="Polarização P")
plt.legend()
plt.show()

Resulting graph:

Expected graph: Obtained from here

Why are the results not the same?

Solution

You can try this. The major changes are:

N2 is given by n + i.k, rather than n + i.n as you had
there is just N2 ** 2 in the expression for delta, not abs(N2)**2, which would lose an important complex part
I think the reflection coefficients should have a '+', not a '-' sign in the denominator (though I could be wrong).

I found a little bit of background (and a few typos) in https://ir.lib.nycu.edu.tw/bitstream/11536/27329/1/000186822000019.pdf

Code:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

## Function definition
def Ang_Refrac(na,nb,angle_a):
    ang_refrac = np.arcsin( (na/nb) * np.sin(angle_a) )
    return ang_refrac

# Data
Thin_film = pd.read_csv("Si_refractive index.txt", delimiter='\t')

d = 200                                         # thin film thickness, [d]: nm
lamb = Thin_film[Thin_film.columns[0]]          # Wave-length

# Complex refractive index       n + ik
num = len(lamb)
N1 = np.ones(num)
N2 = Thin_film[Thin_film.columns[1]] + Thin_film[Thin_film.columns[2]]*1j      # <=========== NOTE ======
N3 = np.ones(num)

# Angle:
ang_1_s = 0
ang_1 = ang_1_s * np.pi / 180
ang_2 = [] #radians
ang_3 = [] #radians
for i in range(num):
    ang_refrac_12 = Ang_Refrac(N1[i],N2[i],ang_1)
    ang_refrac_23 = Ang_Refrac(N2[i],N3[i],ang_refrac_12)
    ang_2.append(ang_refrac_12)
    ang_3.append(ang_refrac_23)

## Reflectance
R_s = np.zeros(num)
R_p = np.zeros(num)

for i in range(num):
    # S-type polarization
    r12_s = (N1[i]*np.cos(ang_1) - N2[i]*np.cos(ang_2[i]))/(N1[i]*np.cos(ang_1) + N2[i]*np.cos(ang_2[i]))
    r23_s = (N2[i]*np.cos(ang_2[i]) - N3[i]*np.cos(ang_3[i]))/(N2[i]*np.cos(ang_2[i]) + N3[i]*np.cos(ang_3[i]))

    # P-type polarization
    r12_p = (N2[i]*np.cos(ang_1) - N1[i]*np.cos(ang_2[i])) / (N2[i]*np.cos(ang_1) + N1[i]*np.cos(ang_2[i]))
    r23_p = (N3[i]*np.cos(ang_2[i]) - N2[i]*np.cos(ang_3[i])) / (N3[i]*np.cos(ang_2[i]) + N2[i]*np.cos(ang_3[i]))

    # Phase shift
    delta = 2 * np.pi * ( d / lamb[i] ) * np.sqrt( N2[i]**2 - np.sin(ang_1)**2)       # <====== NOTE ====

    # Reflection coefficient
    r123_s = (r12_s + r23_s*np.exp(2j*delta)) / (1 + r12_s*r23_s*np.exp(2j*delta))    # <====== NOTE ====
    r123_p = (r12_p + r23_p*np.exp(2j*delta)) / (1 + r12_p*r23_p*np.exp(2j*delta))    # <====== NOTE ====

    # Reflectance
    R_s[i] = abs( r123_s ) ** 2
    R_p[i] = abs( r123_p ) ** 2


# Reflectance normalization
R_s_normalized = ( R_s - min( R_s ) ) / ( max( R_s ) - min( R_s ) )
R_p_normalized = ( R_p - min( R_p ) ) / ( max( R_p ) - min( R_p ) )

# Plotting
plt.title( r"Refletancia $\phi=" + str(ang_1_s) + "$ deg")
plt.xlim(300, 1000)  # Definindo os limites do eixo x
plt.plot(lamb, R_s_normalized, 'go', label="Polarizacao S", markersize=4)
plt.plot(lamb, R_p_normalized, 'b-', label="Polarizacao P")
plt.legend()
plt.show()

Output (note: some systematic difference in the normalisation from your "expected")