#

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