PolarPlotLocalEnergy.py

#!/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()