#

# 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()