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

%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

Populating the interactive namespace from numpy and matplotlib

array([ 0.0004094])

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

1.5746097513521148e-05

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

[<matplotlib.lines.Line2D at 0x7f84ceeb4c10>]

u.micron

1.0

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

(6.010003337302212, 6.010003337302213)
9.87009441159e-15

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

[<matplotlib.lines.Line2D at 0x7f84cf2873d0>]

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

Populating the interactive namespace from numpy and matplotlib

1.2001072726244708e-14

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

s = State1(Nxyz=(64*4,), Lxyz=(23*u.micron,))
s.get_Barrier(1)


AttributeErrorTraceback (most recent call last)
<ipython-input-81-0e8f6ad2b40a> in <module>()
      1 s = State1(Nxyz=(64*4,), Lxyz=(23*u.micron,))
----> 2 s.get_Barrier(1)

<ipython-input-78-660e4ab3070a> in get_Barrier(self, t_)
     28         """Return the time dependece barriry"""
     29         return sum(np.exp(_x - _v * t_)**2
---> 30                              for _v, _x in zip(self.v, self.xyz))
     31 
     32 s = State1(Nxyz=(64,), Lxyz=(23*u.micron,))

/projects/990b1d3c-5a22-40a3-849e-cd68d8e43eb4/.local/lib/python2.7/site-packages/pytimeode/mixins.pyc in __getattr__(self, name)
    118                 "Cannot get attribute `{}`.  Did you mean `writeable`?"
    119                 .format(name))
--> 120         self.__getattribute__(name)
    121 
    122     def __setattr__(self, name, value):

AttributeError: 'State1' object has no attribute 'v'