#
# Shadows_Exercise_3.py
#
# A linear array of N point sources located a specified distance from
# an L-shaped aperture composed of two rectangular apertures,
# one of which is centered on the origin.
# A screen is located a specified distance from
# the aperture.
#
# This file will generate a filled contour plot of
# the irradiance of light reaching the screen
# versus lateral coordinates (x,y)
#
# Written by:
#
# Ernest R. Behringer
# Department of Physics and Astronomy
# Eastern Michigan University
# Ypsilanti, MI 48197
# (734) 487-8799
# ebehringe@emich.edu
#
# 20160112-13 by ERB one aperture
# 20160119 by ERB compound aperture
# 20160609 clean up by ERB
#
# import the commands needed to make the plot
from pylab import axis,xlim,ylim,xlabel,ylabel,show,contourf,colorbar,figure,title
from matplotlib import cm
# import the command needed to make a 1D array
from numpy import meshgrid,absolute,where,zeros,linspace
# inputs
zs = -20.0 # distance between the source and aperture [cm]
zsc = 40.0 # distance between the screen and aperture [cm]
nso = 101 # number of point sources
length_so = 4.0 # length of the array of point sources [cm]
ap1_width = 2.0 # aperture width [cm]
ap1_height = 1.0 # aperture height [cm]
ap2_width = 1.0 # second aperture width [cm]
ap2_height = 3.0 # second aperture height [cm]
ap2_x = -0.5 # second aperture center x-coordinate [cm]
ap2_y = 2.0 # second aperture center y coordinate [cm]
screen_width = 60.0 # screen width [cm]
screen_height = 60.0 # wcreen height [cm]
nw = 240 # Number of screen width intervals
nh = 240 # Number of screen height intervals
# initialize needed arrays to zero
rso = zeros((nso,2)) # (x,y) location of the source [cm,cm]
screen_x = zeros((nw+1,nh+1)) # x-coordinates of screen points
screen_y = zeros((nw+1,nh+1)) # y coordinates of screen points
rsq = zeros((nw+1,nh+1)) # r squared values
irradiance = zeros((nw+1,nh+1)) # irradiance values
x0 = zeros((nw+1,nh+1)) # x coordinate of intersection at aperture plane
y0 = zeros((nw+1,nh+1)) # y coordinate of intersection at aperture plane
# create 1D arrays to create a meshgrid for contour plotting
screen_xx = linspace(-0.5*screen_width,0.5*screen_width,nw+1)
screen_yy = linspace(-0.5*screen_height,0.5*screen_height,nh+1)
# generate the meshgrid
screen_xx, screen_yy = meshgrid(screen_xx,screen_yy)
# Set up point source coordinates
for i in range (0,nso):
for j in range (0,2):
if j==0: # the x-coordinate is
rso[i,j] = -0.5*length_so + i*length_so/(nso-1)
else: # j = 1 and the y-coordinate is
rso[i,j] = 0.0
# Calculate grid increments
deltaw = screen_width/nw # grid increment, width [cm]
deltah = screen_height/nh # grid increment, height [cm]
# Define the array of screen x and screen y values
for i in range (0,nw+1):
for j in range (0,nh+1):
screen_x[i,j] = -0.5*screen_width + deltaw*i
screen_y[i,j] = -0.5*screen_height + deltah*j
# Calculate the irradiance at each screen point
for k in range (0,nso):
for i in range (0,nw+1):
for j in range (0,nh+1):
# First calculate square of distance from source to screen
rsq[i,j] = (screen_x[i,j]-rso[k,0])**2 + (screen_y[i,j]-rso[k,1])**2 + (zsc-zs)**2
# Calculate x and y coordinates at the aperture
x0[i,j] = rso[k,0] + abs(zs)*(screen_x[i,j] - rso[k,0])/(zsc - zs)
y0[i,j] = rso[k,1] + abs(zs)*(screen_y[i,j] - rso[k,1])/(zsc - zs)
# Check if the ray coordinates at the aperture fall within the apertures
maskx1 = where(absolute(x0) < 0.5*ap1_width,1.0,0.0)
masky1 = where(absolute(y0) < 0.5*ap1_height,1.0,0.0)
maskx2 = where((x0 < ap2_x + 0.5*ap2_width)&(x0 > ap2_x - 0.5*ap2_width),1.0,0.0)
masky2 = where((y0 < ap2_y + 0.5*ap2_height)&(y0 > ap2_y - 0.5*ap2_height),1.0,0.0)
# Calculate the irradiance (note that we are accumulating irradiance)
irradiance = irradiance + (maskx1*masky1 + maskx2*masky2)/rsq
# make a filled contour plot of the period vs overlap and length ratios
figure()
contourf(screen_yy,screen_xx,irradiance,100,cmap=cm.bone)
title('Illumination pattern: $w = $%s cm, $h = $%s; $N = $%d'%(ap1_width,ap1_height,nso))
axis('equal')
xlim(-0.5*screen_width,0.5*screen_width)
ylim(-0.5*screen_height,0.5*screen_height)
xlabel("$x$ [cm]")
ylabel("$y$ [cm]")
colorbar().set_label(label='Irradiance [arb. units]',size=16)
show()