Skip to content
Snippets Groups Projects
Plot-AngleCurv-Alpha_intermediateGamma_SubPlots.py 8.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • 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
    
    
    
    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
    
    
    
    # ----- Setup Paths -----
    # InputFile  = "/inputs/cellsolver.parset"
    # OutputFile = "/outputs/output.txt"
    
    InputFile  = "/inputs/computeMuGamma.parset"
    OutputFile = "/outputs/outputMuGamma.txt"
    
    # --------- Run  from src folder:
    path_parent = os.path.dirname(os.getcwd())
    os.chdir(path_parent)
    path = os.getcwd()
    print(path)
    InputFilePath = os.getcwd()+InputFile
    OutputFilePath = os.getcwd()+OutputFile
    print("InputFilepath: ", InputFilePath)
    print("OutputFilepath: ", OutputFilePath)
    print("Path: ", path)
    
    #---------------------------------------------------------------
    print('---- Input parameters: -----')
    mu1 = 1.0  #10.0
    rho1 = 1.0
    lambda1 = 0.0
    
    # # INTERESTING!:
    # alpha = 3
    beta = 2.0
    # theta= 1/8
    
    # alpha = 2.0
    theta = 0.5
    beta = 10.0
    #TEST
    # beta=2
    
    plot_curvature = True
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
    plot_curvature = False
    
    
    
    
    gamma = 'infinity'  #Elliptic Setting
    gamma = '0'       #Hyperbolic Setting
    
    Gamma_Values = [0.5, 0.75, 1.5, 3.0]
    Gamma_Values = [3.0]
    Gamma_Values = ['0', 'infinity','infinity','infinity','infinity','infinity']
    Gamma_Values = ['0', 0.5, 0.75, 1.5,3.0, 'infinity']
    Gamma_Values = ['0', 0.5, 0.75, 1.0, 1.5, 'infinity']
    # Gamma_Values = ['0',0.25, 0.5, 0.75, 1.5, 'infinity']
    print('(Input) Gamma_Values:', Gamma_Values)
    
    # --- 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
    
    
    # xmin = -2.0
    # xmin = -1.0
    # xmax = 1.0
    #
    # xmin = -1.25
    xmin = -1.5
    xmax = 1.0
    
    xmin = -1.25
    xmax = 1.25
    
    numPoints = 2000
    numPoints = 300
    numPoints = 500
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
    numPoints = 50
    numPoints = 100
    
    X_Values = np.linspace(xmin, xmax, num=numPoints)
    X_Values = np.array(X_Values)
    print(X_Values)
    
    Angle_Container = []
    
    for gamma in Gamma_Values:  # COMPUTE DATA
    
        print('mu1: ', mu1)
        print('rho1: ', rho1)
        # print('alpha: ', alpha)
        print('beta: ', beta)
        print('theta: ', theta)
        print('gamma:', gamma)
        print('----------------------------')
    
    
    
        tmp_Container = []
    
    
        for alpha in X_Values:
            print('Situation of Lemma1.4')
            q12 = 0.0
    
            q3 = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath)
    
            G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3,  mu1, rho1)
            if plot_curvature:
                tmp_Container.append(curvature)
            else:
                tmp_Container.append(angle)
    
    
    
    
    
        tmp_Container = np.array(tmp_Container)  #convert the np array
        Angle_Container.append(tmp_Container)
    
    
    
    
    
    ##########################################################
    # ---------------- CREATE PLOT  -------------------
    # Styling
    plt.style.use("seaborn-darkgrid")
    # plt.style.use("seaborn-whitegrid")
    plt.style.use("seaborn")
    
    # 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'] = 0
    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()
    
    # ax = plt.axes((0.25,0.28,0.6,0.6))
    #
    # 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)
    
    
    gs = fig.add_gridspec(nrows=2,ncols=3, hspace=0.25, wspace=0.25)
    
    ax4 = fig.add_subplot(gs[1, 0])
    ax5 = fig.add_subplot(gs[1, 1],sharey=ax4)
    ax6 = fig.add_subplot(gs[1, 2],sharey=ax4)
    plt.setp(ax5.get_yticklabels(), visible=False)
    plt.setp(ax6.get_yticklabels(), visible=False)
    
    ax1 = fig.add_subplot(gs[0, 0],sharex=ax4)
    ax2 = fig.add_subplot(gs[0, 1],sharey=ax1)
    ax3 = fig.add_subplot(gs[0, 2],sharey=ax1)
    plt.setp(ax1.get_xticklabels(), visible=False)
    plt.setp(ax2.get_xticklabels(), visible=False)
    plt.setp(ax3.get_xticklabels(), visible=False)
    plt.setp(ax2.get_yticklabels(), visible=False)
    plt.setp(ax3.get_yticklabels(), visible=False)
    
    
    
    
    ax = [ax1, ax2,ax3,ax4,ax5,ax6]
    
    print('ax:', ax)
    print('ax[0]:', ax[0])
    for i in range(len(ax)):
        print('i=',i)
    
        gamma = Gamma_Values[i]
        # TITLE
        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[i].set_title(Title,pad=5)
    
        # ax[i].xaxis.set_major_locator(MultipleLocator(0.25))
        # ax[i].xaxis.set_minor_locator(MultipleLocator(0.125))
        ax[i].xaxis.set_major_locator(MultipleLocator(1.0))
        ax[i].xaxis.set_minor_locator(MultipleLocator(0.5))
        # ax[i].xaxis.set_major_locator(MultipleLocator(0.5))
        # ax[i].xaxis.set_minor_locator(MultipleLocator(0.25))
    
        #---- print data-types
        print(ax[i].xaxis.get_major_locator())
        print(ax[i].xaxis.get_minor_locator())
        print(ax[i].xaxis.get_major_formatter())
        print(ax[i].xaxis.get_minor_formatter())
    
    
        #----- Fancy Tick Formats
        if plot_curvature == False:
            ax[i].yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4))
            ax[i].yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
            ax[i].yaxis.set_major_formatter(plt.FuncFormatter(format_func))
    
        a=ax[i].yaxis.get_major_locator()
        b=ax[i].yaxis.get_major_formatter()
        c = ax[i].get_xticks()
        d = ax[i].get_xticklabels()
    
    
    
        # ax[i].set_xticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2])
    
    
        # print('xticks:',c)
        # print('xticklabels:',d)
    
        ax[i].grid(True,which='major',axis='both',alpha=0.3)
    
    
    
        # ax[i].set_xlabel(r"$\theta_\rho$",labelpad=0)
        # ax[i].set_ylabel(r"$\alpha$")
        # ax.set_ylabel(r"angle $\alpha$")
        if plot_curvature:
            line= ax[i].scatter(X_Values,Angle_Container[i], color='forestgreen', s=0.15 ,zorder=3, label=r"$\theta = 0.5$")
        else:
            line= ax[i].scatter(X_Values,Angle_Container[i], color='royalblue', s=0.15 ,zorder=3, label=r"$\theta = 0.5$")
    
        if plot_curvature == False:
            # ax[i].set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8  ])
            ax[i].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[i].set_ylim(0-0.1,np.pi/2+0.1)
            ax[i].set_yticklabels(labels)
    
    plt.subplots_adjust(wspace=1, hspace=1)
    fig.subplots_adjust(bottom=0.15)
    # fig.tight_layout()
    fig.text(0.5, 0.04, r"$\theta_\rho$", ha='center')
    
    
    if plot_curvature:
        fig.text(0.03, 0.5, r"curvature $\kappa$", va='center', rotation='vertical')
        pdf_outputName = 'Plot-Curv-Alpha_Gamma'+ str(gamma)+ '_transition'+'.pdf'
    else:
        fig.text(0.03, 0.5, r"angle $\alpha$", va='center', rotation='vertical')
        pdf_outputName = 'Plot-Angle-Alpha_Gamma'+ str(gamma)+ '_transition'+'.pdf'
    
    fig.set_size_inches(width, height)
    fig.savefig(pdf_outputName)
    plt.show()