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
from plotly.subplots import make_subplots
import plotly.io as pio
# 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()


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




# 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 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 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 = 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)

# num_Points = 400
num_Points = 100
# 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(-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

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
height = width
fig = plt.figure()

# ax = plt.axes(projection ='3d', adjustable='box')
ax = plt.axes((0.17,0.21 ,0.75,0.75))
# 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.5,),zorder=5)


## PLOTLY:

fig = make_subplots(specs=[[{"secondary_y": True}]])


data1 = go.Contour(x = x_in, y = y_in, z = Z_in,
		# colorscale='Electric',
		colorscale='agsunset',
		contours=dict(
            coloring ='heatmap',
            showlabels = True, # show labels on contours
            labelfont = dict( # label font properties
                size = 14,
                color = 'white',
            )
		), line_width=2 ,line_smoothing=0.85,
		# colorbar=dict(
        #     title='Color bar title', # title here
        #     titleside='right',
        #     titlefont=dict(
        #         size=14,
        #         family='Arial, sans-serif')
		# )


		)
epsilon = 0.4
epsilon = 0.1
# data2 = go.Scatter(x=[0+epsilon,x.max(),x.max(), 0+epsilon], y=[0-epsilon,0-epsilon,y.min(),y.min()], fill="toself")
data2 = go.Scatter(x=[0+epsilon,x.max(),x.max(), 0+epsilon], y=[0-epsilon,0-epsilon,y.min(),y.min()])

# data_tmp = [data1,data2]
data_tmp = [data1]

fig = go.Figure(data = data_tmp)

# ADD RECTANGLES:
fig.add_shape(type="rect",
    x0=0+epsilon, y0=0-epsilon, x1=x.max(), y1=y.min(),
    line=dict(
        color="White",
        width=1,
    ),
    fillcolor="White",
    opacity=0.75
)
fig.add_shape(type="rect",
    x0=x.min(), y0=0+epsilon, x1=0-epsilon, y1=y.max(),
    line=dict(
        color="White",
        width=1
    ),
    fillcolor="White",
    opacity=0.75
)


# ADD COORDINATE LINES

fig.add_shape(type="line",
    x0=x.min(), y0=0, x1=x.max(), y1=0,
    line=dict(
        color="black",
        width=2,
        dash="dash",
    )
)

fig.add_shape(type="line",
    x0=0, y0=y.min(), x1=0, y1=y.max(),
    line=dict(
        color="black",
        width=2,
        dash="dash",
    )
)

# fig = go.Figure(data = go.Contour(Z_in, a1_in, a2_in))
# fig = go.Figure(data = go.Contour(x = x_in, y = y_in, z = Z_in,
# 		# colorscale='Electric',
# 		colorscale='agsunset',
# 		contours=dict(
#             coloring ='heatmap',
#             showlabels = True, # show labels on contours
#             labelfont = dict( # label font properties
#                 size = 12,
#                 color = 'white',
#             )
# 		), line_width=2 ,line_smoothing=0.85,
# 		# colorbar=dict(
#         #     title='Color bar title', # title here
#         #     titleside='right',
#         #     titlefont=dict(
#         #         size=14,
#         #         family='Arial, sans-serif')
# 		# )
# 		)
#         )

epsilon = 0.4
epsilon = 0.1
# 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


# trace2 = go.Figure(go.Scatter(x=[0+epsilon,x.max(),x.max(), 0+epsilon], y=[0-epsilon,0-epsilon,y.min(),y.min()], fill="toself"))
# fig.add_trace(trace1)

# fig.add_trace(trace1)
# fig.add_trace(trace2.data[0],secondary_y=True)

fig.update_layout(
    autosize=False,
    # width=500,
    # height=500,
    width=300,
    height=300,
    # margin=dict(
    #     l=50,
    #     r=50,
    #     b=100,
    #     t=100,
    #     pad=4
    # ),
    # paper_bgcolor="LightSteelBlue",
)

width = 452.9579/2
size= set_size(width, fraction=0.5)
print('set_size(width, fraction=0.5)', set_size(width, fraction=0.5))
print('size[0]',size[0])



# fig.show()
# show(fig)
# fig.write_image("Plotly-fig1.png",  width=width, height=height, scale=1)
fig.write_image("Plotly-fig1.png")
# pio.write_image(fig, "test.svg", width=6.28 *0.33*300, height=6.28 *0.33*300, scale=1)
# pio.write_image(fig, "test.png", width=6.28*0.33*300, height=6.28*0.33*300)
pio.write_image(fig, "test.pdf", width=6.28*0.33*300, height=6.28*0.33*300)
# pio.write_image(fig, "test.pdf", width=width, height=width)
# 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)


# 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 = cm.brg(Z)
divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=0., vmax=Z.max())

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

plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
                  cmap=cmap, alpha=0.6)





# COLORBAR :
cbar = plt.colorbar()
cbar.ax.tick_params(labelsize=8)




##----- ADD RECTANGLE TO COVER QUADRANT :
epsilon = 0.4
epsilon = 0.1
# 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.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.set_xlabel(r"$a_1$", fontsize=10 ,labelpad=0)
ax.set_ylabel(r"$a_2$", fontsize=10 ,labelpad=0)
# ax.set_ylabel(r"energy")


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