Add Script to plot cylindrical Minimizers!

src/CylindricalMinimizer-Plot.py

+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 matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+from mpl_toolkits.mplot3d import Axes3D
+import pandas as pd
+import matplotlib.colors as mcolors
+from matplotlib import cm
+
+from mpl_toolkits.mplot3d.proj3d import proj_transform
+# from mpl_toolkits.mplot3d.axes3d import Axes3D
+from matplotlib.text import Annotation
+from matplotlib.patches import FancyArrowPatch
+
+# Extra packages :
+# from HelperFunctions import *
+# from ClassifyMin import *
+# from subprocess import Popen, PIPE
+#import sys
+
+###################### Documentation #########################
+
+#..... add description here
+
+###########################################################
+
+
+def rot(v,alpha):
+
+#rotate about axis v with degree deg in radians:
+
+
+    tmp = np.array([ [v[0]**2*(1-np.cos(alpha))+np.cos(alpha), v[0]*v[1]*(1-np.cos(alpha))-v[2]*np.sin(alpha), v[0]*v[2]*(1-np.cos(alpha))+ v[1]*np.sin(alpha) ],
+         [v[0]*v[1]*(1-np.cos(alpha))+v[2]*np.sin(alpha), v[1]**2*(1-np.cos(alpha))+np.cos(alpha), v[1]*v[2]*(1-np.cos(alpha))+v[0]*np.sin(alpha) ],
+         [v[2]*v[0]*(1-np.cos(alpha))-v[1]*np.sin(alpha), v[2]*v[1]*(1-np.cos(alpha))+v[0]*np.sin(alpha) , v[2]**2*(1-np.cos(alpha))+np.cos(alpha) ] ])
+
+    return tmp
+
+
+
+def rotate_data(X, R):
+#rotate about axis v with degree deg in radians:
+# X : DataSet
+# R : RotationMatrix
+    print('ROTATE DATA FUNCTION ---------------')
+
+    rot_matrix = R
+    # print('rot_matrix:', rot_matrix)
+    # print('rot_matrix.shape:', rot_matrix.shape)
+    # print('X', X)
+    # print('shape of X[0]', X.shape[0])
+    B = np.dot(rot_matrix, X.reshape(rot_matrix.shape[1],-1))
+    # print('shape of B', B.shape)
+    # print('B',B)
+    # print('B[0,:]', B[0,:])
+    # print('B[0,:].shape', B[0,:].shape)
+    Out = np.array([B[0,:].reshape(X.shape[1],X.shape[2]), B[1,:].reshape(X.shape[1],X.shape[2]), B[2,:].reshape(X.shape[1],X.shape[2])])
+    print('shape of Out', Out.shape)
+
+    return Out
+
+# def rotate_data(X, v,alpha): #(Old Version)
+# #rotate about axis v with degree deg in radians:
+# # X : DataSet
+#     print('ROTATE DATA FUNCTION ---------------')
+#     # v = np.array([1,0,0])
+#     # rotM = rot(v,np.pi/2)
+#     # print('rotM:', rotM)
+#     rot_matrix = rot(v,alpha)
+#     # print('rot_matrix:', rot_matrix)
+#     # print('rot_matrix.shape:', rot_matrix.shape)
+#
+#     # print('X', X)
+#     # print('shape of X[0]', X.shape[0])
+#     B = np.dot(rot_matrix, X.reshape(rot_matrix.shape[1],-1))
+#
+#     # print('shape of B', B.shape)
+#     # print('B',B)
+#     # print('B[0,:]', B[0,:])
+#     # print('B[0,:].shape', B[0,:].shape)
+#     Out = np.array([B[0,:].reshape(X.shape[1],X.shape[2]), B[1,:].reshape(X.shape[1],X.shape[2]), B[2,:].reshape(X.shape[1],X.shape[2])])
+#     print('shape of Out', Out.shape)
+#
+#     return Out
+
+
+# def translate_data(X, v):  ...
+# #rotate about axis v with degree deg in radians:
+# # X : DataSet
+#     print('ROTATE DATA FUNCTION ---------------')
+#     # v = np.array([1,0,0])
+#     # rotM = rot(v,np.pi/2)
+#     # print('rotM:', rotM)
+#
+#     print('X', X)
+#     print('shape of X[0]', X.shape[0])
+#
+#     Out = X + v
+#     return Out
+
+
+def u(x,kappa,e):
+
+    tmp = (x.dot(e))*kappa
+    # print('tmp for u',tmp)
+    if kappa == 0 :
+        tmp = np.array([0*x[0],  x[0]*e[0] + x[1]*e[1], x[1]*e[0] - x[0]*e[1] ])
+    else :
+        tmp = np.array([-(1/kappa)*np.cos(tmp)+(1/kappa),  (1/kappa)*np.sin(tmp), -x[0]*e[1]+x[1]*e[0] ])
+    return tmp
+
+
+
+
+def grad_u(x,kappa,e):
+
+    tmp = (x.dot(e))*kappa
+    # print('tmp',tmp)
+
+    grad_u = np.array([ [np.sin(tmp)*e[0], np.sin(tmp)*e[1]], [np.cos(tmp)*e[0], np.cos(tmp)*e[1]], [-e[1], e[0]] ])
+    # print('produkt', grad_u.dot(e) )
+    mapped_e = grad_u.dot(e)
+    # print('mapped_e:', mapped_e)
+    # print('siize of mapped_e', mapped_e.shape)
+    # mapped_e = mapped_e.transpose()
+    # print('mapped_e:', mapped_e)
+    # print('siize of mapped_e', mapped_e.shape)
+    return mapped_e
+
+def compute_normal(x,kappa,e):
+    tmp = (x.dot(e))*kappa
+    partial1_u = np.array([ np.sin(tmp)*e[0] ,np.cos(tmp)*e[0], -e[1] ])
+    partial2_u = np.array([ np.sin(tmp)*e[1], np.cos(tmp)*e[1], e[0]      ])
+    normal = np.cross(partial1_u,partial2_u)
+    # print('normal=',normal)
+    return normal
+
+
+
+class Annotation3D(Annotation):
+    def __init__(self, text, xyz, *args, **kwargs):
+        super().__init__(text, xy=(0, 0), *args, **kwargs)
+        self._xyz = xyz
+
+    def draw(self, renderer):
+        x2, y2, z2 = proj_transform(*self._xyz, self.axes.M)
+        self.xy = (x2, y2)
+        super().draw(renderer)
+
+def _annotate3D(ax, text, xyz, *args, **kwargs):
+    '''Add anotation `text` to an `Axes3d` instance.'''
+
+    annotation = Annotation3D(text, xyz, *args, **kwargs)
+    ax.add_artist(annotation)
+
+setattr(Axes3D, 'annotate3D', _annotate3D)
+
+class Arrow3D(FancyArrowPatch):
+
+    def __init__(self, x, y, z, dx, dy, dz, *args, **kwargs):
+        super().__init__((0, 0), (0, 0), *args, **kwargs)
+        self._xyz = (x, y, z)
+        self._dxdydz = (dx, dy, dz)
+
+    def draw(self, renderer):
+        x1, y1, z1 = self._xyz
+        dx, dy, dz = self._dxdydz
+        x2, y2, z2 = (x1 + dx, y1 + dy, z1 + dz)
+
+        xs, ys, zs = proj_transform((x1, x2), (y1, y2), (z1, z2), self.axes.M)
+        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
+        super().draw(renderer)
+
+
+def _arrow3D(ax, x, y, z, dx, dy, dz, *args, **kwargs):
+    '''Add an 3d arrow to an `Axes3D` instance.'''
+
+    arrow = Arrow3D(x, y, z, dx, dy, dz, *args, **kwargs)
+    ax.add_artist(arrow)
+
+setattr(Axes3D, 'arrow3D', _arrow3D)
+################################################################################################################
+################################################################################################################
+################################################################################################################
+
+
+############################################################################################################################################
+####################################################################### KAPPA NEGATIVE ####################################################
+############################################################################################################################################
+kappa = -2
+num_Points = 100
+
+e = np.array([1,0])
+# e = np.array([0,1])
+# e = np.array([1/np.sqrt(2),1/np.sqrt(2)])
+# e = np.array([1/2,np.sqrt(3)/2])
+# e = np.array([np.sqrt(3)/2,1/2])
+# e = np.array([-1,0])
+# e = np.array([0,-1])
+
+# Creating dataset
+x = np.linspace(-1,1,num_Points)
+y = np.linspace(-1/2,1/2,num_Points)
+
+print('type of x', type(x))
+print('max of x:', max(x))
+print('max of y:', max(y))
+# print('x:', x)
+
+x1, x2 = np.meshgrid(x,y)
+zero = 0*x1
+
+if kappa == 0 :
+    u1 = 0*x1
+    u2 = x1*e[0] + x2*e[1]
+    u3 = x2*e[0] - x1*e[1]
+
+else :
+    u1 = -(1/kappa)*np.cos(kappa*(x1*e[0]+x2*e[1])) + (1/kappa)
+    u2 = (1/kappa)*np.sin(kappa*(x1*e[0]+x2*e[1]))
+    u3 = x2*e[0] -x1*e[1]
+
+# print('np.size(u1)',np.size(u1))
+# print('u1.shape',u1.shape)
+# colorfunction=(u1**2+u2**2)
+# print('colofunction',colorfunction)
+
+# print('u1.size:',np.size(u1))
+# tmp = np.ones(np.size(u1))*kappa
+# print('np.size(tmp)',np.size(tmp))
+B = np.full_like(u1, 1)
+# colorfunction=(u3)                                              # TODO Color by angle
+# colorfunction=(np.ones(np.size(u1))*kappa)
+colorfunction=(B*kappa)
+# print('colofunction',colorfunction)
+norm=mcolors.Normalize(colorfunction.min(),colorfunction.max())
+
+# Display the mesh
+fig = plt.figure()
+
+ax = plt.axes(projection ='3d', adjustable='box')
+
+
+###---TEST MAP e-vectprs!
+# e1 = np.array([1,0])
+# e2 = np.array([0,1])
+# e3 = np.array([1/np.sqrt(2),1/np.sqrt(2)])
+# e1 = np.array([0,1])
+# e2 = np.array([-1,0])
+# e3 = np.array([-1/np.sqrt(2),1/np.sqrt(2)])
+# e1_mapped = u(e1,kappa,e1)
+# e2_mapped = u(e2,kappa,e2)
+# e3_mapped = u(e3,kappa,e3)
+# print('e1 mapped:',e1_mapped)
+# print('e2 mapped:',e2_mapped)
+# print('e3 mapped:',e3_mapped)
+### -----------------------------------
+
+#--e1 :
+# Rotation_angle = -np.pi/2
+# Rotation_vector = np.array([0,1,0])
+
+#--e2:
+Rotation_angle = np.pi/2
+Rotation_vector = np.array([1,0,0])
+
+###--e = np.array([1/np.sqrt(2),1/np.sqrt(2)])
+# Rotation_angle = -np.pi/2
+# Rotation_vector = np.array([1,0,0])
+# #2te rotation :
+# Rotation_angle = np.pi/4
+# Rotation_vector = np.array([0,0,1])
+
+
+
+Rotation_angle = -np.pi/2
+Rotation_angle = 0
+# Rotation_angle = np.pi/2
+Rotation_vector = np.array([0,1,0])
+Rotation_vector = np.array([1,0,0])
+
+# rot(np.array([0,1,0]),np.pi/2)
+
+#### if e1 :
+Rotation = rot(np.array([0,1,0]),-np.pi/2)
+#### if e2:
+Rotation = rot(np.array([0,1,0]),-np.pi/2).dot(rot(np.array([1,0,0]),-np.pi/2))
+# #### if e3 :
+Rotation = rot(np.array([0,0,1]),np.pi/4).dot(rot(np.array([0,1,0]),-np.pi/2).dot(rot(np.array([1,0,0]),-np.pi/2)))
+Rotation = rot(np.array([0,0,1]),np.pi/2).dot(rot(np.array([0,1,0]),-np.pi/2).dot(rot(np.array([1,0,0]),-np.pi/2)))
+
+# Rotation = rot(np.array([1,0,0]),np.pi/2)
+
+# Rotation_vector = e3_mapped  #TEST
+# Rotation_vector = np.array([-1/np.sqrt(2),1/np.sqrt(2)])
+# Rotation_vector = np.array([0,0,1])
+
+# v = np.array([1,0,0])
+# X = np.array([u1,u2,u3])
+
+
+
+
+# T = rotate_data(np.array([u1,u2,u3]),Rotation_vector,Rotation_angle)
+T = rotate_data(np.array([u1,u2,u3]),Rotation)
+
+ax.plot_surface(T[0], T[1], T[2], color = 'w', rstride = 1, cstride = 1, facecolors=cm.brg(colorfunction), alpha=.4, zorder=4)
+
+
+
+###---- PLOT PARAMETER-PLANE:
+# ax.plot_surface(x1,x2,zero,color = 'w', rstride = 1, cstride = 1 )
+
+
+#midpoint:
+midpoint = np.array([(max(x)+min(x))/2,(max(y)+min(y))/2])
+print('midpoint',midpoint)
+
+# Map midpoint:
+midpoint_mapped = u(midpoint,kappa,e)
+print('mapped midpoint', midpoint_mapped)
+
+#map origin
+origin = np.array([0,0])
+origin_mapped = u(origin,kappa,e)
+
+
+mapped_e = grad_u(midpoint,kappa,e)
+normal = compute_normal(midpoint,kappa,e)
+
+
+rotM = rot(Rotation_vector,Rotation_angle)
+
+mapped_e = Rotation.dot(mapped_e)
+normal = Rotation.dot(normal)
+
+
+# Plot Mapped_midPoint
+ax.plot(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],    # data
+marker='o',     # each marker will be rendered as a circle
+markersize=4,   # marker size
+markerfacecolor='orange',   # marker facecolor
+markeredgecolor='black',  # marker edgecolor
+markeredgewidth=1,       # marker edge width
+linewidth=1,
+zorder=5)          # line width
+
+# ax.quiver([midpoint_mapped[0]], [midpoint_mapped[1]], [midpoint_mapped[2]], [mapped_e[0]], [mapped_e[1]], [mapped_e[2]], color="red")
+# ax.quiver([midpoint_mapped[0]], [midpoint_mapped[1]], [midpoint_mapped[2]], [normal[0]], [normal[1]], [normal[2]], color="blue")
+
+ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+           mapped_e[0],mapped_e[1],mapped_e[2],
+           mutation_scale=10,
+           arrowstyle="-|>",
+           linestyle='dashed',fc='green',
+           ec ='green',
+           zorder=5)
+
+ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+           normal[0],normal[1],normal[2],
+           mutation_scale=10,
+           arrowstyle="-|>",
+           linestyle='dashed',fc='blue',
+           ec ='blue',
+           zorder = 5)
+
+
+
+###-- TEST Rotation :
+# v = np.array([1,0,0])
+# t = np.array([0,1,0])
+#
+# ax.arrow3D(0,0,0,
+#            t[0],t[1],t[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='blue',
+#            ec ='blue')
+#
+# # e_extend
+#
+# rotM = rot(v,np.pi/2)
+#
+# print('rotM:', rotM)
+#
+# rot_t = rotM.dot(t)
+#
+# print('rot_t:', rot_t)
+#
+# ax.arrow3D(0,0,0,
+#            rot_t[0],rot_t[1],rot_t[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='blue',
+#            ec ='blue')
+
+### -------------------------------------------
+
+############################################################################################################################################
+####################################################################### KAPPA POSITIVE  ####################################################
+############################################################################################################################################
+
+kappa = (-1)*kappa
+
+if kappa == 0 :
+    u1 = 0*x1
+    u2 = x1*e[0] + x2*e[1]
+    u3 = x2*e[0] - x1*e[1]
+else :
+    u1 = -(1/kappa)*np.cos(kappa*(x1*e[0]+x2*e[1])) + (1/kappa)
+    u2 = (1/kappa)*np.sin(kappa*(x1*e[0]+x2*e[1]))
+    u3 = x2*e[0] -x1*e[1]
+# ax.plot_surface(u1, u2, u3, color = 'w', rstride = 1, cstride = 1, facecolors=cm.autumn(colorfunction), alpha=.3)  ##This one!
+
+
+# T = rotate_data(X,Rotation_vector,Rotation_angle)
+
+T = rotate_data(np.array([u1,u2,u3]),Rotation)
+# T = rotate_data(T,np.array([0,1,0]),Rotation_angle)
+# T = rotate_data(T,np.array([0,0,1]),-1*Rotation_angle/2)
+
+ax.plot_surface(T[0], T[1], T[2], color = 'w', rstride = 1, cstride = 1, facecolors=cm.autumn(colorfunction), alpha=.4, zorder=4)
+
+
+# midpoint = np.array([(max(x)+min(x))/2,(max(y)+min(y))/2])
+# print('midpoint',midpoint)
+
+# Map midpoint:
+midpoint_mapped = u(midpoint,kappa,e)
+print('mapped midpoint', midpoint_mapped)
+
+# ax.plot(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2], color='black', markersize=10,marker='o',  zorder=5)
+ax.plot(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],    # data
+marker='o',     # each marker will be rendered as a circle
+markersize=4,   # marker size
+markerfacecolor='orange',   # marker facecolor
+markeredgecolor='black',  # marker edgecolor
+markeredgewidth=1,       # marker edge width
+linewidth=1,
+zorder=5)          # line width
+# ax.scatter3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2], color='black', s=100,  zorder=5)
+
+
+# mapped_e = grad_u(midpoint,kappa,e)
+# normal = compute_normal(midpoint,kappa,e)
+
+
+# ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+#            mapped_e[0],mapped_e[1],mapped_e[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='red',
+#            ec ='red')
+#
+# ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+#            normal[0],normal[1],normal[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='blue',
+#            ec ='blue')
+
+############################################################################################################################################
+####################################################################### KAPPA ZERO #########################################################
+############################################################################################################################################
+kappa = 0
+
+if kappa == 0 :
+    u1 = 0*x1
+    u2 = x1*e[0] + x2*e[1]
+    u3 = x2*e[0] - x1*e[1]
+else :
+    u1 = -(1/kappa)*np.cos(kappa*(x1*e[0]+x2*e[1])) + (1/kappa)
+    u2 = (1/kappa)*np.sin(kappa*(x1*e[0]+x2*e[1]))
+    u3 = x2*e[0] -x1*e[1]
+# ax.plot_surface(u1, u2, u3,  rstride = 1, cstride = 1, color = 'white', alpha=0.85)
+
+# T = rotate_data(np.array([u1,u2,u3]),Rotation_vector,Rotation_angle)
+
+T = rotate_data(np.array([u1,u2,u3]),Rotation)
+# T = rotate_data(T,np.array([0,1,0]),Rotation_angle)
+# T = rotate_data(T,np.array([0,0,1]),-1*Rotation_angle/2)
+
+ax.plot_surface(T[0], T[1], T[2],  rstride = 1, cstride = 1, color = 'white', alpha=0.55, zorder=3)
+
+
+# midpoint = np.array([(max(x)+min(x))/2,(max(y)+min(y))/2])
+mapped_e = grad_u(midpoint,kappa,e)
+normal_zeroCurv = compute_normal(midpoint,kappa,e)
+
+# Map midpoint:
+midpoint_mapped = u(midpoint,kappa,e)
+print('mapped midpoint', midpoint_mapped)
+
+##-----  PLOT MAPPED MIDPOINT :::
+
+ax.plot(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],    # data
+marker='o',     # each marker will be rendered as a circle
+markersize=4,   # marker size
+markerfacecolor='orange',   # marker facecolor
+markeredgecolor='black',  # marker edgecolor
+markeredgewidth=1,       # marker edge width
+# linestyle='--',            # line style will be dash line
+linewidth=1,
+zorder=5)
+
+# ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+#            mapped_e[0],mapped_e[1],mapped_e[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='red',
+#            ec ='red')
+#
+# ax.arrow3D(midpoint_mapped[0],midpoint_mapped[1],midpoint_mapped[2],
+#            normal_zeroCurv[0],normal_zeroCurv[1],normal_zeroCurv[2],
+#            mutation_scale=10,
+#            arrowstyle="-|>",
+#            linestyle='dashed',fc='blue',
+#            ec ='blue')
+
+
+##----------  PLOT MAPPED ORIGIN :::
+# origin = np.array([0,0])
+# origin_mapped = u(origin,kappa,e)
+# print('origin_mapped', origin_mapped)
+#
+# ax.plot(origin_mapped[0],origin_mapped[1],origin_mapped[2],    # data
+# marker='o',     # each marker will be rendered as a circle
+# markersize=4,   # marker size
+# markerfacecolor='green',   # marker facecolor
+# markeredgecolor='black',  # marker edgecolor
+# markeredgewidth=1,       # marker edge width
+# linewidth=1,
+# zorder=5)          # line width
+#
+# # rotate mapped origin
+# # v = np.array([1,0,0])
+# # alpha  = Rotation_angle
+#
+# rotM = rot(Rotation_vector,Rotation_angle)
+# # origin_mRot = rotate_data(origin_mapped,v,alpha)
+# origin_mRot = rotM.dot(origin_mapped)
+# print('origin_mapped Rotated', origin_mRot)
+#
+# # --- Compute Distance to Origin 3D
+# origin_3D=np.array([0,0,0])
+# distance = origin_mapped-origin_3D
+# print('distance', distance)
+
+## --------------------------------------------------------
+
+# COMPUTE ANGLE WITH Z AXIS
+z = np.array([0,0,1])
+print('test', normal_zeroCurv*z)
+angle_z = np.arccos(normal_zeroCurv.dot(z) /( (np.linalg.norm(z)*np.linalg.norm(normal_zeroCurv) ) ))
+print('angle between normal and z-axis', angle_z)
+
+## unfinished...
+
+
+
+
+###----------  PLOT :
+plt.axis('off')
+# plt.axis('tight')
+
+# ADD colorbar
+# scamap = plt.cm.ScalarMappable(cmap='inferno')
+# fig.colorbar(scamap)
+
+# ax.colorbar()
+# ax.axis('auto')
+# ax.set_title(r'Cylindrical minimizer_$\kappa$='+ str(kappa)+ '_$e$=' +  str(e))
+# ax.set_title(r'Cylindrical minimizer' + '_$e$=' +  str(e))
+ax.set_xlabel(r"x-axis")
+ax.set_ylabel(r"y-axis")
+ax.set_zlabel(r"z-axis")
+
+# TEST :
+# ax.annotate3D('point 1', (0, 0, 0), xytext=(3, 3), textcoords='offset points')
+# ax.annotate3D('point 2', (0, 1, 0),
+#               xytext=(-30, -30),
+#               textcoords='offset points',
+#               arrowprops=dict(ec='black', fc='white', shrink=2.5))
+# ax.annotate3D('point 3', (0, 0, 1),
+#               xytext=(30, -30),
+#               textcoords='offset points',
+#               bbox=dict(boxstyle="round", fc="lightyellow"),
+#               arrowprops=dict(arrowstyle="-|>", ec='black', fc='white', lw=5))
+
+#######################################################################################################################
+
+u1 = T[0]
+u2 = T[1]
+u3 = T[2]
+
+max_range = np.array([u1.max()-u1.min(), u2.max()-u2.min(), u3.max()-u3.min()]).max() /3
+mid_u1 = (u1.max()+u1.min()) * 0.5
+mid_u2 = (u2.max()+u2.min()) * 0.5
+mid_u3 = (u3.max()+u3.min()) * 0.5
+
+
+ax.set_xlim(mid_u1 - max_range, mid_u1 + max_range)
+ax.set_ylim(mid_u2 - max_range, mid_u2 + max_range)
+ax.set_zlim(mid_u3 - max_range, mid_u3 + max_range)
+
+# ax.set_xlim3d(-2, 2)
+# ax.set_ylim3d(-1.0,3.0)
+# ax.set_zlim3d(-1.5,2.5)
+
+# ax.set_ylim3d(-10,10)
+# ax.set_xlim(mid_u1 - max_range-0.2, mid_u1 + max_range+0.2)
+# ax.set_zlim(mid_u3 - max_range-0.2, mid_u3 + max_range+0.2)
+# ax.set_ylim(mid_u2 - max_range-0.2, mid_u2 + max_range+0.2)
+
+# Figurename = r'Cylindrical minimizer_$\kappa$='+ str(kappa)+ '_$e$=' +  str(e)
+Figurename = r'Cylindrical minimizer' + '_$e$=' +  str(e)
+# plt.savefig("test.png", bbox_inches='tight')
+plt.savefig(Figurename+".png", bbox_inches='tight')
+plt.show()
+
+
+
+# #---------------------------------------------------------------