User:Thierry Dugnolle/Python/Interference of a particle with itself

print("Single particle interference")

from PIL import Image, ImageDraw
import random
import math
import cmath

imWidth = 800
imHeight = 440
lengthUnit = 96.0 # length unit in number of pixels
a = 1.0  # packet width
m = 100.0  # packet mass
x1Start = -5.0*a
x2Start = x1Start
y1Start = -2.5*a
y2Start = -y1Start
kX1 = 20.0
kY1 = kX1*y1Start/x1Start
kX2 = kX1
kY2 = -kY1
dt = 0.16
imNumber = 310

WnumPixelCenter = imWidth/2.0 - 0.5
HnumPixelCenter = imHeight/2.0 - 0.5

im = Image.new("RGB", (imWidth, imHeight), (0, 0, 0))
px = im.load()
def psi(k0, x, t):
    r1 = pow(pow(a, 4) + 4 * t * t / (m * m), 0.25)
    theta = 0.5 * math.atan(2 * t / (m * a * a))
    phi = -theta - k0 * k0 * t / (2.0 * m)
    c1 = pow(math.e, 1j * (phi + k0 * x) - pow( x - k0 * t / m, 2)/ (a*a+1j*2.0*t/m))
    return pow(2*a*a/math.pi,0.25) * c1 / r1

for imNum in range(imNumber):
    t = imNum * dt
    for i in range(imWidth):
        for j in range(imHeight):
            x = (i - WnumPixelCenter) / lengthUnit
            y = (HnumPixelCenter - j) / lengthUnit
            psiX1 = psi(kX1, x-x1Start, t)
            psiY1 = psi(kY1, y-y1Start, t)
            psiX2 = psi(kX2, x - x2Start, t)
            psiY2 = psi(kY2, y - y2Start, t)
            psiXY = psiX1 * psiY1 + psiX2 * psiY2
            psiPhase = cmath.phase(psiXY)
            psiPhaseNorm= (psiPhase+math.pi)/(2.0*math.pi)
            psiMod = abs(psiXY)
            bright= psiMod*psiMod
            if psiPhaseNorm<1.0/3.0:
                red=1.0-3.0*psiPhaseNorm
                blue=0.0
            else:
                if psiPhaseNorm<2.0/3.0:
                    red=0.0
                    blue=3.0*(psiPhaseNorm-1.0/3.0)
                else:
                    red=3.0*(psiPhaseNorm-2.0/3.0)
                    blue=1.0-3.0*(psiPhaseNorm-2.0/3.0)
            px[i,j] = (math.floor(255.0*red*bright), math.floor(255.0*bright), math.floor(255.0*blue*bright))

    im.save("SPI" + str(1000 + imNum) + ".bmp")

print("Good bye")