Homec4science
Diffusion JuNoTE_git_repository (staging) Browse

R9442/notebooks/FresnelCoefficientsCode.py
2c31c87ef8c2staging

Last Change
Raw File

FresnelCoefficientsCode.py
View Options

import numpy as np
from ipywidgets import IntSlider, interact, FloatSlider, Text, widgets, interact_manual
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
params = {'legend.fontsize': 'x-large',
'figure.figsize': (10, 5),
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
'xtick.labelsize':'x-large',
'ytick.labelsize':'x-large',
'font.size':10}
pylab.rcParams.update(params)
from IPython.display import display
import matplotlib.ticker as tck
def snell(t1, n1=2, n2=1):
sint2 = n1 * np.sin(t1) / n2
if sint2 < 1:
return np.arcsin(sint2)
def plot_snell(n1, n2, t1):
n1 = float(n1)
n2 = float(n2)
t1= np.deg2rad(t1)
t2 = snell(t1, n1, n2)
x = np.linspace(-1, 0, 100)
y = np.linspace(0, 1, 100)
plt.polar(np.pi * np.ones(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.zeros(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.pi/2 * np.ones(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.ones(100) * 3 *np.pi/2, np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.text(np.pi/4, 0.2, '$n_1=%.0f$'%n1, fontsize=15)
plt.text(4.5*np.pi/4, 0.6, '$n_2=%.0f$'%n2, fontsize=15)
try:
plt.polar((np.pi / 2 + t1) * np.ones(100),
y,
c='k',
lw=3,
label='Incident')
plt.polar(np.linspace(np.pi / 2, np.pi / 2 + t1, 10),
0.5 * np.ones(10),
ls=':',
c='k',
lw=3,
label='$\\theta_1=%.2f^\circ$' % np.rad2deg(t1))
plt.polar((3 * np.pi / 2 + t2) * np.ones(100),
y,
c='b',
lw=3,
alpha = 1,
label='Refracted' )
plt.polar(np.linspace(3*np.pi / 2,3*np.pi / 2 + t2,10) ,
0.5 * np.ones(10),
ls=':',
c='b',
lw=3,
label='$\\theta_2=%.2f^\circ$'% np.rad2deg(t2))
except TypeError as te:
tc = np.arcsin(n2/n1)
plt.polar((np.pi / 2 - t1) * np.ones(100),
y,
c='r',
lw=3,
alpha = 1,
label='Reflection')
plt.polar(np.linspace(np.pi / 2,np.pi / 2 - t1,10) ,
0.5 * np.ones(10),
ls=':',
c='r',
lw=3,
label='$\\theta_1\'=%.2f^\circ$'% np.rad2deg(t1))
plt.polar((np.pi / 2 + tc) * np.ones(100),
y,
c='purple',
lw=3, ls='--',
label='Critical angle\n$\\theta_c=%.2f^{\circ}$'%np.rad2deg(tc))
plt.gca().set_yticklabels([])
plt.gca().grid(False)
plt.legend(loc=(1, 1))
from mpl_toolkits.mplot3d import Axes3D
def E(z, n):
'''Returns a static electric field waveform given a polatization vector n'''
return np.exp(1j * z) * n[0], np.exp(1j * z) * n[1]
def alpha_beta_to_theta(alpha, beta):
norm = np.sqrt(alpha**2 + beta**2)
alpha/=norm
beta/=norm
theta = np.arccos(alpha)
return theta
def n(theta, phase1, phase2):
'''Returns a polarization (Jones) vector given an and phases'''
n0 = np.array([
np.cos(theta) * np.exp(1j * phase1), #change by alpha
np.sin(theta) * np.exp(1j * phase2) #change by beta
])
norm = np.sqrt(n0.conj().dot(n0))
n0/=norm
return n0
def plot_polarization(alpha, beta, phase1, phase2):
theta = alpha_beta_to_theta(alpha, beta)
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122, projection='polar')
#theta = np.deg2rad(theta)
phase1 = np.deg2rad(phase1)
phase2 = np.deg2rad(phase2)
z = np.linspace(0, 4 * np.pi, 100)
jones = n(theta, phase1, phase2)
E_field = np.real(E(z, jones))
ax.plot(1.5 * np.ones(100), E_field[1], z, zdir='x', label='$E_2$')
ax.plot(E_field[0], -1 * np.ones(100), z, zdir='x', label='$E_1$')
ax.plot(E_field[0], E_field[1], 0, zdir='x')
ax.plot(E_field[0], E_field[1], z, lw=5, label='$\\vec{E}$', zdir='x')
ax.azim = -35
ax.elev = 20
ax.dist = 12
angle = np.real(np.arctan(jones[1] / jones[0]))
mag = np.sqrt(jones.conj().dot(jones))
ax2.plot(np.linspace(0, angle, 2), np.linspace(0, mag, 2), lw=5, c='k', label = '$\hat{n}$')
ax2.plot(np.linspace(0, angle, 10),
0.5 * mag * np.ones(10),
ls=':',
c='k',
lw=3,
label='$\\theta=%.2f^\circ$' % np.rad2deg(angle))
ax.set_xlabel('$\hat{e}_3$', fontsize=25)
ax.set_ylabel('$\hat{e}_1$', fontsize=25)
ax.set_zlabel('$\hat{e}_2$', fontsize=25)
ax.set_ylim(-1, 1.5) #labeled x
ax.set_xlim(0, max(z)) #labeled z
ax.set_zlim(-1, 2) #labeled y
ax.tick_params(pad=55)
ax.tick_params(axis='y', rotation=-45)
ax2.set_ylim(0, 1)
ax2.plot(np.arctan(E_field[1] / E_field[0]),
np.sqrt(E_field[0]**2 +E_field[1]**2 ), '.r')
ax2.plot(np.arctan(E_field[1] / E_field[0])+np.pi,
np.sqrt(E_field[0]**2 +E_field[1]**2 ), '.r', label='$|E|$')
ax.legend(fontsize=25)
ax2.legend(fontsize=25)
ax.set_title(r'$\vec{E} = \alpha~\exp(i\phi_1)~ E_1 \hat{e}_1 + \beta~\exp(i\phi_2)~ E_2 \hat{e}_2$', fontsize = 20)
plt.tight_layout()
def sen_t2(n1, n2, t1):
sin_t2 = n1 * np.sin(t1) / n2
return sin_t2
def rs(n1, n2, t1, t2=0):
sin_t2 = sen_t2(n1, n2, t1)
cos_t2 = np.sqrt(1 - sin_t2**2)
return (n1 * np.cos(t1) - n2 * cos_t2) / (n1 * np.cos(t1) + n2 * cos_t2)
def rp(n1, n2, t1, t2=0):
sin_t2 = sen_t2(n1, n2, t1)
cos_t2 = np.sqrt(1 - sin_t2**2)
return (n1 * cos_t2 - n2 * np.cos(t1)) / (n2 * np.cos(t1) + n1 * cos_t2)
def ts(n1, n2, t1, t2=0):
sin_t2 = sen_t2(n1, n2, t1)
cos_t2 = np.sqrt(1 - sin_t2**2)
return 2 * n1 * np.cos(t1)/(n1 * np.cos(t1) + n2 * cos_t2)
def tp(n1, n2, t1, t2=0):
sin_t2 = sen_t2(n1, n2, t1)
cos_t2 = np.sqrt(1 - sin_t2**2)
return 2 * n1 * np.cos(t1)/(n2 * np.cos(t1) + n1 * cos_t2)
def plot_fresnel(n1, n2):
fig = plt.figure()
n1 = float(n1)
n2 = float(n2)
if n2 < n1 :
tmax = np.arcsin(n2/n1)
plt.axvline(np.rad2deg(tmax), label = 'Critical angle')
else:
tmax = np.pi/2
t1 = np.linspace(0,tmax, 100)
t1_deg = np.rad2deg(t1)
plt.plot(t1_deg, rs(n1, n2, t1), label = '$r_s$', c = 'b', ls = '-')
plt.plot(t1_deg, ts(n1, n2, t1), label = '$t_s$', c = 'b', ls = '-.')
plt.plot(t1_deg, rp(n1, n2, t1), label = '$r_p$', c = 'r', ls = '-')
plt.plot(t1_deg, tp(n1, n2, t1), label = '$t_p$', c = 'r', ls = '-.')
tb = np.rad2deg(np.arctan(n2/n1))
plt.axvline(tb, c= 'k', ls=':', label = '$\\theta_B = %.1f ^\circ$'%tb)
plt.legend(loc=0)
plt.grid()
plt.xlim(0,90)
plt.xlabel('$\\theta_1$ [$^\circ$]')
def plot_fresnel_RT(n1, n2):
n1 = float(n1)
n2 = float(n2)
if n2 < n1:
tmax = np.arcsin(n2 / n1)
plt.axvline(np.rad2deg(tmax), label='Critical angle')
else:
tmax = np.pi / 2
t1 = np.linspace(0, tmax, 100)
t1_deg = np.rad2deg(t1)
cos_t2 = np.sqrt(1 - sen_t2(n1, n2, t1)**2)
cos_t1 = np.cos(t1)
T_s = (n2 * cos_t2 / (n1 * cos_t1)) * np.abs(ts(n1, n2, t1))**2
T_p = (n2 * cos_t2 / (n1 * cos_t1)) * np.abs(tp(n1, n2, t1))**2
R_s = np.abs(rs(n1, n2, t1))**2
R_p = np.abs(rp(n1, n2, t1))**2
plt.plot(t1_deg, R_s, label='$R_s$', c='b', ls='-')
plt.plot(t1_deg, T_s,
label='$T_s$',
c='b',
ls='-.')
plt.plot(t1_deg, R_p, label='$R_p$', c='r', ls='-')
plt.plot(t1_deg, T_p,
label='$T_p$',
c='r',
ls='-.')
tb = np.rad2deg(np.arctan(n2 / n1))
plt.axvline(tb, c='k', ls=':', label='$\\theta_B = %.1f ^\circ$' % tb)
plt.legend(loc=0)
plt.grid()
plt.xlim(0, 90)
plt.ylim(0, 1)
plt.xlabel('$\\theta_1$ [$^\circ$]')
def plot_snell_ref(n1, n2, t1, alpha_refl, alpha_refr):
n1 = float(n1)
n2 = float(n2)
x = np.linspace(-1, 0, 100)
y = np.linspace(0, 1, 100)
t1= np.deg2rad(t1)
t2 = snell(t1, n1, n2)
plt.polar(np.pi * np.ones(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.zeros(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.pi/2 * np.ones(100), np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.polar(np.ones(100) * 3 *np.pi/2, np.linspace(0, 1, 100), c = 'k', ls = '--')
plt.text(np.pi/4, 0.2, '$n_1=%.0f$'%n1, fontsize=15)
plt.text(4.5*np.pi/4, 0.6, '$n_2=%.0f$'%n2, fontsize=15)
try:
plt.polar((np.pi / 2 + t1) * np.ones(100),
y,
c='k',
lw=5,
label='Incident')
plt.polar(np.linspace(np.pi / 2, np.pi / 2 + t1, 10),
0.5 * np.ones(10),
ls=':',
c='k',
lw=5,
label='$\\theta_1=%.2f^\circ$' % np.rad2deg(t1))
lrefr, = plt.polar((3 * np.pi / 2 + t2) * np.ones(100),
y,
c='b',
lw=5,
alpha = alpha_refr,
label='Refracted' )
lrefr.set_alpha(alpha_refr)
plt.polar(np.linspace(3*np.pi / 2,3*np.pi / 2 + t2,10) ,
0.5 * np.ones(10),
ls=':',
c='b',
lw=5,
label='$\\theta_2=%.2f^\circ$'% np.rad2deg(t2))
lrefl, = plt.polar((np.pi / 2 - t1) * np.ones(100),
y,
c='r',
lw=5,
alpha = alpha_refr,
label='Reflected' )
lrefl.set_alpha(alpha_refl)
plt.polar(np.linspace(np.pi / 2,np.pi / 2 - t1,10) ,
0.5 * np.ones(10),
ls=':',
c='r',
lw=5,
label='$\\theta_1\'=%.2f^\circ$'% np.rad2deg(t1))
except TypeError as te:
tc = np.arcsin(n2/n1)
plt.polar((np.pi / 2 - t1) * np.ones(100),
y,
c='r',
lw=5,
alpha = 1,
label='Reflection')
plt.polar(np.linspace(np.pi / 2,np.pi / 2 - t1,10) ,
0.5 * np.ones(10),
ls=':',
c='r',
lw=5,
label='$\\theta_2=%.2f^\circ$'% np.rad2deg(t1))
plt.polar((np.pi / 2 + tc) * np.ones(100),
y,
c='b',
lw=5, ls='--',
label='Critical angle\n$\\theta_c=%.2f^{\circ}$'%np.rad2deg(tc))
plt.gca().set_yticklabels([])
plt.gca().grid(False)
plt.legend(loc=(1, 1))
def wrap_plot_snell_ref(n1, n2, t1, sorp = 's'):
n1, n2, t1 = float(n1), float(n2), float(t1)
t1_rad = np.deg2rad(t1)
cos_t2 = np.sqrt(1 - sen_t2(n1, n2, t1_rad)**2)
cos_t1 = np.cos(t1_rad)
T_p = (n2 * cos_t2 / (n1 * cos_t1)) * np.abs(tp(n1, n2, t1_rad))**2
R_p = np.abs(rp(n1, n2, t1_rad))**2
if sorp == 's':
T = (n2 * cos_t2 / (n1 * cos_t1)) * np.abs(ts(n1, n2, t1_rad))**2
R = np.abs(rs(n1, n2, t1_rad))**2
elif sorp == 'p':
T = (n2 * cos_t2 / (n1 * cos_t1)) * np.abs(tp(n1, n2, t1_rad))**2
R = np.abs(rp(n1, n2, t1_rad))**2
else:
raise NotImplementedError
plot_snell_ref(n1, n2, t1, R,T)
c4science · Help