In [71]:
import math
import numpy as np
from mmfutils.performance.fft import fftn, ifftn
from pytimeode import interfaces, mixins, evolvers
import States
from States import State, u
In [72]:
%pylab inline --no-import-all
from IPython.display import display, clear_output
s = State(Nxyz=(128,), Lxyz=(40*u.micron,),ws = np.array([26.0, ])*u.Hz, N=1e5)
s.plot()
s.ws
Out[72]:
In [73]:
np.array([np.sqrt(8)*126.0, 126.0, 126.0])*u.Hz
u.Hz
Out[73]:
In [74]:
def get_err(N=128, L=24*u.micron):
s = State(Nxyz=(N,), Lxyz=(L,), N=1e5)
s.g = 0
a = np.sqrt(u.hbar/u.m/s.ws[0])
x = s.xyz[0]
psi_0 = np.exp(-(x/a)**2/2.0)
s[...] = psi_0
s.normalize()
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=True)
return abs(dy[...]).max()
Ns = np.arange(2,128,2)
#Ns = 2**np.arange(2,8)
errs = map(get_err, Ns)
plt.semilogy(Ns, errs, '+-')
Out[74]:
In [75]:
u.micron
Out[75]:
In [76]:
s = State(Nxyz=(46,), Lxyz=(23*u.micron,))
a = np.sqrt(u.hbar/u.m/s.ws[0])
L, N = s.Lxyz[0], s.Nxyz[0]
k_max = np.pi*(N-2)/L # For Khalid...
print (k_max, s.kxyz[0].max())
print(np.exp(-(L/2/a)**2/2)) # Wavefunction drops by factor of macheps
In [77]:
psi_0 = s.xyz[0]*np.exp(-(s.xyz[0]/a)**2/2)
plt.semilogy(np.fft.fftshift(s.kxyz[0]),
np.fft.fftshift(abs(np.fft.fft(psi_0))), '-+')
Out[77]:
I got error for this part, the error message is shown in the output, That comfuses me for a while.
Exact Solution with Interactions¶
I am trying to figure out which part is the interactions term. The formula in the code looks strange to me, I know there is an interacton term $gn_0$, and it's added to the homarmnic poterntial term. But it also includes a $2k$ term. $V_{ext} = \frac{\hbar^2}{2m} [4(kx)^2 - 2k) - gn_0]$
OK, I remember what we did in the meeting¶
In [78]:
%pylab inline --no-import-all
from IPython.display import display, clear_output
import States
from States import State, u
s = State(Nxyz=(64,), Lxyz=(23*u.micron,), N=1e5)
a = np.sqrt(u.hbar/u.m/s.ws[0])
x = s.xyz[0]
psi_0 = np.exp(-(x/a)**2/2)
#print(x)
#print(psi_0)
class State1(State):
def __init__(self, *v, **kw):
State.__init__(self, *v, **kw)
a = np.sqrt(u.hbar/u.m/self.ws[0])
x = self.xyz[0]
k = 1./2./a**2
psi_0 = 4.0*np.exp(-(x/a)**2/2)
n_0 = abs(psi_0)**2
self._V_ext = (u.hbar**2/2.0/u.m*(4*(k*x)**2 - 2*k) - self.g*n_0)
self.data[...] = psi_0
self.get_Vext = lambda: self._V_ext
self.pre_evolve_hook()
def get_Barrier(self,t_):
"""Return the time dependece barriry"""
return sum(np.exp(_x - _v * t_)**2
for _v, _x in zip(self.v, self.xyz))
s = State1(Nxyz=(64,), Lxyz=(23*u.micron,))
s.plot()
plt.plot(x, s.get_Vext())
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=False)
abs(dy[...]).max()
Out[78]:
In [79]:
#print(s.kxyz)
K1=sum((u.hbar*_k)**2/2.0/u.m for _k in s.kxyz)
#print(K1)
K=(u.hbar*s.kxyz[0])**2/u.m
#print(K)
In [80]:
from mmfutils.contexts import NoInterrupt
from pytimeode.evolvers import EvolverSplit, EvolverABM
from IPython.display import display, clear_output
s = State1(Nxyz=(64*4,), Lxyz=(23*u.micron,))
assert np.allclose(s._N, s.get_N())
s[...] = 1.0
s.normalize()
s.cooling_phase = 1j
E_max = u.hbar**2*np.abs(s.kxyz).max()**2/2.0/u.m
#e = EvolverSplit(s, dt=0.01*u.hbar/E_max, normalize=True)
e = EvolverABM(s, dt=0.1*u.hbar/E_max, normalize=False)
e.y.plot()
with NoInterrupt(ignore=True) as interrupted: while e.y.t < 4*u.ms and not interrupted: e.evolve(100) plt.clf() e.y.plot() display(plt.gcf()) clear_output(wait=True)
In [81]:
s = State1(Nxyz=(64*4,), Lxyz=(23*u.micron,))
s.get_Barrier(1)
In [ ]: