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, )