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
from HelperFunctions import *
from ClassifyMin import *

import matplotlib.ticker as tickers
import matplotlib as mpl
from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
import pandas as pd

# import tikzplotlib
# # from pylab import *
# from tikzplotlib import save as tikz_save


# Needed ?
mpl.use('pdf')

# from subprocess import Popen, PIPE
#import sys

###################### makePlot.py #########################
#  Generalized Plot-Script giving the option to define
#  quantity of interest and the parameter it depends on
#  to create a plot
#
#  Input: Define y & x for "x-y plot" as Strings
#  - Run the 'Cell-Problem' for the different Parameter-Points
#  (alternatively run 'Compute_MuGamma' if quantity of interest
#   is q3=muGamma for a significant Speedup)

###########################################################



# figsize argument takes inputs in inches
# and we have the width of our document in pts.
# To set the figure size we construct a function
# to convert from pts to inches and to determine
# an aesthetic figure height using the golden ratio:
# 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 == 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 == 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



# TODO
# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
# - Also Add option to plot Minimization Output


# ----- Setup Paths -----
# InputFile  = "/inputs/cellsolver.parset"
# OutputFile = "/outputs/output.txt"

InputFile  = "/inputs/computeMuGamma.parset"
OutputFile = "/outputs/outputMuGamma.txt"

# path = os.getcwd()
# InputFilePath = os.getcwd()+InputFile
# OutputFilePath = os.getcwd()+OutputFile
# --------- 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: -----')
mu1 = 1.0  #10.0
# lambda1 = 10.0
rho1 = 1.0
# alpha = 5.0
# beta = 10.0
# alpha = 2.0
# beta = 2.0
# theta = 1.0/8.0  #1.0/4.0

lambda1 = 0.0
# gamma = 1.0/4.0

# TEST:
# alpha=3.0;




# # INTERESTING!:
# alpha = 3
beta = 2.0
# theta= 1/8




#TEST
# beta=2



gamma = 'infinity'  #Elliptic Setting
gamma = '0'       #Hyperbolic Setting
# gamma = 0.01
# # gamma= 3.0
gamma = 0.5
gamma = 0.75
# # gamma = 100.0
# gamma = 3.0

Gamma_Values = [0.5, 0.75, 1.5, 3.0]
# Gamma_Values = [ 1.5, 3.0]
Gamma_Values = [3.0]
Gamma_Values = ['infinity']
print('(Input) Gamma_Values:', Gamma_Values)
# #
for gamma in Gamma_Values:

    print('mu1: ', mu1)
    print('rho1: ', rho1)
    # print('alpha: ', alpha)
    print('beta: ', beta)
    # print('theta: ', theta)
    print('gamma:', gamma)
    print('----------------------------')



    # --- define Interval of x-va1ues:
    xmin = -2.0
    # xmax = 0.41
    xmax = 3.0


    xmin = -1.5
    xmax = 2.0

    xmin = -1.0
    xmax = -0.5


    Jumps = False


    numPoints = 2000
    numPoints = 300
    # numPoints = 30
    X_Values = np.linspace(xmin, xmax, num=numPoints)
    print(X_Values)


    Y_Values = []





    Angle_Theta01 = []
    Angle_Theta025 = []
    Angle_Theta05 = []

    Angle_Theta075 = []
    Angle_Theta09 = []






    for alpha in X_Values:
        print('Situation of Lemma1.4')
        q12 = 0.0
        # q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
        # q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
        # b1 = prestrain_b1(rho1, beta, alpha,theta)
        # b2 = prestrain_b2(rho1, beta, alpha,theta)
        # b3 = 0.0
        # q3_Theta01 = GetMuGamma(beta,0.1,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
        # q3_Theta025 = GetMuGamma(beta,0.25,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
        q3_Theta05 = GetMuGamma(beta,0.5,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
        # q3_Theta075 = GetMuGamma(beta,0.75,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
        # q3_Theta09 = GetMuGamma(beta,0.9,gamma,mu1,rho1,InputFilePath ,OutputFilePath)


        # G, angle, Type, curvature = classifyMin_ana(alpha,beta,0.1, q3_Theta01,  mu1, rho1)
        # Angle_Theta01.append(angle)
        #
        # G, angle, Type, curvature = classifyMin_ana(alpha,beta,0.25, q3_Theta025,  mu1, rho1)
        # Angle_Theta025.append(angle)

        G, angle, Type, curvature = classifyMin_ana(alpha,beta,0.5,  q3_Theta05,  mu1, rho1)
        Angle_Theta05.append(angle)

        # G, angle, Type, curvature = classifyMin_ana(alpha,beta,0.75, q3_Theta075,  mu1, rho1)
        # Angle_Theta075.append(angle)
        #
        # G, angle, Type, curvature = classifyMin_ana(alpha,beta,0.9, q3_Theta09,  mu1, rho1)
        # Angle_Theta09.append(angle)




        #
        # G, angle, Type, curvature = classifyMin_ana(-0.5,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg05 .append(angle)
        # G, angle, Type, curvature = classifyMin_ana(-0.25,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg025.append(angle)
        # G, angle, Type, curvature = classifyMin_ana(3.0,beta,theta, q3,  mu1, rho1)
        # Angle_alpha3.append(angle)
        # G, angle, Type, curvature = classifyMin_ana(-1.0,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg075.append(angle)
        # G, angle, Type, curvature = classifyMin_ana(0,beta,theta, q3,  mu1, rho1)
        # Angle_alpha0.append(angle)
        # # G, angle, Type, curvature = classifyMin_ana(-0.125,beta,theta, q3,  mu1, rho1)
        # # Angle_alphaNeg0125.append(angle)
        # G, angle, Type, curvature = classifyMin_ana(-0.7,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg0125.append(angle)
        #
        # G, angle, Type, curvature = classifyMin_ana(-0.625,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg0625.append(angle)
        # G, angle, Type, curvature = classifyMin_ana(-0.875,beta,theta, q3,  mu1, rho1)
        # Angle_alphaNeg0875.append(angle)

    #
    #
    # print("(Output) Values of angle: ", Y_Values)
    #
    #
    # idx = find_nearestIdx(Y_Values, 0)
    # print(' Idx of value  closest to 0', idx)
    # ValueClose = Y_Values[idx]
    # print('GammaValue(Idx) with mu_gamma closest to q_3^*', ValueClose)
    #
    #
    #
    # # Find Indices where the difference between the next one is larger than epsilon...
    # jump_idx = []
    # jump_xValues = []
    # jump_yValues = []
    # tmp = X_Values[0]
    # for idx, x in enumerate(X_Values):
    #     print(idx, x)
    #     if idx > 0:
    #         if abs(Y_Values[idx]-Y_Values[idx-1]) > 1:
    #             print('jump candidate')
    #             jump_idx.append(idx)
    #             jump_xValues.append(x)
    #             jump_yValues.append(Y_Values[idx])
    #
    #
    #
    #
    #
    #
    #
    # print("Jump Indices", jump_idx)
    # print("Jump X-values:", jump_xValues)
    # print("Jump Y-values:", jump_yValues)
    #
    # y_plotValues = [Y_Values[0]]
    # x_plotValues = [X_Values[0]]
    # # y_plotValues.extend(jump_yValues)
    # for i in jump_idx:
    #     y_plotValues.extend([Y_Values[i-1], Y_Values[i]])
    #     x_plotValues.extend([X_Values[i-1], X_Values[i]])
    #
    #
    # y_plotValues.append(Y_Values[-1])
    # # x_plotValues = [X_Values[0]]
    # # x_plotValues.extend(jump_xValues)
    # x_plotValues.append(X_Values[-1])
    #
    #
    # print("y_plotValues:", y_plotValues)
    # print("x_plotValues:", x_plotValues)


    # Y_Values[np.diff(y) >= 0.5] = np.nan


    #get values bigger than jump position
    # gamma = infty
    # x_rest = X_Values[X_Values>x_plotValues[1]]
    # Y_Values = np.array(Y_Values)  #convert the np array
    # y_rest = Y_Values[X_Values>x_plotValues[1]]
    #
    #
    # # gamma = 0
    # x_rest = X_Values[X_Values>x_plotValues[3]]
    # Y_Values = np.array(Y_Values)  #convert the np array
    # y_rest = Y_Values[X_Values>x_plotValues[3]]

    # gamma between
    # Y_Values = np.array(Y_Values)  #convert the np array
    # X_Values = np.array(X_Values)  #convert the np array
    #
    # x_one = X_Values[X_Values>x_plotValues[3]]
    # # ax.scatter(X_Values, Y_Values)
    # y_rest = Y_Values[X_Values>x_plotValues[3]]
    # ax.plot(X_Values[X_Values>0.135], Y_Values[X_Values<0.135])
    #
    #
    #


    # y_rest = Y_Values[np.nonzero(X_Values>x_plotValues[1]]
    # print('X_Values:', X_Values)
    # print('Y_Values:', Y_Values)
    # print('x_rest:', x_rest)
    # print('y_rest:', y_rest)
    # print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) )




    # --- Convert to numpy array
    Y_Values = np.array(Y_Values)
    X_Values = np.array(X_Values)

    Angle_Theta01 = np.array(Angle_Theta01)
    Angle_Theta025 = np.array(Angle_Theta025)

    Angle_Theta05 = np.array(Angle_Theta05)

    Angle_Theta075 = np.array(Angle_Theta075)
    Angle_Theta09 = np.array(Angle_Theta09)
    # ---------------- Create Plot -------------------

    # mpl.rcParams['text.usetex'] = True
    # mpl.rcParams["font.family"] = "serif"
    # mpl.rcParams["font.size"] = "9"

    # 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})
    mpl.rcParams['axes.labelpad'] = 0






    #---- Scale Figure apropriately to fit tex-File Width
    # width = 452.9679

    # width as measured in inkscape
    width = 6.28 *0.5
    width = 6.28 *0.333
    # width = 6.28
    height = width / 1.618

    #setup canvas first
    fig = plt.figure()      #main
    # fig, ax = plt.subplots()
    # fig, (ax, ax2) = plt.subplots(ncols=2)
    # fig,axes = plt.subplots(nrows=1,ncols=2,figsize=(width,height)) # more than one plot


    # fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97)  #TEST


    # TEST
    # mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1)
    # fig = plt.figure(figsize=(width+0.1,height+0.1))


    # mpl.rcParams['figure.figsize'] = (width,height)
    # fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
    # fig = plt.figure(figsize=(width,height)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
    # fig = plt.figure(figsize=set_size(width))
    # fig = plt.subplots(1, 1, figsize=set_size(width))

    # --- To create a figure half the width of your document:#
    # fig = plt.figure(figsize=set_size(width, fraction=0.5))



    #--- You must select the correct size of the plot in advance
    # fig.set_size_inches(3.54,3.54)

    # ax = plt.axes((0.15,0.18,0.8,0.8))
    # ax = plt.axes((0.15,0.18,0.6,0.6))
    # ax = plt.axes((0.15,0.2,0.75,0.75))
    # ax = plt.axes((0.18,0.2,0.75,0.75))


    # ax = plt.axes((0.28,0.3,0.65,0.65))  # This one!

    # ax = plt.axes((0.25,0.25,0.6,0.6))
    ax = plt.axes((0.25,0.28,0.6,0.6))
    # ax = plt.axes((0.1,0.1,0.5,0.8))
    # ax = plt.axes((0.1,0.1,1,1))
    # ax = plt.axes()

    # ax.spines['right'].set_visible(False)
    # ax.spines['left'].set_visible(False)
    # ax.spines['bottom'].set_visible(False)
    # ax.spines['top'].set_visible(False)
    # ax.tick_params(axis='x',which='major',direction='out',length=10,width=5,color='red',pad=15,labelsize=15,labelcolor='green',
    #                labelrotation=15)
    # ax.tick_params(axis='x',which='major', direction='out',pad=5,labelsize=10)
    # ax.tick_params(axis='y',which='major', length=5, width=1, direction='out',pad=5,labelsize=10)
    ax.tick_params(axis='x',which='major', direction='out',pad=3)
    ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)


    if gamma == '0':
        Title = r'$0< \gamma \ll 1$'
        # ax.set_title(Title)

    elif gamma == 'infinity':
        print('THIS CASE')
        Title = r'$\gamma \gg 1$'
        # ax.set_title(Title)

    else:
        Title = r'$ \gamma =$' + str(gamma)

    ax.set_title(Title,pad=3)

    # ax.xaxis.set_major_locator(MultipleLocator(0.05))
    # ax.xaxis.set_minor_locator(MultipleLocator(0.025))
    # ax.xaxis.set_major_locator(MultipleLocator(0.1))
    # ax.xaxis.set_minor_locator(MultipleLocator(0.05))
    ax.xaxis.set_major_locator(MultipleLocator(0.25))
    ax.xaxis.set_minor_locator(MultipleLocator(0.125))
    # ax.xaxis.set_major_locator(MultipleLocator(0.5))
    # ax.xaxis.set_minor_locator(MultipleLocator(0.25))
    #---- print data-types
    print(ax.xaxis.get_major_locator())
    print(ax.xaxis.get_minor_locator())
    print(ax.xaxis.get_major_formatter())
    print(ax.xaxis.get_minor_formatter())

    #---- Hide Ticks or Labels
    # ax.yaxis.set_major_locator(plt.NullLocator())
    # ax.xaxis.set_major_formatter(plt.NullFormatter())

    #---- Reducing or Increasing the Number of Ticks
    # ax.xaxis.set_major_locator(plt.MaxNLocator(3))
    # ax.yaxis.set_major_locator(plt.MaxNLocator(3))


    #----- Fancy Tick Formats
    ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
    ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
    ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))







    # --- manually change ticks&labels:
    # ax.set_xticks([0.2,1])
    # ax.set_xticklabels(['pos1','pos2'])

    # ax.set_yticks([0, np.pi/8, np.pi/4 ])
    # labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
    # ax.set_yticklabels(labels)

    a=ax.yaxis.get_major_locator()
    b=ax.yaxis.get_major_formatter()
    c = ax.get_xticks()
    d = ax.get_xticklabels()
    print('xticks:',c)
    print('xticklabels:',d)

    ax.grid(True,which='major',axis='both',alpha=0.3)






    # plt.figure()

    # f,ax=plt.subplots(1)

    # plt.title(r''+ yName + '-Plot')
    # plt.plot(X_Values, Y_Values,linewidth=2, '.k')
    # plt.plot(X_Values, Y_Values,'.k',markersize=1)
    # plt.plot(X_Values, Y_Values,'.',markersize=0.8)

    # plt.plot(X_Values, Y_Values)

    # ax.plot([[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])



    # Gamma = '0'
    # ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b')
    #
    # ax.plot([x_plotValues[1],x_plotValues[3]], [y_plotValues[2],y_plotValues[3]] , 'b')
    #
    # ax.plot(x_rest, y_rest, 'b')


    # Gamma between

    # x jump values (gamma 0): [0.13606060606060608, 0.21090909090909093]

    # ax.plot([[0,jump_xValues[0]], [0, 0]] , 'b')
    # ax.plot([jump_xValues[0],xmin], [y_plotValues[2],y_plotValues[2]] , 'b')

    # ax.plot([[0,0.13606060606060608], [0, 0]] , 'b')
    # ax.plot([[0.13606060606060608,xmin], [(math.pi/2),(math.pi/2)]], 'b')

    # jump_xValues[0]



    # --- leave out jumps:
    # ax.scatter(X_Values, Y_Values)

    # ax.set_xlabel(r"prestrain ratio $\theta_\rho$",labelpad=0)
    ax.set_xlabel(r"$\theta_\rho$",labelpad=0)
    ax.set_ylabel(r"$\alpha$")
    # ax.set_ylabel(r"angle $\alpha$")

    # ax.set_title(r"$\gamma = \ $"+str(gamma), fontsize=10)


    if Jumps:

        # --- leave out jumps:
        if gamma == 'infinity':
            ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]] , 'royalblue')
            ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue')



            # ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]])
            # ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]])




        # ax.plot(X_Values[X_Values>0.136], Y_Values[X_Values>0.136])
        # ax.plot(X_Values[X_Values<0.135], Y_Values[X_Values<0.135])
        # ax.scatter(X_Values, Y_Values)
        # ax.plot(X_Values, Y_Values)

        # plt.plot(x_plotValues, y_plotValues,'.')
        # plt.scatter(X_Values, Y_Values, alpha=0.3)
        # plt.scatter(X_Values, Y_Values)
        # plt.plot(X_Values, Y_Values,'.')
        # plt.plot([X_Values[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]])
        # plt.axis([0, 6, 0, 20])

        # ax.set_xlabel(r"volume fraction $\theta$", size=11)
        # ax.set_ylabel(r"angle $\angle$",  size=11)
        # ax.set_xlabel(r"volume fraction $\theta$")
        # # ax.set_ylabel(r"angle $\angle$")
        # ax.set_ylabel(r"angle $\alpha$")
        # plt.ylabel('$\kappa$')

        # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
        # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))




        # Plot every other line.. not the jumps...

        if gamma == '0':
            tmp = 1
            for idx, x in enumerate(x_plotValues):
                if idx > 0 and tmp == 1:
                    # plt.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]] )
                    ax.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]], 'royalblue', zorder=2)
                    tmp = 0
                else:
                    tmp = 1

        # plt.plot([x_plotValues[0],x_plotValues[1]] ,[y_plotValues[0],y_plotValues[1]] )
        # plt.plot([x_plotValues[2],x_plotValues[3]] ,[y_plotValues[2],y_plotValues[3]] )
        # plt.plot([x_plotValues[4],x_plotValues[5]] ,[y_plotValues[4],y_plotValues[5]] )
        # plt.plot([x_plotValues[6],x_plotValues[7]] ,[y_plotValues[6],y_plotValues[7]] )


        for x in jump_xValues:
            plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1, zorder=1)
            # plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed',  label=r'$\theta_*$')

        # plt.axvline(x_plotValues[1],ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')

        # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
        # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
        # plt.legend()


        # -- SETUP LEGEND
        # ax.legend(prop={'size': 11})
        # ax.legend()

        # ------------------ SAVE FIGURE
        # tikzplotlib.save("TesTout.tex")
        # plt.close()
        # mpl.rcParams.update(mpl.rcParamsDefault)

        # plt.savefig("graph.pdf",
        #             #This is simple recomendation for publication plots
        #             dpi=1000,
        #             # Plot will be occupy a maximum of available space
        #             bbox_inches='tight',
        #             )
        # plt.savefig("graph.pdf")



        # ---- ADD additional scatter:
        # ax.scatter(X_Values,Y_Values,s=1,c='black',zorder=4)

        # Find transition point
        lastIdx = len(Y_Values)-1

        for idx, y in enumerate(Y_Values):
            if idx != lastIdx:
                if abs(y-0) < 0.01 and abs(Y_Values[idx+1] - 0) > 0.05:
                    transition_point1 = X_Values[idx+1]
                    print('transition point1:', transition_point1 )
                if abs(y-0.5*np.pi) < 0.01 and abs(Y_Values[idx+1] -0.5*np.pi)>0.01:
                    transition_point2 = X_Values[idx]
                    print('transition point2:', transition_point2 )
                if abs(y-0) > 0.01 and abs(Y_Values[idx+1] - 0) < 0.01:
                    transition_point3 = X_Values[idx+1]
                    print('transition point3:', transition_point3 )

        # Add transition Points:
        if gamma == '0':
            ax.scatter([transition_point1, transition_point2],[np.pi/2,np.pi/2],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
                                      edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)

            ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
                               )

            ax.text(transition_point2+0.012 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
                               )
        else:
            ax.scatter([transition_point1, transition_point2, transition_point3 ],[np.pi/2,np.pi/2,0 ],s=6, marker='o', cmap=None, norm=None, facecolor = 'black',
                                      edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)

            ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
                               )

            ax.text(transition_point2 +0.011 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
                               )

            ax.text(transition_point3 +0.009 , 0+0.08, r"$3$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5)
                                   )

    else:
            # ax.scatter(X_Values,Y_Values,s=1, marker='o', cmap=None, norm=None, facecolor = 'blue',
            #                           edgecolor = 'none', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
             # ---------------------------------------------------------------
            # l1 = ax.scatter(X_Values,Angle_alpha0,s=1, marker='o', edgecolor = 'black',cmap=None, norm=None, vmin=None, vmax=None, alpha=0.75, linewidths=None, zorder=4)
            # l6 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1, label=r"$\theta_\rho = -1.0$")
            # l4 = ax.scatter(X_Values,Angle_alphaNeg05,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l3 = ax.scatter(X_Values,Angle_alphaNeg025,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
            # # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
            # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l2 = ax.scatter(X_Values,Angle_alphaNeg0125,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            #
            # line_labels = [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = 3.0$"]
            # ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8  ])
            # labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$',r'$5\pi/8$']
            # ax.set_yticklabels(labels)
            #
            # ax.legend(handles=[l1,l2,l3,l4, l5, l6, l7],
            #           labels= [ r"$\theta_\rho = 0$", r"$\theta_\rho = -0.125$", r"$\theta_\rho = -0.25$", r"$\theta_\rho = -0.5$", r"$\theta_\rho = -0.75$",  r"$\theta_\rho = -1.0$",  r"$\theta_\rho = 3.0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))
           # ---------------------------------------------------------------
            # l1 = ax.scatter(X_Values,Angle_alphaNeg1,s=2, marker='s', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=1)
            # l2 = ax.scatter(X_Values,Angle_alphaNeg0875,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
            # l3 = ax.scatter(X_Values,Angle_alphaNeg075,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l4 = ax.scatter(X_Values,Angle_alphaNeg0625,s=2, marker='o',cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)
            # l5 = ax.scatter(X_Values,Angle_alphaNeg05,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l6 = ax.scatter(X_Values,Angle_alphaNeg025,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l7 = ax.scatter(X_Values,Angle_alphaNeg0125,s=2, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)
            # l8 = ax.scatter(X_Values,Angle_alpha0,s=2, marker='s', edgecolor = 'black', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=4)

            # l1 = ax.plot(X_Values,Angle_alphaNeg05, color='blue', linewidth=1.5, zorder=3, label=r"$\theta_\rho=-0.5$")
            # l2 = ax.plot(X_Values,Angle_alphaNeg055, linewidth=1.5, linestyle = '--', zorder=3,label=r"$\theta_\rho=-0.55$")
            # l3 = ax.plot(X_Values,Angle_alphaNeg06,color='orangered', linewidth=1.5 ,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.6$")
            # l4 = ax.plot(X_Values,Angle_alphaNeg065, linewidth=1.5 ,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.65$")
            # l5 = ax.plot(X_Values,Angle_alphaNeg07,color='orange', linewidth=1.5 ,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.7$")
            # l6 = ax.plot(X_Values,Angle_alphaNeg075, linewidth=1.5,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.75$")
            # l7 = ax.plot(X_Values,Angle_alphaNeg08, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.8$")
            # l8 = ax.plot(X_Values,Angle_alphaNeg085, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.85$")
            # l9 = ax.plot(X_Values,Angle_alphaNeg09, color='teal',linestyle = '--', linewidth=1.5 ,  zorder=3, label=r"$\theta_\rho=-0.9$")
            # l10 = ax.plot(X_Values,Angle_alphaNeg095, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.95$")
            # l11 = ax.plot(X_Values,Angle_alphaNeg1, color='red', linewidth=1.5 ,zorder=1, label=r"$\theta_\rho=-1.0$")


            # l1 = ax.plot(X_Values,Angle_Theta01, color='blue', linewidth=1.5, zorder=3, label=r"$\theta=0.1$")
            # # l2 = ax.plot(X_Values,Angle_alphaNeg055, linewidth=1.5, linestyle = '--', zorder=3,label=r"$\theta_\rho=-0.55$")
            # l3 = ax.plot(X_Values,Angle_Theta025,color='orangered', linewidth=1.5  ,zorder=3, label=r"$\theta = 0.25$")
            # # l4 = ax.plot(X_Values,Angle_alphaNeg065, linewidth=1.5 ,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.65$")
            # l5 = ax.plot(X_Values,Angle_Theta05,color='orange', linewidth=1.5  ,zorder=3, label=r"$\theta = 0.5$")
            # # l6 = ax.plot(X_Values,Angle_alphaNeg075, linewidth=1.5,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.75$")
            # l7 = ax.plot(X_Values,Angle_Theta075, linewidth=1.5 ,  zorder=3, label=r"$\theta = 0.75$")
            # # l8 = ax.plot(X_Values,Angle_alphaNeg085, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.85$")
            # l9 = ax.plot(X_Values,Angle_Theta09, color='teal', linewidth=1.5 ,  zorder=3, label=r"$\theta =0.9$")
            # # l10 = ax.plot(X_Values,Angle_alphaNeg095, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.95$")


            # l1 = ax.scatter(X_Values,Angle_Theta01, color='blue', s=2, zorder=3, label=r"$\theta=0.1$")
            # l2 = ax.plot(X_Values,Angle_alphaNeg055, linewidth=1.5, linestyle = '--', zorder=3,label=r"$\theta_\rho=-0.55$")
            # l3 = ax.scatter(X_Values,Angle_Theta025,color='orangered', s=2  ,zorder=3, label=r"$\theta = 0.25$")
            # l4 = ax.plot(X_Values,Angle_alphaNeg065, linewidth=1.5 ,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.65$")
            l5 = ax.scatter(X_Values,Angle_Theta05, color='royalblue', s=0.15 ,zorder=3, label=r"$\theta = 0.5$")
            # l5 = ax.scatter(X_Values,Angle_Theta05, color='royalblue' ,zorder=3, label=r"$\theta = 0.5$")
            # l6 = ax.plot(X_Values,Angle_alphaNeg075, linewidth=1.5,linestyle = '--' ,zorder=3, label=r"$\theta_\rho=-0.75$")
            # l7 = ax.scatter(X_Values,Angle_Theta075, s=2, zorder=3, label=r"$\theta = 0.75$")
            # l8 = ax.plot(X_Values,Angle_alphaNeg085, linewidth=1.5,linestyle = '--' ,  zorder=3, label=r"$\theta_\rho=-0.85$")
            # l9 = ax.scatter(X_Values,Angle_Theta09, color='teal',s=2,  zorder=3, alpha=1.0, label=r"$\theta =0.9$")

            # ax.legend(handles=[l1[0],l2[0],l3[0],l4[0], l5[0], l6[0], l7[0], l8[0], l9[0], l10[0], l11[0]],
            #           # labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))

            # ax.legend(handles=[l1[0],l3[0], l5[0], l7[0], l9[0]] ,
            #           # labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))

            # ax.legend(handles=[l1,l3, l5, l7, l9] ,
            #           # labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))

            # ax.legend(handles=[l5] ,
            #           # labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))

    # ax.plot(X_Values, Y_Values,   marker='o',  markerfacecolor='orange', markeredgecolor='black', markeredgewidth=1,  linewidth=1, zorder=3)
            # l7 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', facecolor = 'none',edgecolor = 'forestgreen', cmap=None, norm=None, vmin=None, vmax=None, alpha=1.0, linewidths=None, zorder=5)
            # l4 = ax.scatter(X_Values,Angle_alpha3,s=1, marker='o', markerfacecolor='red',markeredgecolor='black',markeredgewidth=2, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.5, linewidths=None, zorder=3)
            # l5 = ax.scatter(X_Values,Angle_alphaNeg075,s=1, marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3)


            # line_labels = [r"$\theta_\rho = -1.0$",r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -  \frac{3}{4}$", r"$\theta_\rho = -  \frac{5}{8}$",r"$\theta_\rho = - 0.5 $" ,r"$\theta_\rho = -  0.25", r"$\theta_\rho = -  \frac{1}{8}" , r"$\theta_\rho = 0$"]
            ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8  ])
            ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2])
            labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$',r'$5\pi/8$']
            labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']

            ax.set_ylim(0,np.pi/2+0.05)
            ax.set_yticklabels(labels)



            # ax.legend(handles=[l1,l2,l3,l4, l5, l6, l7, l8],
            #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))

            #
            # ax.legend(handles=[l1,l3, l5, l7, l8],
            #           labels= [r"$\theta_\rho = -1.0$", r"$\theta_\rho = -  \frac{7}{8}$", r"$\theta_\rho = -\frac{3}{4}$" , r"$\theta_\rho = -  \frac{5}{8}$", r"$\theta_\rho = - \frac{1}{2} $" , r"$\theta_\rho = - \frac{1}{4}$", r"$\theta_\rho = -  \frac{1}{8}$" , r"$\theta_\rho = 0$"],
            #           loc='upper left',
            #           bbox_to_anchor=(1,1))
            #





            # fig.legend([l1, l2, l3, l4],     # The line objects
            #            labels=line_labels,   # The labels for each line
            #            # loc="upper center",   # Position of legend
            #            loc='upperleft', bbox_to_anchor=(1,1),
            #            borderaxespad=0.15    # Small spacing around legend box
            #            # title="Legend Title"  # Title for the legend
            #            )

    pdf_outputName = 'Plot-Angle-Alpha_Gamma'+ str(gamma)+ '_transition'+'.pdf'

    fig.set_size_inches(width, height)
    # fig.savefig('Plot-Angle-Theta.pdf')
    fig.savefig(pdf_outputName)




    # tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm')

    # tikz_save('fig.tikz',
    #            figureheight = '\\figureheight',
    #            figurewidth = '\\figurewidth')

    # ----------------------------------------


# plt.show()
# #---------------------------------------------------------------