Add ps14

modern-hw7-2d.py

@@ -0,0 +1,17 @@
+import numpy as np
+import pylab as pl
+
+C = np.sqrt(12)
+x_vals = np.linspace(0, 1, 1000)
+psi_sqrd = C ** 2 * (x_vals ** 2 - x_vals ** 3)
+
+pl.rcParams['figure.dpi'] = 300
+pl.plot(x_vals, psi_sqrd)
+pl.title("|Psi|^2")
+pl.xlabel("x")
+pl.ylabel("|Psi|^2")
+pl.vlines(0.6, ymin=0, ymax=1.8, label="<x>", color="red")
+pl.vlines(2 / 3, ymin=0, ymax=1.8, label="x_mp", color="green")
+pl.legend()
+pl.grid()
+pl.show()

ps14-1.py

@@ -0,0 +1,74 @@
+# P372-PS14-FFT
+# Matthew Oros
+"""
+This script introduces the Real FFT and Inverse Real FFT capabilities.
+"""
+
+import numpy as np
+import pylab as pl
+import numpy.fft as fft
+
+# Some constants, I left these all the same
+numPoints = 1000  # number of points in x-waveform, how much resolution
+numFreqs = numPoints / 2 + 1  # (numPoints/2 -1) complex & 2 real k-amplitudes
+# This is still same numPoints of bits of information!
+# The wavenumber k's are also called "spatial frequencies"
+# We are going to look at the region between 0 and 2 in x-space,
+# but suppose we don't want anything to repeat within twice this, so
+windowLength = 4.0  # effective window length
+dx = windowLength / numPoints  # Step size in x-space
+x = np.arange(0.0, windowLength, dx)  # All the x-positions
+
+# Now prepare the spectrum and k wavenumbers
+dk = 2 * np.pi / windowLength  # Step size in k-space
+ks = np.arange(numFreqs) * dk  # All the k's, highest k found at ks[-1]
+
+# HERE IS WHERE THE REAL WORK GETS DONE, the above is all just setting up
+# Choose the <k>
+k0 = 0.1 * ks[-1]  # I chose 0.1 to have lower frequency inside the wave packet
+alpha = 0.1  # I chose 0.1 to make the wave packet wider
+# Now define the k-spectrum amplitudes
+phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)  # the spectral amplitudes (complex)
+# And generate the wavefunction from phi(k) using Inverse Real FFT
+psi_x = fft.irfft(phi_k)  # inverse real FFT gets the real wavefunction
+
+# I will now normalize the wave function,
+# which was not done in the original script
+norm_factor = 1 / np.sqrt(np.sum(psi_x * psi_x))
+psi_x *= norm_factor
+
+pl.figure(1)
+pl.plot(ks, phi_k)  # phi_k could be complex, which would throw a warning here
+# if so, could use: pl.plot(ks, np.abs(phi_k) )
+pl.title(r"Spectrum $\phi(k)$")
+pl.xlabel("wavenumber k")
+pl.ylabel("spectral amplitude")
+
+pl.figure(2)
+pl.plot(x, psi_x)
+pl.title(r"Wavefunction $\psi(x)$")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\psi(x)$")
+pl.annotate("This piece also appears to left of x=0!", (3.8, -0.010),
+            xytext=(1.0, -0.03), arrowprops=dict(facecolor='black', shrink=0.05))
+
+# Here's a fancier figure, without recalculating anything but noticing
+# that the wavepacket right now is centered around x=0.  Some of that is
+# showing up on the right side of the window, because the whole signal is
+# repeating with a period of 1 windowLength, so we will copy that and stitch
+# it on the left side in a different color.  (If you look closely, you will
+# see that the last point of the "orange" needs to be connected to the first
+# point of the "blue")
+# In the NEXT script P372-PS14-FFT-shift.py, the wavepacket is shifted in x-space
+pl.figure(3)
+cutPoint = (3 * numPoints) // 4  # 3/4 of the way through the window
+pl.plot(x[0:cutPoint], psi_x[0:cutPoint])  # as before, but only to cutPoint
+pl.plot(x[cutPoint:] - windowLength, psi_x[cutPoint:])  # put last 1/4 on left side
+pl.plot([-dx, 0], [psi_x[-1], psi_x[0]], 'k')  # Final stitch
+pl.title(r"Wavefunction $\psi(x)$")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\psi(x)$")
+
+pl.show()  # show the graph windows and loop here until they are closed
+
+print("Done.")

ps14-2.py

@@ -0,0 +1,79 @@
+# P372-PS14-FFT
+# Matthew Oros
+"""
+This script introduces the Real FFT and Inverse Real FFT capabilities.
+"""
+
+import numpy as np
+import pylab as pl
+import numpy.fft as fft
+
+# Some constants, I left these all the same
+numPoints = 1000  # number of points in x-waveform, how much resolution
+numFreqs = numPoints / 2 + 1  # (numPoints/2 -1) complex & 2 real k-amplitudes
+# This is still same numPoints of bits of information!
+# The wavenumber k's are also called "spatial frequencies"
+# We are going to look at the region between 0 and 2 in x-space,
+# but suppose we don't want anything to repeat within twice this, so
+windowLength = 4.0  # effective window length
+dx = windowLength / numPoints  # Step size in x-space
+x = np.arange(0.0, windowLength, dx)  # All the x-positions
+
+# Now prepare the spectrum and k wavenumbers
+dk = 2 * np.pi / windowLength  # Step size in k-space
+ks = np.arange(numFreqs) * dk  # All the k's, highest k found at ks[-1]
+
+# HERE IS WHERE THE REAL WORK GETS DONE, the above is all just setting up
+# Choose the <k>
+k0 = 0.1 * ks[-1]  # I chose 0.1 to have lower frequency inside the wave packet
+alpha = 20  # 0.1  # I chose 0.1 to make the wave packet wider
+# Now define the k-spectrum amplitudes
+# phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)  # the spectral amplitudes (complex)
+# Upside-down parabola!
+phi_k = np.piecewise(ks,
+                     [(ks - k0 < -alpha), (-alpha <= ks - k0) & (ks - k0 <= alpha), (ks - k0 > alpha)],
+                     [0, lambda k: (1 / alpha ** 2) * (k - k0) ** 2, 0])
+
+# And generate the wavefunction from phi(k) using Inverse Real FFT
+psi_x = fft.irfft(phi_k)  # inverse real FFT gets the real wavefunction
+
+# I will now normalize the wave function,
+# which was not done in the original script
+norm_factor = 1 / np.sqrt(np.sum(psi_x * psi_x))
+psi_x *= norm_factor
+
+pl.figure(1)
+pl.plot(ks, phi_k)  # phi_k could be complex, which would throw a warning here
+# if so, could use: pl.plot(ks, np.abs(phi_k) )
+pl.title(r"Spectrum $\phi(k)$")
+pl.xlabel("wavenumber k")
+pl.ylabel("spectral amplitude")
+
+pl.figure(2)
+pl.plot(x, psi_x)
+pl.title(r"Wavefunction $\psi(x)$")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\psi(x)$")
+pl.annotate("This piece also appears to left of x=0!", (3.8, -0.010),
+            xytext=(1.0, -0.03), arrowprops=dict(facecolor='black', shrink=0.05))
+
+# Here's a fancier figure, without recalculating anything but noticing
+# that the wavepacket right now is centered around x=0.  Some of that is
+# showing up on the right side of the window, because the whole signal is
+# repeating with a period of 1 windowLength, so we will copy that and stitch
+# it on the left side in a different color.  (If you look closely, you will
+# see that the last point of the "orange" needs to be connected to the first
+# point of the "blue")
+# In the NEXT script P372-PS14-FFT-shift.py, the wavepacket is shifted in x-space
+pl.figure(3)
+cutPoint = (3 * numPoints) // 4  # 3/4 of the way through the window
+pl.plot(x[0:cutPoint], psi_x[0:cutPoint])  # as before, but only to cutPoint
+pl.plot(x[cutPoint:] - windowLength, psi_x[cutPoint:])  # put last 1/4 on left side
+pl.plot([-dx, 0], [psi_x[-1], psi_x[0]], 'k')  # Final stitch
+pl.title(r"Wavefunction $\psi(x)$")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\psi(x)$")
+
+pl.show()  # show the graph windows and loop here until they are closed
+
+print("Done.")

ps14-3-a.py

@@ -0,0 +1,68 @@
+# P372-PS14-FFT-shift.py
+"""
+A wavepacket is made to move.  Based on P372-PS14-FFT.py
+"""
+
+import numpy as np
+import pylab as pl
+import numpy.fft as fft
+
+# Some basic numbers, first.  These are all design decisions you can change!
+numPoints = 1000  # number of points in x-waveform, how much resolution
+numFreqs = numPoints / 2 + 1  # (numPoints/2 -1) complex & 2 real k-amplitudes
+# This is still same numPoints of bits of information!
+# The wavenumber k's are also called "spatial frequencies"
+# We are going to look at the region between 0 and 2 in x-space,
+# but suppose we don't want anything to repeat within twice this, so
+windowLength = 4.0  # effective window length
+dx = windowLength / numPoints  # Step size in x-space
+x = np.arange(0.0, windowLength, dx)  # All the x-positions
+
+# Now prepare the spectrum and k wavenumbers
+dk = 2 * np.pi / windowLength  # Step size in k-space
+ks = np.arange(numFreqs) * dk  # All the k's, highest k found at ks[-1]
+
+# HERE IS WHERE THE REAL WORK GETS DONE, the above is all just setting up
+# Choose the <k>
+k0 = 0.1 * ks[-1]  # 25% of highest k value <-- PICK THIS!
+alpha = 0.1  # smaller alpha means wider phi_k, narrower psi_x <-- AND THIS!
+# Now define the k-spectrum amplitudes
+phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)  # the spectral amplitudes (complex)
+# And generate the wavefunction from phi(k) using Inverse Real FFT
+psi_x = fft.irfft(phi_k)  # inverse real FFT gets the real wavefunction
+
+# Now we want the wavepacket to move!  Let's define a velocity:
+v = 0.5
+# First let's have it move just 1/4 cycle to the right
+# Thus k0 v dt = pi/2
+dt = 0.5 * np.pi / v / k0
+# Calculate the phase shifts
+phaseShift = np.exp(-1j * ks * v * dt)
+# And apply them before doing the F.T.
+psi_x_shift = fft.irfft(phi_k * phaseShift)
+
+# And now let's get them to shift all the way to x=1.0, which means vt = 1.0
+t1 = 1.0 / v
+# Calculate the phase shifts
+phaseShift1 = np.exp(-1j * ks * v * t1)
+# And apply them before doing the F.T.
+psi_x_shift1 = fft.irfft(phi_k * phaseShift1)
+
+pl.figure(1)
+pl.plot(ks, phi_k)  # phi_k could be complex, which would throw a warning here
+# if so, could use: pl.plot(ks, np.abs(phi_k) )
+pl.title(r"Spectrum $\phi(k)$")
+pl.xlabel("wavenumber k")
+pl.ylabel("spectral amplitude")
+
+pl.figure(2)
+pl.plot(x[:numPoints // 2], psi_x[:numPoints // 2])  # only plot half of the points
+pl.plot(x[:numPoints // 2], psi_x_shift[:numPoints // 2])  # These are the shifted wavepacket
+pl.plot(x[:numPoints // 2], psi_x_shift1[:numPoints // 2])  # These are the shifted wavepacket
+pl.title(r"Wavefunction $\psi(x)$ shifted")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\psi(x)$")
+
+pl.show()  # show the graph windows and loop here until they are closed
+
+print("Done.")

ps14-3-b.py

@@ -0,0 +1,70 @@
+# P372-PS14-FFT-animate.py
+"""
+A wavepacket is made to move.  Based on P372-PS14-FFT.py
+and uses matplotlib.animation library
+"""
+
+import numpy as np
+import pylab as pl
+import matplotlib.animation as animation
+import numpy.fft as fft
+
+# Some basic numbers, first.  These are all design decisions you can change!
+numPoints = 1000  # number of points in x-waveform, how much resolution
+numFreqs = numPoints / 2 + 1  # (numPoints/2 -1) complex & 2 real k-amplitudes
+# This is still same numPoints of bits of information!
+# The wavenumber k's are also called "spatial frequencies"
+# We are going to look at the region between 0 and 2 in x-space,
+# but suppose we don't want anything to repeat within twice this, so
+windowLength = 4.0  # effective window length
+dx = windowLength / numPoints  # Step size in x-space
+x = np.arange(0.0, windowLength, dx)  # All the x-positions
+
+# Now prepare the spectrum and k wavenumbers
+dk = 2 * np.pi / windowLength  # Step size in k-space
+ks = np.arange(numFreqs) * dk  # All the k's, highest k found at ks[-1]
+
+# HERE IS WHERE THE REAL WORK GETS DONE, the above is all just setting up
+# Choose the <k>
+k0 = 0.1 * ks[-1]  # 25% of highest k value <-- PICK THIS!
+alpha = 0.1  # smaller alpha means wider phi_k, narrower psi_x <-- AND THIS!
+# Now define the k-spectrum amplitudes
+phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)  # the spectral amplitudes (complex)
+# And generate the wavefunction from phi(k) using Inverse Real FFT
+psi_x = fft.irfft(phi_k)  # inverse real FFT gets the real wavefunction
+
+# Now we want the wavepacket to move!  Let's define a velocity:
+v = 1.0
+
+fig, ax = pl.subplots()  # Create the graphics canvas
+pl.title(r"Wavefunction $\Psi(x,t)$ animated")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\Psi(x,t)$")
+
+T = 1.0  # total time duration
+Nsteps = 200  # number of steps of animation
+times = np.linspace(0, T, Nsteps + 1)  # the times at which to calculate the shifts
+# Draw the first plot and get the lines object
+lines = ax.plot(x[:numPoints // 2], psi_x[:numPoints // 2], lw=1)[0]
+
+
+def update(frame):
+    # This calculates the new values for the plot, called by the animation
+    time_now = times[frame]
+    # Velocity changes with time!
+    v = np.cos(time_now * 10)
+    # Calculate the phase shifts
+    phase_shift = np.exp(-1j * ks * v * time_now)
+    # And apply them before doing the F.T.
+    psi_x_shift = fft.irfft(phi_k * phase_shift)
+    lines.set_ydata(psi_x_shift[:numPoints // 2])
+    return lines
+
+
+# Connect the update function to the animation
+ani = animation.FuncAnimation(fig=fig, func=update, frames=Nsteps, interval=30)
+
+# Show the animation and run it until the window is closed
+pl.show()
+
+print("Done.")

ps14-3-c.py

@@ -0,0 +1,70 @@
+# P372-PS14-FFT-animate.py
+"""
+A wavepacket is made to move.  Based on P372-PS14-FFT.py
+and uses matplotlib.animation library
+"""
+
+import numpy as np
+import pylab as pl
+import matplotlib.animation as animation
+import numpy.fft as fft
+
+# Some basic numbers, first.  These are all design decisions you can change!
+numPoints = 1000  # number of points in x-waveform, how much resolution
+numFreqs = numPoints / 2 + 1  # (numPoints/2 -1) complex & 2 real k-amplitudes
+# This is still same numPoints of bits of information!
+# The wavenumber k's are also called "spatial frequencies"
+# We are going to look at the region between 0 and 2 in x-space,
+# but suppose we don't want anything to repeat within twice this, so
+windowLength = 4.0  # effective window length
+dx = windowLength / numPoints  # Step size in x-space
+x = np.arange(0.0, windowLength, dx)  # All the x-positions
+
+# Now prepare the spectrum and k wavenumbers
+dk = 2 * np.pi / windowLength  # Step size in k-space
+ks = np.arange(numFreqs) * dk  # All the k's, highest k found at ks[-1]
+
+# HERE IS WHERE THE REAL WORK GETS DONE, the above is all just setting up
+# Choose the <k>
+k0 = 0.1 * ks[-1]  # 25% of highest k value <-- PICK THIS!
+alpha = 0.1  # smaller alpha means wider phi_k, narrower psi_x <-- AND THIS!
+# Now define the k-spectrum amplitudes
+phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)  # the spectral amplitudes (complex)
+# And generate the wavefunction from phi(k) using Inverse Real FFT
+psi_x = fft.irfft(phi_k)  # inverse real FFT gets the real wavefunction
+
+# Now we want the wavepacket to move!  Let's define a velocity:
+v = 1.0
+
+fig, ax = pl.subplots()  # Create the graphics canvas
+pl.title(r"Wavefunction $\Psi(x,t)$ animated")
+pl.xlabel("position x")
+pl.ylabel(r"wavefunction $\Psi(x,t)$")
+
+T = 1.0  # total time duration
+Nsteps = 200  # number of steps of animation
+times = np.linspace(0, T, Nsteps + 1)  # the times at which to calculate the shifts
+# Draw the first plot and get the lines object
+lines = ax.plot(x[:numPoints // 2], psi_x[:numPoints // 2], lw=1)[0]
+
+
+def update(frame):
+    # This calculates the new values for the plot, called by the animation
+    time_now = times[frame]
+    # Calculate the phase shifts
+    phase_shift = np.exp(-1j * ks * 2 * v * time_now)
+    alpha = time_now
+    phi_k = np.exp(-0.5 * alpha * (ks - k0) ** 2)
+    # And apply them before doing the F.T.
+    psi_x_shift = fft.irfft(phi_k * phase_shift)
+    lines.set_ydata(psi_x_shift[:numPoints // 2])
+    return lines
+
+
+# Connect the update function to the animation
+ani = animation.FuncAnimation(fig=fig, func=update, frames=Nsteps, interval=30)
+
+# Show the animation and run it until the window is closed
+pl.show()
+
+print("Done.")