Forked from
Klaus Böhnlein / dune-microstructure
474 commits behind, 209 commits ahead of the upstream repository.
-
Klaus Böhnlein authoredKlaus Böhnlein authored
Energy_ContourG+_v2.py 26.12 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
import time
import matplotlib.ticker as ticker
import matplotlib as mpl
from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
import pandas as pd
import seaborn as sns
import matplotlib.colors as mcolors
from chart_studio import plotly
import plotly.graph_objs as go
import plotly.express as px
import plotly.colors
# from matplotlib import rc
# rc('text', usetex=True) # Use LaTeX font
#
# import seaborn as sns
# sns.set(color_codes=True)
def show(fig):
import io
import plotly.io as pio
from PIL import Image
buf = io.BytesIO()
pio.write_image(fig, buf)
img = Image.open(buf)
img.show()
# set the colormap and centre the colorbar
class MidpointNormalize(mcolors.Normalize):
"""
Normalise the colorbar so that diverging bars work there way either side from a prescribed midpoint value)
e.g. im=ax1.imshow(array, norm=MidpointNormalize(midpoint=0.,vmin=-100, vmax=100))
"""
def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
self.midpoint = midpoint
mcolors.Normalize.__init__(self, vmin, vmax, clip)
def __call__(self, value, clip=None):
# I'm ignoring masked values and all kinds of edge cases to make a
# simple example...
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
return np.ma.masked_array(np.interp(value, x, y), np.isnan(value))
def set_size(width, fraction=1):
"""Set figure dimensions to avoid scaling in LaTeX.
Parameters
----------
width: float
Document textwidth or columnwidth in pts
fraction: float, optional
Fraction of the width which you wish the figure to occupy
Returns
-------
fig_dim: tuple
Dimensions of figure in inches
"""
# Width of figure (in pts)
fig_width_pt = width * fraction
# Convert from pt to inches
inches_per_pt = 1 / 72.27
# Golden ratio to set aesthetic figure height
# https://disq.us/p/2940ij3
golden_ratio = (5**.5 - 1) / 2
# Figure width in inches
fig_width_in = fig_width_pt * inches_per_pt
# Figure height in inches
fig_height_in = fig_width_in * golden_ratio
fig_dim = (fig_width_in, fig_height_in)
return fig_dim
def format_func(value, tick_number):
# find number of multiples of pi/2
# N = int(np.round(2 * value / np.pi))
# if N == 0:
# return "0"
# elif N == 1:
# return r"$\pi/2$"
# elif N == -1:
# return r"$-\pi/2$"
# elif N == 2:
# return r"$\pi$"
# elif N % 2 > 0:
# return r"${0}\pi/2$".format(N)
# else:
# return r"${0}\pi$".format(N // 2)
##find number of multiples of pi/2
N = int(np.round(4 * value / np.pi))
if N == 0:
return "0"
elif N == 1:
return r"$\pi/4$"
elif N == -1:
return r"$-\pi/4$"
elif N == 2:
return r"$\pi/2$"
elif N == -2:
return r"$-\pi/2$"
elif N % 2 > 0:
return r"${0}\pi/2$".format(N)
else:
return r"${0}\pi$".format(N // 2)
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return array[idx]
def find_nearestIdx(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
def check_case(q1,q2,q12,q3,b1,b2):
# a = np.array([a1,a2])
b = np.array([b1,b2])
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
print('det(H)=', np.linalg.det(H))
# check if g* is in G^*_R^2
tmp = A.dot(b)
## compute inverse of H :
inv_H = np.linalg.inv(H)
g_star = inv_H.dot(tmp)
print('g_star=', g_star)
return g_star
def energy(a1,a2,q1,q2,q12,q3,b1,b2):
a = np.array([a1,a2])
b = np.array([b1,b2])
H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
tmp = H.dot(a)
# print('H',H)
# print('A',A)
# print('b',b)
# print('a',a)
# print('tmp',tmp)
tmp = (1/2)*a.dot(tmp)
# print('tmp',tmp)
tmp2 = A.dot(b)
# print('tmp2',tmp2)
tmp2 = 2*a.dot(tmp2)
# print('tmp2',tmp2)
energy = tmp - tmp2
# print('energy',energy)
# energy_axial1.append(energy_1)
return energy
# def energy(a1,a2,q1,q2,q12,q3,b1,b2):
#
#
# b = np.array([b1,b2])
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
#
#
# tmp = H.dot(a)
#
# print('H',H)
# print('A',A)
# print('b',b)
# print('a',a)
# print('tmp',tmp)
#
# tmp = (1/2)*a.dot(tmp)
# print('tmp',tmp)
#
# tmp2 = A.dot(b)
# print('tmp2',tmp2)
# tmp2 = 2*a.dot(tmp2)
#
# print('tmp2',tmp2)
# energy = tmp - tmp2
# print('energy',energy)
#
#
# # energy_axial1.append(energy_1)
#
# return energy
#
################################################################################################################
################################################################################################################
################################################################################################################
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)
print('---- Input parameters: -----')
# q1=1;
# q2=2;
# q12=1/2;
# q3=((4*q1*q2)**0.5-q12)/2;
# # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
#
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
# abar = np.array([q12+2*q3, 2*q2])
# abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar
#
# print('abar:',abar)
#
# b = np.linalg.lstsq(A, abar)[0]
# print('b',b)
#
#
# # print('abar:',np.shape(abar))
# # print('np.transpose(abar):',np.shape(np.transpose(abar)))
# sstar = (1/(q1+q2))*abar.dot(A.dot(b))
# # sstar = (1/(q1+q2))*abar.dot(tmp)
# print('sstar', sstar)
# abarperp= np.array([abar[1],-abar[0]])
# print('abarperp:',abarperp)
# -------------------------- Input Parameters --------------------
mu1 = 1.0
rho1 = 1.0
alpha = 5.0
theta = 1.0/2
# theta= 0.1
beta = 5.0
# mu1 = 1.0
# rho1 = 1.0
# alpha = -0.75
# theta = 1.0/2
# # theta= 0.1
# beta = 5.0
# mu1 = 1.0
# rho1 = 1.0
# alpha = 2.0
# theta = 1.0/2
# # theta= 0.1
# beta = 5.0
#Figure3:
# mu1 = 1.0
# rho1 = 1.0
# alpha = 2.0
# theta = 1.0/8
# # theta= 0.1
# beta = 2.0
# alpha= -5
#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
gamma = '0'
# gamma = 'infinity'
lambda1 = 0.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('----------------------------')
# ----------------------------------------------------------------
print('----------------------------')
# ----------------------------------------------------------------
q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
q12 = 0.0
q3 = GetMuGamma(beta, theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
b1 = prestrain_b1(rho1,beta, alpha, theta )
b2 = prestrain_b2(rho1,beta, alpha, theta )
## ---- 1-ParameterFamilyCase:
# q1=1;
# q2=2;
# q12=1/2;
# q3=((4*q1*q2)**0.5-q12)/2;
# # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
# abar = np.array([q12+2*q3, 2*q2])
# abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar
# print('abar:',abar)
# b = np.linalg.lstsq(A, abar)[0]
# print('b',b)
# b1=b[0]
# b2=b[1]
## ---------------
print('q1 = ', q1)
print('q2 = ', q2)
print('q3 = ', q3)
print('q12 = ', q12)
print('b1 = ', b1)
print('b2 = ', b2)
## --- CHECK CASE ---
g_star = check_case(q1,q2,q12,q3,b1,b2)
## -------------------------------------
num_Points = 400
num_Points = 200
# num_Points = 20
# Creating dataset
x = np.linspace(-5,5,num_Points)
y = np.linspace(-5,5,num_Points)
x = np.linspace(-20,20,num_Points)
y = np.linspace(-20,20,num_Points)
x = np.linspace(-2,2,num_Points)
y = np.linspace(-2,2,num_Points)
# x = np.linspace(-10,10,num_Points)
# y = np.linspace(-10,10,num_Points)
# x = np.linspace(-60,60,num_Points)
# y = np.linspace(-60,60,num_Points)
#
#
# x = np.linspace(-40,40,num_Points)
# y = np.linspace(-40,40,num_Points)
a1, a2 = np.meshgrid(x,y)
# geyser = sns.load_dataset("geyser")
# print('type of geyser:', type(geyser))
# print('geyser:',geyser)
# ContourRange=20
ContourRange=2
x_in = np.linspace(-ContourRange,ContourRange,num_Points)
y_in = np.linspace(-ContourRange,ContourRange,num_Points)
a1_in, a2_in = np.meshgrid(x_in,y_in)
print('a1:', a1)
print('a2:',a2 )
print('a1.shape', a1.shape)
#-- FILTER OUT VALUES for G+ :
# tmp1 = a1[np.where(a1*a2 >= 0)]
# tmp2 = a2[np.where(a1*a2 >= 0)]
#
# np.take(a, np.where(a>100)[0], axis=0)
# tmp1 = np.take(a1, np.where(a1*a2 >= 0)[0], axis=0)
# tmp2 = np.take(a1, np.where(a1*a2 >= 0)[0], axis=0)
# tmp2 = a2[np.where(a1*a2 >= 0)]
tmp1 = a1[a1*a2 >= 0]
tmp2 = a2[a1*a2 >= 0]
tmp1 = tmp1.reshape(-1,5)
tmp2 = tmp2.reshape(-1,5)
# tmp1_pos = a1[np.where(a1*a2 >= 0) ]
# tmp2_pos = a2[np.where(a1*a2 >= 0) ]
# tmp1_pos = tmp1_pos[np.where(tmp1_pos >= 0)]
# tmp2_pos = tmp2_pos[np.where(tmp2_pos >= 0)]
#
# tmp1_neg = a1[a1*a2 >= 0 ]
# tmp2_neg = a2[a1*a2 >= 0 ]
# tmp1_neg = tmp1_neg[tmp1_neg < 0]
# tmp2_neg = tmp2_neg[tmp2_neg < 0]
# a1 = tmp1
# a2 = tmp2
#
# a1 = a1.reshape(-1,5)
# a2 = a2.reshape(-1,5)
#
# tmp1_pos = tmp1_pos.reshape(-1,5)
# tmp2_pos = tmp2_pos.reshape(-1,5)
# tmp1_neg = tmp1_neg.reshape(-1,5)
# tmp2_neg = tmp2_neg.reshape(-1,5)
# print('a1:', a1)
# print('a2:',a2 )
# print('a1.shape', a1.shape)
energyVec = np.vectorize(energy)
# Z = energyVec(np.array([a1,a2]),q1,q2,q12,q3,b1,b2)
Z = energyVec(a1,a2,q1,q2,q12,q3,b1,b2)
Z_in = energyVec(a1_in,a2_in,q1,q2,q12,q3,b1,b2)
print('Z:', Z)
print('any', np.any(Z<0))
#
# negZ_a1 = a1[np.where(Z<0)]
# negZ_a2 = a2[np.where(Z<0)]
# negativeValues = Z[np.where(Z<0)]
# print('negativeValues:',negativeValues)
#
# print('negZ_a1',negZ_a1)
# print('negZ_a2',negZ_a2)
#
#
# negZ_a1 = negZ_a1.reshape(-1,5)
# negZ_a2 = negZ_a2.reshape(-1,5)
# negativeValues = negativeValues.reshape(-1,5)
#
# Z_pos = energyVec(tmp1_pos,tmp2_pos,q1,q2,q12,q3,b1,b2)
# Z_neg = energyVec(tmp1_neg,tmp2_neg,q1,q2,q12,q3,b1,b2)
# print('Test energy:' , energy(np.array([1,1]),q1,q2,q12,q3,b1,b2))
# print('Z_pos.shape', Z_pos.shape)
## -- PLOT :
mpl.rcParams['text.usetex'] = True
mpl.rcParams["font.family"] = "serif"
mpl.rcParams["font.size"] = "9"
label_size = 8
mpl.rcParams['xtick.labelsize'] = label_size
mpl.rcParams['ytick.labelsize'] = label_size
# plt.style.use('seaborn')
# plt.style.use('seaborn-whitegrid')
# sns.set()
# plt.style.use('seaborn-whitegrid')
label_size = 9
mpl.rcParams['xtick.labelsize'] = label_size
mpl.rcParams['ytick.labelsize'] = label_size
width = 6.28 *0.5
# width = 6.28 *0.33
# width = 6.28
height = width #/ 1.618
# width = 452.9579/2
# size= set_size(width, fraction=0.5)
# print('set_size(width, fraction=0.5)', set_size(width, fraction=1))
# print('size[0]',size[0])
f_size = 8
# fig= plt.figure()
fig, ax = plt.subplots()
# fig.set_size_inches(width, height)
# fig.set_size_inches(set_size(width, fraction=0.5))
# ax = plt.axes(projection ='3d', adjustable='box')
# ax = plt.axes((0.17,0.21 ,0.75,0.75))
# ax = plt.axes((0.17,0.23 ,0.7,0.7))
# ax = plt.axes((0.17,0.23 ,1.0,1.0))
# ax=plt.axes()
# ax = plt.axes((0.15,0.18,0.8,0.8))
# ax.tick_params(axis='x',which='major', direction='out',pad=5)
# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
# ax.xaxis.set_major_locator(MultipleLocator(0.1))
# ax.xaxis.set_minor_locator(MultipleLocator(0.05))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
# ax.grid(True,which='major',axis='both',alpha=0.3)
# colorfunction=(B*kappa)
# print('colofunction',colorfunction)
#translate Data
# Z = Z - (Z.max()-Z.min())/2
# Z = Z - 50
# Z = Z - 500
#
# Z = Z.T
# Substract constant:
# c = (b1**2)*q1+b1*b2*q12+(b2**2)*q2
# Z = Z-c
#
# print('Value of c:', c)
print('Z.min()', Z.min())
print('Z.max()', Z.max())
norm=mcolors.Normalize(Z.min(),Z.max())
# facecolors=cm.brg(norm)
print('norm:', norm)
print('type of norm', type(norm))
print('norm(0):', norm(0))
print('norm(Z):', norm(Z))
# ax.plot(theta_rho, theta_values, 'royalblue', zorder=3, )
# ax.scatter(a1,a2, s=0.5)
# ax.scatter(tmp1_pos,tmp2_pos, s=0.5)
# ax.scatter(tmp1_neg,tmp2_neg, s=0.5)
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, levels=100 )
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, levels=20 )
# sns.kdeplot(np.array([a1, a2, Z]))
# sns.kdeplot(tmp1_pos,tmp2_pos,Z_pos)
# levels = [-5.0, -4, -3, 0.0, 1.5, 2.5, 3.5]
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, corner_mask=True,levels=levels)
# CS = ax.contour(a1, a2, Z, cmap=plt.cm.gnuplot(norm(Z)), corner_mask=True)
# CS = ax.contour(a1, a2, Z, cm.brg(norm(Z)), levels=20)
# CS = ax.contour(a1, a2, Z, cmap=plt.cm.gnuplot, levels=20)
# CS = ax.contour(a1, a2, Z, colors='k', levels=14, linewidths=(0.5,))
# CS = ax.contour(a1, a2, Z, colors='k', levels=18, linewidths=(0.5,))
# ax.contour(negZ_a1, negZ_a2, negativeValues, colors='k', linewidths=(0.5,))
# CS = ax.contour(a1_in, a2_in, Z_in, colors='k', linewidths=(0.75) ,zorder=5)
step = (-1)*Z.min()
levels = np.arange(Z.min(),Z.max()/2,step)
levels = np.arange(Z.min(),Z.max()/12,step)
# levels = np.arange(0,Z.max()/2,200)
# CS = ax.contour(a1_in, a2_in, Z_in, colors='k', linewidths=(0.75), levels=levels ,zorder=5)
# CS = ax.contour(a1, a2, Z, colors='k', linewidths=(0.75), levels=levels ,zorder=5)
CS = ax.contour(a1_in, a2_in, Z_in, colors='k', linewidths=(1) , vmin= Z.min()+0.04, zorder=5)
# df = pd.DataFrame(data=Z_in, columns=a1_in, index=a2_in)
# df2 = pd.DataFrame(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
# columns=['a', 'b', 'c'])
# sns.kdeplot(data=df2, x="waiting", y="duration")
# sns.kdeplot(data=df2)
# CS = ax.contour(a1, a2, Z, colors='k', linewidths=(0.5,))
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, extend='both', levels=50)
# CS = ax.contourf(a1, a2, Z,10, colors='k', extend='both', levels=50)
# CS = ax.contourf(a1, a2, Z,10, colors='k')
#
# # CS = ax.contour(tmp1_pos,tmp2_pos, Z_pos,10, cmap=plt.cm.gnuplot, levels=10 )
# # CS = ax.contour(tmp1_pos,tmp2_pos, Z_pos,10, cmap=plt.cm.gnuplot, corner_mask=True)
#
# CS = ax.contour(a1, a2, Z,10, colors = 'k')
# ax.clabel(CS, inline=True, fontsize=4)
# ax.clabel(CS, inline=True, fontsize=4)
# 1- ParameterFAMILY CASE :
# manual_locations = [
# (1, 1), (-5, -5), (-10, -10 ),(-12.5,-12.5),(-14,-14), (-15,-15),
# (5, 5), (10, 10 ),(12.5,12.5), (15,15), (17,17)]
# GAMMA = inf
manual_locations = [
(1, 1), (-5, -5), (-10, -10 ), (-15,-15),
(5, 5), (10, 10 ),(12.5,12.5), (15,15), (17,17)]
#
#
# # GAMMA = 0
# manual_locations = [
# (1, 1), (-5, -5), (-10, -10 ), (-15,-15), (-15,7),
# (5, 5), (10, 10 ), (15,15), (17,17)]
# ax.clabel(CS, inline=True, fontsize=f_size, colors='black', manual=manual_locations)
ax.clabel(CS, inline=True, fontsize=f_size, colors='black')
# cmap = cm.brg(norm(Z))
#
# C_map = cm.inferno(norm(Z))
# ax.imshow(Z, cmap=C_map, extent=[-20, 20, -20, 20], origin='lower', alpha=0.5)
# ax.imshow(norm(Z), extent=[-20, 20, -20, 20], origin='lower',
# cmap='bwr', alpha=0.8)
# ax.imshow(norm(Z), extent=[-20, 20, -20, 20],origin='lower', vmin=Z.min(), vmax=Z.max(),
# cmap='bwr', alpha=0.6)
# ax.imshow(norm(Z), extent=[-20, 20, -20, 20],origin='lower', norm = norm,
# cmap='coolwarm', alpha=0.6)
cmap=mpl.cm.RdBu_r
cmap=mpl.cm.viridis_r
# cmap=mpl.cm.bwr
cmap=mpl.cm.coolwarm
# cmap=mpl.cm.gnuplot
cmap=mpl.cm.magma_r
# cmap=mpl.cm.inferno_r
# cmap=mpl.cm.plasma
# cmap=mpl.cm.plasma_r
# cmap=mpl.cm.cividis_r
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["blue","violet","red"])
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["blue","orange"])
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", [(0,"red"), (.1,"violet"), (.5, "blue"), (1.0, "green")])
# make a colormap that has land and ocean clearly delineated and of the
# same length (256 + 256)
#
# colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 256))
# colors_land = plt.cm.terrain(np.linspace(0.25, 1, 256))
# all_colors = np.vstack((colors_undersea, colors_land))
# # cmap = mcolors.LinearSegmentedColormap.from_list(
# # 'terrain_map', all_colors)
# cmap = px.colors.sequential.agsunset
# cmap = plotly.colors.PLOTLY_SCALES["Viridis"]
# cmap = cm.brg(Z)
divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=0., vmax=Z.max())
divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=(Z.max()+Z.min())/2, vmax=Z.max())
# divnorm=mcolors.TwoSlopeNorm(vmin=-500, vcenter=0, vmax=Z.max())
# divnorm=mcolors.TwoSlopeNorm(vmin=-10, vcenter=0. ,vmax=10)
# divnorm=mcolors.TwoSlopeNorm(vmin=-10, vcenter=0., vmax=Z.max())
# divnorm=mcolors.LogNorm(vmin=Z.min(), vmax=Z.max()) #Test LogNorm
# cmap = cm.brg(divnorm(Z))
# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower', norm = norm,
# cmap='coolwarm', alpha=0.6)
# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower',
# cmap='coolwarm', alpha=0.6)
# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower',
# cmap=cmap, alpha=0.6)
# divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=0., vmax=Z.max())
# plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower',
# cmap=cmap, alpha=0.6)
# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
# cmap=cmap, alpha=0.6)
plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
cmap=cmap, alpha=0.9)
# plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
# cmap=cmap, alpha=0.6)
# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower',
# cmap=cmap, alpha=0.6)
# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = mcolors.CenteredNorm(),
# cmap=cmap, alpha=0.6)
# COLORBAR :
# cbar = plt.colorbar()
# cbar.ax.tick_params(labelsize=f_size)
# fig.colorbar(I)
##----- ADD RECTANGLE TO COVER QUADRANT :
epsilon = 0.4
epsilon = 0.2
# ax.axvspan(0, x.max(), y.min(), 0, alpha=1, color='yellow', zorder=5)#yellow
# ax.fill_between([0, x.max()], y.min(), 0, alpha=0.3, color='yellow', zorder=5)#yellow
# ax.fill_between([x.min(), 0], 0, y.max(), alpha=0.3, color='yellow', zorder=5)#yellow
fillcolor = 'white'
# ax.fill_between([0+epsilon, x.max()], y.min(), 0-epsilon, alpha=0.7, color=fillcolor, zorder=4)#yellow
# ax.fill_between([x.min(), 0-epsilon], 0+epsilon, y.max(), alpha=0.7, color=fillcolor, zorder=4)#yellow
ax.fill_between([0+epsilon, x.max()], y.min(), 0-epsilon, alpha=1.0, color=fillcolor, zorder=4)#yellow
ax.fill_between([x.min(), 0-epsilon], 0+epsilon, y.max(), alpha=1.0, color=fillcolor, zorder=4)#yellow
# ax.plot_surface(a1,a2, Z, cmap=cm.coolwarm,
# linewidth=0, antialiased=False)
# ax.plot(theta_rho, energy_axial1, 'royalblue', zorder=3, label=r"axialMin1")
# ax.plot(theta_rho, energy_axial2, 'forestgreen', zorder=3, label=r"axialMin2")
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )
# lg = ax.legend(bbox_to_anchor=(0.0, 0.75), loc='upper left')
### PLot x and y- Axes
ax.plot(ax.get_xlim(),[0,0],'k--', linewidth=0.5)
ax.plot([0,0],ax.get_ylim(), 'k--', linewidth=0.5)
ax.scatter(g_star[0],g_star[1], s=20, zorder=5)
ax.text(g_star[0]+1,g_star[1]-1, r"$g_*$", color='royalblue', size=15, zorder = 5)
ax.set_xlabel(r"$a_1$", fontsize=f_size ,labelpad=0)
ax.set_ylabel(r"$a_2$", fontsize=f_size ,labelpad=0)
# ax.set_ylabel(r"energy")
ax.tick_params(axis='both', which='major', labelsize=f_size)
ax.tick_params(axis='both', which='minor', labelsize=f_size)
# ax.set_xticks([-np.pi/2, -np.pi/4 ,0, np.pi/4, np.pi/2 ])
# labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']
# ax.set_yticklabels(labels)
# ax.legend(loc='upper right')
fig.set_size_inches(width, height)
# fig.set_size_inches(set_size(width, fraction=0.33))
fig.savefig('Energy_ContourG+.pdf')
plt.show()
#
#
#
# # Curve parametrised by \theta_rho = alpha in parameter space
# N=100;
# theta_rho = np.linspace(1, 3, num=N)
# print('theta_rho:', theta_rho)
#
#
# theta_values = []
#
#
# for t in theta_rho:
#
# s = (1.0/10.0)*t+0.1
# theta_values.append(s)
#
#
#
#
#
# theta_rho = np.array(theta_rho)
# theta_values = np.array(theta_values)
#
# betas_ = 2.0
#
# alphas, betas, thetas = np.meshgrid(theta_rho, betas_, theta_values, indexing='ij')
#
#
# harmonicMeanVec = np.vectorize(harmonicMean)
# arithmeticMeanVec = np.vectorize(arithmeticMean)
# prestrain_b1Vec = np.vectorize(prestrain_b1)
# prestrain_b2Vec = np.vectorize(prestrain_b2)
#
# GetMuGammaVec = np.vectorize(GetMuGamma)
# muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
#
# q1_vec = harmonicMeanVec(mu1, betas, thetas)
# q2_vec = arithmeticMeanVec(mu1, betas, thetas)
#
# b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
# b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)
# special case: q12 == 0!! .. braucht eigentlich nur b1 & b2 ...
# print('type b1_values:', type(b1_values))
# print('size(q1)',q1.shape)
#
#
# energy_axial1 = []
# energy_axial2 = []
#
# # for b1 in b1_values:
# for i in range(len(theta_rho)):
# print('index i:', i)
#
# print('theta_rho[i]',theta_rho[i])
# print('theta_values[i]',theta_values[i])
#
# q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta_values[i])
# q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta_values[i])
# q12 = 0.0
# q3 = GetMuGamma(beta, theta_values[i],gamma,mu1,rho1,InputFilePath ,OutputFilePath )
# b1 = prestrain_b1(rho1,beta, theta_rho[i],theta_values[i] )
# b2 = prestrain_b2(rho1,beta, theta_rho[i],theta_values[i] )
#
#
# # q2_vec = arithmeticMean(mu1, betas, thetas)
# #
# # b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
# # b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)
# print('q1[i]',q1)
# print('q2[i]',q2)
# print('q3[i]',q3)
# print('b1[i]',b1)
# print('b2[i]',b2)
# # print('q1[i]',q1[0][i])
# # print('q2[i]',q2[i])
# # print('b1[i]',b1[i])
# # print('b2[i]',b2[i])
# #compute axial energy #1 ...
#
# a_axial1 = np.array([b1,0])
# a_axial2 = np.array([0,b2])
# b = np.array([b1,b2])
#
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
#
#
# tmp = H.dot(a_axial1)
#
# print('H',H)
# print('A',A)
# print('b',b)
# print('a_axial1',a_axial1)
# print('tmp',tmp)
#
# tmp = (1/2)*a_axial1.dot(tmp)
# print('tmp',tmp)
#
# tmp2 = A.dot(b)
# print('tmp2',tmp2)
# tmp2 = 2*a_axial1.dot(tmp2)
#
# print('tmp2',tmp2)
# energy_1 = tmp - tmp2
# print('energy_1',energy_1)
#
#
# energy_axial1.append(energy_1)
#
#
# tmp = H.dot(a_axial2)
#
# print('H',H)
# print('A',A)
# print('b',b)
# print('a_axial2',a_axial2)
# print('tmp',tmp)
#
# tmp = (1/2)*a_axial2.dot(tmp)
# print('tmp',tmp)
#
# tmp2 = A.dot(b)
# print('tmp2',tmp2)
# tmp2 = 2*a_axial2.dot(tmp2)
#
# print('tmp2',tmp2)
# energy_2 = tmp - tmp2
# print('energy_2',energy_2)
#
#
# energy_axial2.append(energy_2)
#
#
#
#
#
# print('theta_values', theta_values)
#
#
#
#
#
#
#
# kappas = []
# alphas = []
# # G.append(float(s[0]))
#
#
#
#
# for t in T :
#
# abar_current = sstar*abar+t*abarperp;
# # print('abar_current', abar_current)
# abar_current[abar_current < 1e-10] = 0
# # print('abar_current', abar_current)
#
# # G = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# G = [abar_current[0], abar_current[1] , (2*abar_current[0]*abar_current[1])**0.5 ]
#
# e = [(abar_current[0]/(abar_current[0]+abar_current[1]))**0.5, (abar_current[1]/(abar_current[0]+abar_current[1]))**0.5]
# kappa = abar_current[0]+abar_current[1]
# alpha = math.atan2(e[1], e[0])
#
# print('angle current:', alpha)
#
# kappas.append(kappa)
# alphas.append(alpha)
#
#
#
# alphas = np.array(alphas)
# kappas = np.array(kappas)
#
#
# print('kappas:',kappas)
# print('alphas:',alphas)
# print('min alpha:', min(alphas))
# print('min kappa:', min(kappas))
#
# mpl.rcParams['text.usetex'] = True
# mpl.rcParams["font.family"] = "serif"
# mpl.rcParams["font.size"] = "9"
# width = 6.28 *0.5
# height = width / 1.618
# fig = plt.figure()
# # ax = plt.axes((0.15,0.21 ,0.75,0.75))
# ax = plt.axes((0.15,0.21 ,0.8,0.75))
# ax.tick_params(axis='x',which='major', direction='out',pad=5)
# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
# # ax.xaxis.set_major_locator(MultipleLocator(0.1))
# # ax.xaxis.set_minor_locator(MultipleLocator(0.05))
# # ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# # ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
# ax.grid(True,which='major',axis='both',alpha=0.3)
#
#
#
#
# ax.plot(alphas, kappas, 'royalblue', zorder=3, )
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )