#from pylab import *
from scipy.special import iv
from scipy import optimize
import itertools
from matplotlib import animation
from matplotlib.colors import LogNorm
from .science_basics import *
from numpy import array,sqrt,zeros
[docs]class Morse:
def __init__(self,crystal,h,xi,D):
self.crystal = crystal
self.a = crystal.a1[0]
self.h = h
self.xi = xi
self.D = D
[docs] def calc_morse_potential(self, a, h, xi, D, z_list, gmax):
cc00 = self.calc_coupl_2d_hex(a, h, xi, 0, 0)
gg_range = xrange(-2 * gmax , 2 * gmax + 1)
pot = zeros([len(gg_range), len(gg_range), len(z_list)])
for i, gg in enumerate(itertools.product(gg_range, gg_range)):
gg_x = gg[0]
gg_y = gg[1]
for z_ind, z in enumerate(z_list):
if gg_x == 0 and gg_y == 0:
pot[2 * gmax, 2 * gmax, z_ind] = D * (exp(-2 * xi * z) - 2 * exp(-xi * z))
else:
cc = self.calc_coupl_2d_hex(a, h, xi, gg_x, gg_y, cc00)
pot[gg_x + 2 * gmax, gg_y + 2 * gmax, z_ind] = cc * D * exp(-2 * xi * z)
return pot
[docs] def calc_coupl_2d_hex(self, a, h, xi, G1, G2, div=1.0):
"""
Calculate coupling factor
Args:
a: Edge length of the primitive cell hexagon (in meters)
h: Corrugation parameter of the cosine corrugation function
xi: Stiffnes parameter of the corrugated morse parameter (in m**-1)
G1: First index of G-Vector to the channel
G2: Second index of G-Vector to the channel
div: divison parameter - The solution will be divided by this number
returns:
scalar: Coupling factor from (0,0) to (G1,G2)
"""
kmax = 10
k = arange(-1 * kmax, kmax + 1, 1)
alpha = 2 * xi * h * a / 3.0
qu = iv(k, alpha) * (iv(k + G2, alpha) * iv(k - G1, alpha) + iv(k - G2, alpha) * iv(k + G1, alpha))
return sum(qu) / div
[docs] def show(self, pot):
fig = figure()
a=pot
quad = pcolormesh(a[:, :, 7], norm=LogNorm(vmin=10E-100, vmax=a.max()))
l = len(a[:, 0, 0])
xlim([0, l])
ylim([0, l])
title("Potential")
colorbar()
def animate(i):
title("z-index:{}".format(i))
quad.set_array(a[:, :, i].ravel())
anim = animation.FuncAnimation(fig, animate, frames=len(a[0, 0, :]), blit=False)
show()
[docs] @staticmethod
def boundstate(alpha, D, xi, m):
omg = xi * sqrt(2 * D / m)
gam = 2 * D / (hb * omg)
E_alpha = -D + hb * omg * (alpha + 0.5) * (1 - (alpha + 0.5) / (2 * gam))
return E_alpha
[docs] @staticmethod
def fit(bse,index,m):
# Definiere Parameter
param_0 = array([5.0 * meV, 0.2 / A])
def func(a,b,c):
return Morse.boundstate(a,b,c,m)
popt, pcov = optimize.curve_fit(func, index, bse, p0=param_0)
D,xi = popt
omg = xi * sqrt(2 * D / m)
gam = 2 * D / (hb * omg)
print ('maxbound = ' + str(gam + 0.5))
print ('D = ' + str(D / meV)[0:5] + '+/-' + str(sqrt(pcov[0][0]) / meV)[0:5] + ' meV')
print ('xi = ' + str(xi * A)[0:5] + '+/-' + str(sqrt(pcov[1][1]) * A)[0:5] + ' A**-1')
for i in range(int(gam + 0.5)):
print (Morse.boundstate(i, D, xi,m) / meV)
cont_E = zeros(len(index))
for n in range(len(index)):
cont_E[n] = Morse.boundstate(n, D, xi,m)
print (sqrt(sum((bse / meV - cont_E / meV) ** 2)) / len(index))