Question: 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?

Answer: 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")