Forked from
Klaus Böhnlein / dune-microstructure
594 commits behind, 147 commits ahead of the upstream repository.
-
Klaus Böhnlein authoredKlaus Böhnlein authored
PhaseDiagram_Contour_v2.py 24.03 KiB
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import math
import os
import subprocess
import fileinput
import re
import matlab.engine
import sys
from ClassifyMin import *
from HelperFunctions import *
# from CellScript import *
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
from vtk.util import numpy_support
from pyevtk.hl import gridToVTK
from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import matplotlib as mpl
import seaborn as sns
import matplotlib.colors as mcolors
import time
from scipy.ndimage.filters import gaussian_filter
# from scipy import ndimage
# print(sys.executable)
# --------------------------------------------------------------------
# START :
# INPUT (Parameters): alpha, beta, theta, gamma, mu1, rho1
#
# -Option 1 : (Case lambda = 0 => q12 = 0)
# compute q1,q2,b1,b2 from Formula
# Option 1.1 :
# set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
# Option 1.2 :
# compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
# -Option 2 :
# compute Q_hom & B_eff by running 'Cell-Problem'
#
# -> CLASSIFY ...
#
# OUTPUT: Minimizer G, angle , type, curvature
# -----------------------------------------------------------------------
#
#
# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
# OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
# # ------------------------------------ get mu_gamma ------------------------------
# # ---Scenario 1.1: extreme regimes
# if gamma == '0':
# print('extreme regime: gamma = 0')
# mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
# print("mu_gamma:", mu_gamma)
# elif gamma == 'infinity':
# print('extreme regime: gamma = infinity')
# mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta) # = q1
# print("mu_gamma:", mu_gamma)
# else:
# # --- Scenario 1.2: compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
# # print("Run computeMuGamma for Gamma = ", gamma)
# with open(InputFilePath, 'r') as file:
# filedata = file.read()
# filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
# # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
# filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
# filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
# filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
# filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
# f = open(InputFilePath,'w')
# f.write(filedata)
# f.close()
# # --- Run Cell-Problem
#
# # Check Time
# # t = time.time()
# # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
# # capture_output=True, text=True)
# # --- Run Cell-Problem_muGama -> faster
# # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
# # capture_output=True, text=True)
# # --- Run Compute_muGamma (2D Problem much much faster)
#
# subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
# capture_output=True, text=True)
# # print('elapsed time:', time.time() - t)
#
# #Extract mu_gamma from Output-File TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
# with open(OutputFilePath, 'r') as file:
# output = file.read()
# tmp = re.search(r'(?m)^mu_gamma=.*',output).group() # Not necessary for Intention of Program t output Minimizer etc.....
# s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
# mu_gamma = float(s[0])
# # print("mu_gamma:", mu_gammaValue)
# # --------------------------------------------------------------------------------------
# return mu_gamma
#
# ----------- SETUP PATHS
# InputFile = "/inputs/cellsolver.parset"
# OutputFile = "/outputs/output.txt"
InputFile = "/inputs/computeMuGamma.parset"
OutputFile = "/outputs/outputMuGamma.txt"
# --------- Run from src folder:
path_parent = os.path.dirname(os.getcwd())
os.chdir(path_parent)
path = os.getcwd()
print(path)
InputFilePath = os.getcwd()+InputFile
OutputFilePath = os.getcwd()+OutputFile
print("InputFilepath: ", InputFilePath)
print("OutputFilepath: ", OutputFilePath)
print("Path: ", path)
# -------------------------- Input Parameters --------------------
# mu1 = 10.0 # TODO : here must be the same values as in the Parset for computeMuGamma
mu1 = 1.0
rho1 = 1.0
alpha = 2.0
beta = 2.0
# beta = 5.0
theta = 1.0/4.0
#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
gamma = '0'
# gamma = 'infinity'
# gamma = 0.5
# gamma = 0.25
# gamma = 1.0
# gamma = 5.0
#added
# lambda1 = 10.0
lambda1 = 0.0
#Test:
# rho1 = -1.0
print('---- Input parameters: -----')
print('mu1: ', mu1)
print('rho1: ', rho1)
print('alpha: ', alpha)
print('beta: ', beta)
print('theta: ', theta)
print('gamma:', gamma)
print('lambda1: ', lambda1)
print('----------------------------')
# ----------------------------------------------------------------
#
# gamma_min = 0.5
# gamma_max = 1.0
#
# # gamma_min = 1
# # gamma_max = 1
# Gamma_Values = np.linspace(gamma_min, gamma_max, num=3)
# # #
# # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13) # TODO variable Input Parameters...alpha,beta...
# print('(Input) Gamma_Values:', Gamma_Values)
print('type of gamma:', type(gamma))
# # #
Gamma_Values = ['0', 'infinity']
# Gamma_Values = ['infinity']
# Gamma_Values = ['0']
print('(Input) Gamma_Values:', Gamma_Values)
for gamma in Gamma_Values:
print('Run for gamma = ', gamma)
print('type of gamma:', type(gamma))
# muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
# # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
# print('Test MuGamma:', muGamma)
# ------- Options --------
# print_Cases = True
# print_Output = True
#TODO
# generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
generalCase = False
# make_3D_plot = True
# make_3D_PhaseDiagram = True
make_2D_plot = False
make_2D_PhaseDiagram = False
make_3D_plot = False
make_3D_PhaseDiagram = False
make_2D_plot = True
make_2D_PhaseDiagram = True
#
# --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
# q1 = harmonicMean(mu1, beta, theta)
# q2 = arithmeticMean(mu1, beta, theta)
# --- Set q12
# q12 = 0.0 # (analytical example) # TEST / TODO read from Cell-Output
# b1 = prestrain_b1(rho1, beta, alpha, theta)
# b2 = prestrain_b2(rho1, beta, alpha, theta)
#
# print('---- Input parameters: -----')
# print('mu1: ', mu1)
# print('rho1: ', rho1)
# print('alpha: ', alpha)
# print('beta: ', beta)
# print('theta: ', theta)
# print("q1: ", q1)
# print("q2: ", q2)
# print("mu_gamma: ", mu_gamma)
# print("q12: ", q12)
# print("b1: ", b1)
# print("b2: ", b2)
# print('----------------------------')
# print("machine epsilon", sys.float_info.epsilon)
# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
# print("Test", Test)
# ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
# SamplePoints_3D = 10 # Number of sample points in each direction
# SamplePoints_2D = 10 # Number of sample points in each direction
SamplePoints_3D = 300 # Number of sample points in each direction
# SamplePoints_3D = 150 # Number of sample points in each direction
# SamplePoints_3D = 100 # Number of sample points in each direction
# SamplePoints_3D = 200 # Number of sample points in each direction
# SamplePoints_3D = 400 # Number of sample points in each direction
# SamplePoints_2D = 7500 # Number of sample points in each direction
# SamplePoints_2D = 4000 # 4000 # Number of sample points in each direction
SamplePoints_2D = 400 # 4000 # Number of sample points in each direction
# SamplePoints_2D = 500 # 4000 # Number of sample points in each direction
# SamplePoints_2D = 100 # 4000 # Number of sample points in each direction
# SamplePoints_2D = 2000 # 4000 # Number of sample points in each direction
# SamplePoints_2D = 1000 # 4000 # Number of sample points in each direction
# SamplePoints_2D = 1500 # 4000 # Number of sa
if make_3D_PhaseDiagram:
alphas_ = np.linspace(-20, 20, SamplePoints_3D)
# alphas_ = np.linspace(-10, 10, SamplePoints_3D)
# betas_ = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
# betas_ = np.linspace(0.01,20.01,SamplePoints_3D) # FULL Range
# betas_ = np.linspace(0.01,0.99,SamplePoints_3D) # weird part
betas_ = np.linspace(1.01,40.01,SamplePoints_3D) #TEST !!!!! For Beta <1 weird tings happen...
thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
# TEST
# alphas_ = np.linspace(-2, 2, SamplePoints_3D)
# betas_ = np.linspace(1.01,10.01,SamplePoints_3D)
# print('betas:', betas_)
# TEST :
# alphas_ = np.linspace(-40, 40, SamplePoints_3D)
# betas_ = np.linspace(0.01,80.01,SamplePoints_3D) # Full Range
# print('type of alphas', type(alphas_))
# print('Test:', type(np.array([mu_gamma])) )
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
classifyMin_anaVec = np.vectorize(classifyMin_ana)
# Get MuGamma values ...
GetMuGammaVec = np.vectorize(GetMuGamma)
muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
# Classify Minimizers....
G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
# G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas, mu1, rho1, True, True)
# print('size of G:', G.shape)
# print('G:', G)
# Option to print angles
# print('angles:', angles)
# Out = classifyMin_anaVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
GammaString = str(gamma)
VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
print('Written to VTK-File:', VTKOutputName )
if make_2D_PhaseDiagram:
# alphas_ = np.linspace(-20, 20, SamplePoints_2D)
# alphas_ = np.linspace(0, 1, SamplePoints_2D)
thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
alphas_ = np.linspace(-5, 5, SamplePoints_2D)
# alphas_ = np.linspace(-5, 15, SamplePoints_2D)
# thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
# good range:
# alphas_ = np.linspace(9, 10, SamplePoints_2D)
# thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
# range used:
# alphas_ = np.linspace(8, 10, SamplePoints_2D)
# thetas_ = np.linspace(0.05,0.16,SamplePoints_2D)
# alphas_ = np.linspace(8, 12, SamplePoints_2D)
# thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
# betas_ = np.linspace(0.01,40.01,1)
#fix to one value:
betas_ = 2.0;
# betas_ = 10.0;
# betas_ = 5.0;
# betas_ = 0.5;
#intermediate Values
alphas_ = np.linspace(-2, 1, SamplePoints_2D)
# alphas_ = np.linspace(-1.5, 1, SamplePoints_2D)
# thetas_ = np.linspace(0.4,0.6,SamplePoints_2D)
# betas_ = 10.0;
# TEST
# alphas_ = np.linspace(-8, 8, SamplePoints_2D)
# thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
# betas_ = 1.0; #TEST Problem: disvison by zero if alpha = 9, theta = 0.1 !
# betas_ = 0.9;
# betas_ = 0.5; #TEST!!!
# alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
betas = betas_
alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
if generalCase:
classifyMin_matVec = np.vectorize(classifyMin_mat)
GetCellOutputVec = np.vectorize(GetCellOutput, otypes=[np.ndarray, np.ndarray])
Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,lambda1, InputFilePath ,OutputFilePath )
# print('type of Q:', type(Q))
# print('Q:', Q)
G, angles, Types, curvature = classifyMin_matVec(Q,B)
else:
classifyMin_anaVec = np.vectorize(classifyMin_ana)
GetMuGammaVec = np.vectorize(GetMuGamma)
# muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
# G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
if gamma == '0':
G, angles_0, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
if gamma == 'infinity':
G, angles_inf, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!!
# print('size of G:', G.shape)
# print('G:', G)
# print('Types:', Types)
# Out = classifyMin_anaVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
# VTKOutputName = + path + "./PhaseDiagram2DNEW"
# print('angles:',angles)
# GammaString = str(gamma)
# VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
# gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
# print('Written to VTK-File:', VTKOutputName )
# --- Make 3D Scatter plot
if(make_3D_plot or make_2D_plot):
# Styling
plt.style.use("seaborn-darkgrid")
plt.style.use("seaborn-whitegrid")
plt.style.use("seaborn")
# plt.style.use("seaborn-paper")
# plt.style.use('ggplot')
# plt.rcParams["font.family"] = "Avenir"
# plt.rcParams["font.size"] = 16
# plt.style.use("seaborn-darkgrid")
mpl.rcParams['text.usetex'] = True
mpl.rcParams["font.family"] = "serif"
mpl.rcParams["font.size"] = "10"
# mpl.rcParams['xtick.labelsize'] = 16mpl.rcParams['xtick.major.size'] = 2.5
# mpl.rcParams['xtick.bottom'] = True
# mpl.rcParams['ticks'] = True
mpl.rcParams['xtick.bottom'] = True
mpl.rcParams['xtick.major.size'] = 3
mpl.rcParams['xtick.minor.size'] = 1.5
mpl.rcParams['xtick.major.width'] = 0.75
mpl.rcParams['ytick.left'] = True
mpl.rcParams['ytick.major.size'] = 3
mpl.rcParams['ytick.minor.size'] = 1.5
mpl.rcParams['ytick.major.width'] = 0.75
mpl.rcParams.update({'font.size': 10})
### ADJUST GRID:
mpl.rcParams['axes.labelpad'] = 5
mpl.rcParams['grid.linewidth'] = 0.25
mpl.rcParams['grid.alpha'] = 0.9 # 0.75
mpl.rcParams['grid.linestyle'] = '-'
mpl.rcParams['grid.color'] = 'gray'#'black'
# mpl.rcParams['axes.axisbelow'] = True
# mpl.rcParams.update({"axes.grid" : True, "grid.color": "gray"})
# mpl.rcParams["axes.grid"] = False
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# colors = cm.plasma(Types)
colors = cm.coolwarm(angles_inf)
width = 6.28*0.5
# height = width / 1.618
# height = width / 2.5
height = width
# fig, ax = plt.subplots()
# fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(width,height), sharey=True)
fig,ax = plt.subplots(nrows=1,ncols=1,figsize=(width,height))
# ax = plt.axes((0.15,0.21 ,0.8,0.75))
# if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
# if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
#
# if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
# if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
# pnt=ax.scatter(alphas,thetas,c=angles,cmap='coolwarm')
# # ax.colorbar()
# CS = ax.contourf(alphas, thetas, angles,6, cmap=plt.cm.coolwarm, linestyle=dashed)
# # CS = ax.contour(alphas, thetas, angles,6, colors='k')
# ax.clabel(CS, inline=True, fontsize=7.5)
# # ax.set_title('Simplest default with labels')
# matplotlib.rcParams['contour.linestyles'] = 'dashed'
cmap=mpl.cm.coolwarm
# cmap = sns.color_palette("flare", as_cmap=True)
# cmap = sns.color_palette('species')
### GET COLORS :
deep_colors = sns.color_palette("pastel")
print('deep_colors.as_hex():',deep_colors.as_hex())
diverging_colors = sns.color_palette("RdBu", 10)
print('diverging_colors.as_hex():',diverging_colors.as_hex())
pal = sns.color_palette("Blues")
pal = sns.color_palette()
print(pal.as_hex())
# flatui = ["#9b59b6", "#3498db", "#95a5a6", "#e74c3c", "#34495e", "#2ecc71"]
flatui = ["coral","white", "cornflowerblue"]
flatui = ["cornflowerblue", "coral"]
flatui = ['#4c72b0','white', '#c44e52']
flatui = ['#4c72b0','white', '#8de5a1']
flatui = ['#a1c9f4', '#ffb482','#ff9f9b'] #Test colors
flatui = ['#4c72b0','white', '#ffb482']
flatui = ['#4c72b0','white', '#ff9f9b']
flatui = ['#4c72b0','white', '#ab162a']
# flatui = ['#4c72b0','white', '#eb9172']
# flatui = ['#4c72b0','white', '#64b5cd']
cmap = mpl.colors.ListedColormap(sns.color_palette(flatui).as_hex())
cmap = mpl.colors.ListedColormap(sns.color_palette(flatui).as_hex())
cmap = mpl.colors.ListedColormap(sns.color_palette("RdBu_r", 10).as_hex())
cmap = mpl.colors.ListedColormap(sns.color_palette("coolwarm", 10).as_hex()) #Discrete CMAP
cmap = sns.color_palette("coolwarm", as_cmap=True)
# cmap = sns.color_palette("vlag", as_cmap=True)
# cmap = sns.color_palette("icefire", as_cmap=True)
# cmap = sns.color_palette("Spectral_r", as_cmap=True)
# cmap = sns.color_palette("flare_r", as_cmap=True)
# cmap = sns.diverging_palette(220, 20, as_cmap=True)
# cmap = sns.diverging_palette(250, 30, l=65, center="dark", as_cmap=True)
# cmap = mpl.colors.ListedColormap(sns.color_palette().as_hex())
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", sns.color_palette(flatui).as_hex())
# choose PLOT
gamma = '0'
# gamma = 'infinity'
### REVERSE COLORMAP :
# cmap = mpl.colors.ListedColormap(cmap.colors[::-1])
# cmap = cmap.reverse()
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["blue","violet","red"])
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["cornflowerblue","coral"])
# plt.imshow(angles_0, extent=[-2, 1, 0, 1],origin='lower',
# cmap=cmap, alpha=1.0)
if gamma == '0':
divnorm=mcolors.TwoSlopeNorm(vmin=angles_0.min(), vcenter=(angles_0.max()+angles_0.min())/2, vmax=angles_0.max())
ax.imshow(angles_0.T, extent=[-2, 1, 0, 1], origin='lower', norm = divnorm,
cmap=cmap, alpha=0.9, aspect=2.5)
levels = np.arange(0.0, 1.58, 1)
CS = ax.contour(alphas, thetas, angles_0, levels=1, colors='black',linewidths=(0.5), extent=(-2, 1, 0, 1), zorder=5)
# manual_locations = [
# (-0.4, 0.2),(-0.6, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
# ax.clabel(CS2, inline=True, fontsize=6, colors='black', manual=manual_locations)
# ax.clabel(CS2, inline=True, fontsize=6, colors='black')
# ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
# ax.clabel(CS, fontsize=5, colors='black')
# cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
# cbar = fig.colorbar(CS_0, ticks=[0, np.pi/2 ])
# cbar.ax.set_yticklabels(['$0$', r'$\pi/2$'])
# cbar.ax.set_title(r'angle $\alpha$')
if gamma == 'infinity':
divnorm=mcolors.TwoSlopeNorm(vmin=angles_inf.min(), vcenter=(angles_inf.max()+angles_inf.min())/2, vmax=angles_inf.max())
Im = ax.imshow(angles_inf.T, extent=[-2, 1, 0, 1], origin='lower', norm = divnorm,
cmap=cmap, alpha=0.9, aspect=2.5)
levels = np.arange(0.25, 1.6, 0.25)
CS_1 = ax.contour(alphas, thetas, angles_inf,levels, colors='black',linewidths=(0.5), extent=(-2, 1, 0, 1), zorder=5)
manual_locations = [
(-0.4, 0.1), (-0.6,0.30), (-0.7, 0.45), (-0.8, 0.6), (-0.9, 0.75), (-1,0.90)]
ax.clabel(CS_1, inline=True, fontsize=10, colors='black', manual=manual_locations)
axins1 = inset_axes(ax,
width="5%", # width = 5% of parent_bbox width
height="100%", # height : 50%
loc='lower left',
bbox_to_anchor=(1.05, 0., 1, 1),
bbox_transform=ax.transAxes,
borderpad=0,
)
cbar = fig.colorbar(Im, cax=axins1, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
cbar.ax.set_yticklabels(['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$'])
cbar.ax.set_title(r'angle $\alpha$', fontsize=10)
if gamma == '0':
ax.set_xlabel(r'$\theta_\rho$',fontsize=10)
# ax.yaxis.set_major_locator(MultipleLocator(0.1))
# ax.xaxis.set_major_locator(MultipleLocator(1))
ax.yaxis.set_major_locator(MultipleLocator(0.1))
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.set_ylabel(r'$\theta$ ',fontsize=10, rotation=0, labelpad=10)
ax.tick_params(axis='x', labelsize=10 )
ax.tick_params(axis='y', labelsize=10)
plt.subplots_adjust(left = 0.2)
if gamma == 'infinity':
ax.set_xlabel(r'$\theta_\rho$',fontsize=10)
# ax.xaxis.set_minor_locator(MultipleLocator(0.5))
# ax.yaxis.set_major_locator(MultipleLocator(0.1))
# ax.xaxis.set_major_locator(MultipleLocator(1))
ax.yaxis.set_major_locator(MultipleLocator(0.1))
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.tick_params(axis='x', labelsize=10 )
# ax.tick_params(axis='y', labelsize=10 )
ax.yaxis.set_ticklabels([])
plt.subplots_adjust(right= 0.75)
plt.subplots_adjust(left = 0.05)
# ax.set_ylabel('beta')
# ax.set_ylabel(r'$\theta$ ',fontsize=10, rotation=0)
# if make_3D_plot: ax.set_zlabel('theta')
# plt.subplots_adjust(bottom=0.2)
# plt.subplots_adjust(wspace=0.22, hspace=0.1)
# plt.subplots_adjust(hspace=0.15, wspace=0.1)
# plt.subplots_adjust(right= 0.75)
plt.subplots_adjust(bottom=0.2)
# fig.subplots_adjust(right=0.75)
# ax.grid( linestyle = '-', linewidth = 0.25, alpha=0.5, zorder=1)
# ax.grid( linestyle = '--', linewidth = 0.25, zorder=1)
# ax.set_axisbelow(True)
# ax.yaxis.grid(color='gray', linestyle='dashed')
# # Hide grid lines
# ax.grid(False)
fig.set_size_inches(width, height)
# outputName = '2D-PhaseDiagram-Angle.pdf'
outputName = '2D-PhaseDiagram-Angle_' + gamma + '.pdf'
fig.savefig(outputName, dpi=300, format='pdf')
# fig.savefig('Plot-Contour.pdf')
plt.show()
# plt.savefig('common_labels.png', dpi=300)
# print('T:', T)
# print('Type 1 occured here:', np.where(T == 1))
# print('Type 2 occured here:', np.where(T == 2))
# print(alphas_)
# print(betas_)
# ALTERNATIVE
# colors = ("red", "green", "blue")
# groups = ("Type 1", "Type2", "Type3")
#
# # Create plot
# fig = plt.figure()
# ax = fig.add_subplot(1, 1, 1)
#
# for data, color, group in zip(Types, colors, groups):
# # x, y = data
# ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
#
# plt.title('Matplot scatter plot')
# plt.legend(loc=2)
# plt.show()