In [2]:
import numpy as np
import matplotlib.pyplot as plt
from Loader2D import *
%matplotlib qt5
In [3]:
filepath = "/Users/max/git/Ahti/cmake-build-default/bin/init.h5"
nt = 0
nx = 500
ny = 500
xmax = 6.0
xmin = -6.0
ymax = xmax
ymin = xmin
w1 = 2.0 * np.pi / 1000.0
a1 = 0.01
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
x_mesh, y_mesh = np.meshgrid(x,y)
dt = 0.001
In [4]:
class W_2D:
"""
This class is an implementation of a
time only disturbance in order to further
"""
def __init__(self,a1,a2,w1,w2):
self.a1 = a1
self.a2 = a2
self.w1 = w1
self.w2 = w2
def __call__(self,t, x=0, y=0):
return self.a1 * np.cos(self.w1 * t) + self.a2 * np.cos(self.w2 * t)
class VHarm_2D:
"""
This class implements the harmonic potential
in two dimensions.
"""
def __init__(self):
pass
def __call__(self, x, y):
return 0.5 * (x**2 + y**2)
In [5]:
class U_lin_dist:
"""
Implementation of the time evolution
operator using the above method
"""
def __init__(self, E0, W, dt):
self.E0 = E0
self.W = W
self.dt = dt
def __call__(self, t):
w_val = self.W(t)
return np.exp(-1j * (self.E0 + np.trapz(w_val, dx=self.dt)))
In [6]:
%%latex
$| \, \psi(t_0) \, \rangle$
In [7]:
Psi0 = lambda x: np.pi**(-0.25) * np.exp(-x**2*0.5)
Psi2D = lambda x,y: Psi0(x) * Psi0(y)
In [8]:
psi_an = Psi2D(x_mesh, y_mesh) + 0 * 1j
w = W_2D(a1, a1, w1, w1)
U = U_lin_dist(1.0, w, dt)
t = np.linspace(0, nt*dt,nt)
psi_an *= np.exp(- 1j * 1.0 * nt * dt)
In [11]:
file = CLoader2D(nx, ny, filepath)
psi_num = file.get_complex_data("/real","/imag")
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel(r"x $(a.u.)$")
ax.set_ylabel(r"y $(a.u.)$")
ax.set_zlabel(r"$|\psi|^2$")
surf = ax.plot_surface(x_mesh,
y_mesh,
np.abs(psi_an.real - psi_num.real),
cmap=cm.magma_r,
linewidth=1,
antialiased=False)
print("max. rel. error real part {}".format(np.max(np.abs(psi_an.real - psi_num.real))/np.max(np.abs(psi_an.real))))
#print("max. rel. error imaginary part {}".format(np.max(np.abs(psi_an.imag - psi_num.imag))/np.max(np.abs(psi_an.imag))))
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
In [ ]:
In [ ]: