SOC-Barrier.ipynb

Kernel: Python 2 (Ubuntu Linux)

import math
import numpy as np
from mmfutils.performance.fft import fftn, ifftn
from pytimeode import interfaces, mixins, evolvers
import State
from State import State, u

%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

np.array([np.sqrt(8)*126.0, 126.0, 126.0])*u.Hz
u.Hz

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, '+-')

u.micron

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

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))), '-+')

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

%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()
    
        

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()

#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)

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=True)
e.y.plot()

e.y.t=0
with NoInterrupt(ignore=True) as interrupted:
    while e.y.t < 1*u.ms and not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)

b=s.get_Barrier(20)
s.barrierFlag=False
plt.plot(s.xyz[0],s.get_V())

e.y.cooling_phase = 1
e.y.barrierFlag = True
e.y.t = 0;
e.y.barrierOffset = 3.5
e.y.barrierIensity = 10.0
e.y.barrierWidth = np.array([0.3,])
e.y.barrierVelocity = np.array([.1275,])
with NoInterrupt(ignore=True) as interrupted:
    while e.y.t < 2*u.ms and not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        plt.plot(e.y.xyz[0],e.y.get_V())
        display(plt.gcf())
        clear_output(wait=True)