#
# Rainbows_Exercise_5.py
#
# Plot the irradiance as a function of
# deflection angle for rays undergoing
# a single internal reflection,
# assuming a spherical water drop
# and two wavelengths.
#
# Here, as a 'zeroeth' approximation,
# the Fresnel coefficients are NOT taken into account.
#
# 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,arcsin,sin,exp,zeros,argmin,argmax
# 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 = 801 # the number of scaled impact parameters, i.e., rays
width = 1.0 # angular width of ray [deg]
ndetpts = 501 # the number of detector points
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
# Generate the array of uniformly distributed scaled impact parameters
scaled_b = linspace(0.0,1.0,npts)
# Generate the corresponding array of incident angle values [rad]
theta_i = arcsin(scaled_b)
# Generate the corresponding array of incident angle values [deg]
theta_i_deg = 180.0*theta_i/pi
# 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
theta1r = 2.0*(theta_i - theta_r) + 1.0*(pi - 2.0*theta_r) # First deflection function
theta1r_deg = theta1r*180.0/pi
# Calculate the deflection angle for each incident angle for violet light
theta1v = 2.0*(theta_i - theta_v) + 1.0*(pi - 2.0*theta_v) # First deflection function
theta1v_deg = theta1v*180.0/pi
print("The primary rainbow deflection angle for red light is ",min(theta1r_deg))
# The index of this value is
index_r = argmin(theta1r_deg)
print("This primary 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 primary rainbow appears at ",180.0-min(theta1r_deg)," above the horizontal.")
print(" ")
print("The primary rainbow deflection angle for violet light is ",min(theta1v_deg))
# The index of this value is
index_v = argmin(theta1v_deg)
print("This primary 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 primary rainbow appears at ",180.0-min(theta1v_deg)," above the horizontal.")
print(" ")
# Set up the array of irradiance versus deflection angle
deflection_angles_deg = linspace(130,180,ndetpts)
irradiance_r = zeros(ndetpts)
irradiance_v = zeros(ndetpts)
# Accummulate irradiance from each ray
# that corresponds to one impact parameter
for i in range(0,ndetpts):
# for every detection angle...
for j in range(0,npts):
# accumulate the irradiance into a particular angle from each ray
irradiance_r[i] = irradiance_r[i] + exp(-((deflection_angles_deg[i] - theta1r_deg[j])/width)**2)
irradiance_v[i] = irradiance_v[i] + exp(-((deflection_angles_deg[i] - theta1v_deg[j])/width)**2)
# Scale the intensities
scaled_irradiance_r = irradiance_r/max(irradiance_r)
scaled_irradiance_v = irradiance_v/max(irradiance_v)
# The index of this value at which the maximum scaled irradiance occurs is
index_r_irradiance = argmax(scaled_irradiance_r)
print("The detection angle at which ")
print("the red irradiance is maximum = ",deflection_angles_deg[index_r_irradiance]," deg.")
print("The maximum appears ",180.0-deflection_angles_deg[index_r_irradiance], " deg. above the horizontal.")
print(" ")
# The index of this value is
index_v_irradiance = argmax(scaled_irradiance_v)
print("The detection angle at which ")
print("the violet irradiance is maximum = ",deflection_angles_deg[index_v_irradiance]," deg.")
print("The maximum appears ",180.0-deflection_angles_deg[index_v_irradiance], " deg. above the horizontal.")
print(" ")
#----------------------------------------------------------------------
# Start a new figure. This will be a plot of the
# scaled irradiance versus deflection angles
figure()
# Set the limits of the horizontal axis
xlim(130,180)
# Label the horizontal axis
xlabel("Detection angle relative to incident ray $\\phi$ [deg]", size = 16)
# Set the limits of the vertical axis
ylim(0.0,1.0)
# Label the vertical axis
ylabel("Scaled irradiance $I/I_{max}$", size = 16)
# Draw a grid
grid(True)
# Plot the deflection functions
plot(deflection_angles_deg,scaled_irradiance_r,"r-",label="$\lambda = 650 \\, {\\rm nm}$")
plot(deflection_angles_deg,scaled_irradiance_v,"b-",label="$\lambda = 400 \\, {\\rm nm}$")
# Generate the legend
legend(loc=1)
# Generate the title
title("Crude estimate of irradiance versus deflection angle\n for one internal reflection")
show()