phaseMatchingDemo.py

# -*- coding: utf-8 -*-
"""
Created on Tue Sep 10 19:32:39 2013
 
@author: joey
"""
 
#this program will create a two-dimensional demonstration of the phase matching condition
import time
 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
 
startTime = time.time()
 
radius_spotsize = 5.0
 
numberOfEmitters = 100
 
k_1 = np.array([1.0, 0.0]) 
e_1 = np.array([0.0, 1.0]) 
 
k_2 = np.array([0.0, 1.0]) 
e_2 = np.array([1.0, 0.0]) 
 
k_output = k_2 + k_1
AbsK_output = np.sqrt(np.dot(k_output,k_output))
 
transitionDiploeMagniture = .5
 
spaceQuantizationPointsPerDimension = 200
 
xMax = 20
xMin = -xMax
 
yMax = 20
yMin = -yMax
 
xValues = np.linspace(xMin, xMax, spaceQuantizationPointsPerDimension)
yValues = np.linspace(yMin, yMax, spaceQuantizationPointsPerDimension)
XVAL, YVAL = np.meshgrid(xValues, yValues)
 
 
def emptyLightEmissionGrid():
    return np.zeros((spaceQuantizationPointsPerDimension, spaceQuantizationPointsPerDimension), dtype=np.complex)
 
def randomUnitVector():
    theta = np.random.rand() * 2.0 * np.pi
    x = np.cos(theta)
    y = np.sin(theta)
    return np.array([x, y])
     
def newEmitterCoordinate():
    r = np.random.rand() * radius_spotsize
    return r * randomUnitVector()
 
def newEmitterTransitionDipole():
    r = transitionDiploeMagniture
    return r * randomUnitVector()
     
def newEmitterEmissionVectorFunction(amplitude, waveVectorNorm, phase):
    return lambda x, y: amplitude * np.exp(1.0j * (waveVectorNorm * np.sqrt(x**2 + y**2) + phase) ) 
#    return lambda x, y: amplitude * np.cos( (waveVectorNorm * np.sqrt(x**2 + y**2) + phase) )
     
def plotWavevector(k):
    k = k * xMax / 2.0
    plt.plot([0, k[0]], [0, k[1]], "Black", linewidth=5)
 
#totalXAmplitude = emptyLightEmissionGrid()
#totalYAmplitude = emptyLightEmissionGrid()
totalZAmplitude = emptyLightEmissionGrid()
 
emitterXValues = []
emitterYValues = []
 
cumulativeIntensities = []
 
for emitterIndex in range(numberOfEmitters):
     
    newCoordinates = newEmitterCoordinate()
    emitterXValues.append(newCoordinates[0])
    emitterYValues.append(newCoordinates[1])
     
    newDipole = newEmitterTransitionDipole()
     
    newPhase = np.dot(k_output, newCoordinates)
     
    emissionAmplitude = np.dot(e_1, newDipole) * np.dot(e_2, newDipole)
     
    newAmplitudeFunction = newEmitterEmissionVectorFunction(emissionAmplitude, AbsK_output, newPhase)
     
    newOscillations = newAmplitudeFunction(XVAL- newCoordinates[0], YVAL- newCoordinates[1])
     
    totalZAmplitude = totalZAmplitude + newOscillations / emissionAmplitude
     
    cumulativeIntensities.append(np.abs(totalZAmplitude)**2 / np.max(np.abs(totalZAmplitude)**2))
 
intensity = np.abs(totalZAmplitude)**2
scaledIntensity = intensity / np.max(intensity)
plt.figure()
 
plt.plot(emitterXValues, emitterYValues, "r.")
plt.contourf(XVAL, YVAL, scaledIntensity, 100)
plt.title("Scaled Intensity")
plotWavevector(k_1)
plotWavevector(k_2)
plotWavevector(k_output)
plt.colorbar()
 
 
#plt.figure()
#plt.plot(emitterXValues, emitterYValues, "r.")
#plt.contourf(XVAL, YVAL, np.real(totalZAmplitude), 100)
#plt.title("Real")
#plotWavevector(k_1)
#plotWavevector(k_2)
#plotWavevector(k_output)
#plt.colorbar()
 
 
print "Time Elapsed",  time.time() - startTime
 
#plt.figure()
     
fig = plt.figure()
im = plt.contourf(XVAL, YVAL, cumulativeIntensities[0], 100)
plt.title("Scaled Intensity, n = 1")
plotWavevector(k_1)
plotWavevector(k_2)
plotWavevector(k_output)
plt.plot(emitterXValues[0], emitterYValues[0], "r.")
plt.colorbar()
 
plt.legend()
ax = fig.gca()
 
def animate(i, data,  ax, fig):
    fig.clf()
    try:
        im = plt.contourf(XVAL, YVAL, cumulativeIntensities[i], 100)
    except:
        return None
    plt.title("Scaled Intensity, n = %s" % str(i))
    plotWavevector(k_1)
    plotWavevector(k_2)
    plotWavevector(k_output)
    plt.plot(emitterXValues[0:i], emitterYValues[0:i], "r.")
    plt.colorbar()
    return im,
 
#anim = animation.FuncAnimation(fig, animate, frames = numberOfEmitters, fargs=(cumulativeIntensities, ax, fig))
#anim.save("phaseMatchTest.mp4", bitrate = 500,  fps=5)