Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 6 13:17:28 2022
@author: stefan
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import codecs
def energy(kappa,alpha,Q,B) :
G=kappa*np.array([[np.cos(alpha)**2],[np.sin(alpha)**2],[np.sqrt(2)*np.cos(alpha)*np.sin(alpha)]])-B
return np.matmul(np.transpose(G),np.matmul(Q,G))[0,0]
def xytokappaalpha(x,y):
if y>0:
return [np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
else:
return [-np.sqrt(x**2+y**2), np.abs(np.arctan2(y,x))]
# Read effective quantites
def ReadEffectiveQuantities(QFilePath, BFilePath):
# Read Output Matrices (effective quantities)
# From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
# -- Read Matrix Qhom
X = []
# with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
with codecs.open(QFilePath, encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
Q = np.array([[X[0][2], X[1][2], X[2][2]],
[X[3][2], X[4][2], X[5][2]],
[X[6][2], X[7][2], X[8][2]] ])
# -- Read Beff (as Vector)
X = []
# with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
with codecs.open(BFilePath, encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
B = np.array([X[0][2], X[1][2], X[2][2]])
return Q, B
number=7
kappa=np.zeros(number)
for n in range(0,number):
# Read from Date
print(str(n))
DataPath = './experiment/wood-bilayer/results/'+str(n)
QFilePath = DataPath + '/QMatrix.txt'
BFilePath = DataPath + '/BMatrix.txt'
Q, B = ReadEffectiveQuantities(QFilePath,BFilePath)
# Q=0.5*(np.transpose(Q)+Q) # symmetrize
B=np.transpose([B])
#
r, theta = np.meshgrid(np.linspace(0,length,N),np.radians(np.linspace(0, 360, N)))
E=np.zeros(np.shape(r))
for i in range(0,N):
for j in range(0,N):
if theta[i,j]<np.pi:
E[i,j]=energy(r[i,j],theta[i,j],Q,B)
else:
E[i,j]=energy(-r[i,j],theta[i,j],Q,B)
# Compute Minimizer
[imin,jmin]=np.unravel_index(E.argmin(),(N,N))
kappamin=r[imin,jmin]
alphamin=theta[imin,jmin]
kappa[n]=kappamin
fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(projection='polar'))
levs=np.geomspace(E.min(),E.max(),400)
pcm=ax.contourf(theta, r, E, levs, norm=colors.PowerNorm(gamma=0.2), cmap='brg')
ax.set_xticks([0,np.pi/2])
ax.set_yticks([kappamin])
colorbarticks=np.linspace(E.min(),E.max(),6)
plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1)
plt.show()
f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
f.write(str(kappa))
f.close()