#
# Rainbows_Exercise_4.py
#
# Plot the deflection function
# for the secondary (second-order) rainbow
# assuming a spherical water drop
# and a single wavelength,
# which implies a single index of refraction
#
# The refractive index of water is taken from
# Eq. 3 of P.D.T Huibers, Applied Optics,
# Vol. 36, No. 16, pp. 3785-3787 (1997).
#
# The refractive index of air is taken from
# Eq. 1 of P.E. Ciddor, Applied Optics,
# Vol. 35, No. 9, pp. 1566-1573 (1996).
#
# Written by:
#
# Ernest R. Behringer
# Department of Physics and Astronomy
# Eastern Michigan University
# Ypsilanti, MI 48197
# (734) 487-8799
# ebehringe@emich.edu
#
# 20160109 by ERB
#
from __future__ import print_function
from pylab import figure,plot,xlim,xlabel,ylim,ylabel,grid,legend,title,show
from math import pi,asin
from numpy import linspace,sin,zeros,argmin
# Define the index of refraction function for water
def water_index(wavelength):
n_H2O = 1.31279 + 15.762/wavelength - 4382.0/(wavelength**2) + 1.1455e6/(wavelength**3)
return(n_H2O)
# Define the index of refraction function for air
# Note: wavelength is supposed to be in micrometers
def air_index_minus_one(wavelength):
term1 = 0.05792105/(238.0185 - 1.0e6/(wavelength**2))
term2 = 0.00167917/(57.362 - 1.0e6/(wavelength**2))
return(term1+term2)
# Inputs
wavelength_r = 650 # vacuum wavelength in nm ("red")
wavelength_v = 400 # vacuum wavelength in nm ("violet")
n_r = water_index(wavelength_r) # refractive index of water at wavelength_r
n_v = water_index(wavelength_v) # refractive index of water at wavelength_v
npts = 901
n_ar = 1.0 + air_index_minus_one(wavelength_r) # refractive index of air at wavelength_r
n_av = 1.0 + air_index_minus_one(wavelength_v) # refractive index of air at wavelength_v
theta_i_deg = linspace(0.0,90.0,npts) # array of incident angle values [deg]
theta_i = theta_i_deg*pi/180.0 # array of incident angle values [rad]
# Set up the arrays of deflection angles
theta_r = zeros(npts)
theta_v = zeros(npts)
# Calculate the refraction angle for each incident angle
# and for red and for violet light
for j in range(0,npts):
theta_r[j] = asin((n_ar/n_r)*sin(theta_i[j]))
theta_v[j] = asin((n_av/n_v)*sin(theta_i[j]))
# Calculate the deflection angle for each incident angle for red light
theta2r = 2.0*(theta_i - theta_r) + 2.0*(pi - 2.0*theta_r) # First deflection function
theta2r_deg = theta2r*180.0/pi
# Calculate the deflection angle for each incident angle for violet light
theta2v = 2.0*(theta_i - theta_v) + 2.0*(pi - 2.0*theta_v) # First deflection function
theta2v_deg = theta2v*180.0/pi
print("The secondary rainbow deflection angle for red light is ",min(theta2r_deg))
# The index of this value is
index_r = argmin(theta2r_deg)
print("This secondary rainbow deflection angle for red light corresponds to")
print("an incident angle of ",theta_i_deg[index_r]," deg.")
print("The red band of the secondary band appears at ",min(theta2r_deg)-180.0," above the horizontal.")
print(" ")
print("The secondary rainbow deflection angle for violet light is ",min(theta2v_deg))
# The index of this value is
index_v = argmin(theta2v_deg)
print("This secondary rainbow deflection angle for violet light corresponds to")
print("an incident angle of ",theta_i_deg[index_v]," deg.")
print("The violet band of the secondary band appears at ",min(theta2v_deg)-180.0," above the horizontal.")
#----------------------------------------------------------------------
# Start a new figure. This will be a plot of the deflection functions,
# i.e., deflection angle versus incident angle
figure()
# Set the limits of the horizontal axis
xlim(0,90)
# Label the horizontal axis
xlabel("$\\theta_i$ [deg]", size = 16)
# Set the limits of the vertical axis
ylim(220,360)
# Label the vertical axis
ylabel("$\\delta$ [deg]", size = 16)
# Draw a grid
grid(True)
# Plot the deflection functions
plot(theta_i_deg,theta2r_deg,"r-",label="$\lambda = 650 \\, {\\rm nm}$")
plot(theta_i_deg,theta2v_deg,"b-",label="$\lambda = 400 \\, {\\rm nm}$")
# Generate the legend
legend(loc=1)
# Generate the title
title("Deflection function for two internal reflections")
show()