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

# dataset_numbers = [0, 1, 2, 3, 4, 5]
dataset_numbers = [0]

for dataset_number in dataset_numbers:

    kappa=np.zeros(number)
    for n in range(0,number):
        #   Read from Date
        print(str(n))
        DataPath = './experiment/wood-bilayer/results_'+ str(dataset_number) + '/' +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=400
        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_dataset_' +str(dataset_number) + '_exp' +str(n) + '.pdf')

    # f = open("./experiment/wood-bilayer/results/kappa_simulation.txt", "w")
    f = open("./experiment/wood-bilayer/results_" + str(dataset_number) +  "/kappa_simulation.txt", "w")
    f.write(str(kappa))         
    f.close()