#!/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 import re 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 of experiments / folders number=7 show_plot = False kappa=np.zeros(number) for n in range(0,number): # Read from Date print(str(n)) DataPath = './experiment/perforated-bilayer/results/'+str(n) DataPath = './experiment/perforated-bilayer/results/'+str(n) QFilePath = DataPath + '/QMatrix.txt' BFilePath = DataPath + '/BMatrix.txt' ParameterPath = DataPath + '/parameter.txt' # Read Thickness from parameter file (needed for energy scaling) with open(ParameterPath , 'r') as file: parameterFile = file.read() thickness = float(re.findall(r'(?m)h = (\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',parameterFile)[0]) Q, B = ReadEffectiveQuantities(QFilePath,BFilePath) # Q=0.5*(np.transpose(Q)+Q) # symmetrize B=np.transpose([B]) # N=500 length=12 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) * (thickness**2) else: E[i,j]=energy(-r[i,j],theta[i,j],Q,B) * (thickness**2) # 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) cbar = plt.colorbar(pcm, extend='max', ticks=colorbarticks, pad=0.1) cbar.ax.tick_params(labelsize=6) if (show_plot): plt.show() # Save Figure as .pdf width = 5.79 height = width / 1.618 # The golden ratio. fig.set_size_inches(width, height) fig.savefig('Plot_PerforatedBilayer.pdf') f = open("./experiment/perforated-bilayer/results/kappa_simulation.txt", "w") f.write(str(kappa)) f.close()