#

# Shadows_Exercise_1.py

#

# A point source is located a specified distance from

# a rectangular aperture that 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 the 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 original code by ERB

# 20160609 clean up by ERB

#

# import the commands needed to make the plot

from pylab import xlabel,ylabel,show,contourf,axis,colorbar,figure,title,xlim,ylim

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]

rso = [5.0,0.0] # (x,y) location of the source [cm,cm]

ap_width = 4.0 # aperture width [cm]

ap_height = 3.0 # aperture height [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

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)) # source-to-screen point distance, squared

x0 = zeros((nw+1,nh+1)) # x coordinates of ray intercepts at aperture plane

y0 = zeros((nw+1,nh+1)) # y coordinates of ray intercepts at aperture plane

# initialize values

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)

# create the meshgrid on which to evaluate the irradiance

screen_xx, screen_yy = meshgrid(screen_xx,screen_yy)

# 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 intensity at each screen point

for i in range (0,nw+1):

for j in range (0,nh+1):

# First calculate the square of the distance from source to screen

rsq = (screen_x[i,j]-rso[0])**2 + (screen_y[i,j]-rso[1])**2 + (zsc - zs)**2

# Calculate x and y coordinates at the aperture

x0[i,j] = rso[0] + abs(zs)*(screen_x[i,j] - rso[0])/(zsc - zs)

y0[i,j] = rso[1] + abs(zs)*(screen_y[i,j] - rso[1])/(zsc - zs)

# Check if the coordinates fall within the aperture

maskx = where(absolute(x0) < 0.5*ap_width,1.0,0.0)

masky = where(absolute(y0) < 0.5*ap_height,1.0,0.0)

# Calculate the intensity

irradiance = maskx*masky/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)

axis('equal')

xlim(-0.5*screen_width,0.5*screen_width)

ylim(-0.5*screen_height,0.5*screen_height)

xlabel("$x$ [cm]",size = 16)

ylabel("$y$ [cm]",size = 16)

title('Illumination pattern: $w = $%s cm, $h = $%s'%(ap_width,ap_height),size = 16)

colorbar().set_label(label='Irradiance [arb. units]',size=16)

show()