Skip to content
Snippets Groups Projects
Commit 7f4fe101 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Update ClassifyMin.py & rewrite PhaseDiagram Code in Python

parent be6bdf92
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ClassifyMin.py
+ 300
− 126
View file @ 7f4fe101
......@@ -8,192 +8,366 @@ import fileinput
import re
import matlab.engine
import sys
# print(sys.executable)
# from subprocess import Popen, PIPE
def harmonicMean(mu_1,beta,theta):
def harmonicMean(mu_1, beta, theta):
return mu_1*(beta/(theta+(1-theta)*beta))
def arithmeticMean(mu_1,beta,theta):
def arithmeticMean(mu_1, beta, theta):
return mu_1*((1-theta)+theta*beta)
def prestrain_b1(rho_1,beta,alpha):
return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)));
def prestrain_b1(rho_1, beta, alpha, theta):
return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)))
def prestrain_b2(rho_1,beta,alpha):
return (3.0*rho_1/(4.0*((1.0-theta)+ theta*beta)))*(1-theta*(1+beta*alpha));
def prestrain_b2(rho_1, beta, alpha, theta):
return (3.0*rho_1/(4.0*((1.0-theta) + theta*beta)))*(1-theta*(1+beta*alpha))
# define function to be minimized
def f(a1,a2,q1,q2,q3,q12,b1,b2):
# Define function to be minimized
def f(a1, a2, q1, q2, q3, q12, b1, b2):
A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
print(A)
B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
print(B)
a = np.array([a1, a2])
print(a)
tmp = np.dot(A,a)
print(tmp)
tmp2 = np.dot(a,tmp)
print(tmp2)
tmpB = np.dot(B,a)
print(tmpB)
# print(A)
# print(B)
# print(a)
tmp = np.dot(A, a)
# print(tmp)
tmp2 = np.dot(a, tmp)
# print(tmp2)
tmpB = np.dot(B, a)
# print(tmpB)
# print(q1*(b1**2))
# print(q2*(b2**2))
# print(q12*b1*b2)
return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
# Classify Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
def classifyMin(q1, q2, q3, q12, b1, b2, print_output = "True"):
print("Run ClassifyMin...")
# ---- Alternative Version using alpha,beta,theta ,mu_1,rho_1
def classifyMin_geo(alpha,beta,theta,q3,q12,mu_1,rho_1,print_Cases=False, print_Output=False):
q1 = harmonicMean(mu_1, beta, theta)
q2 = arithmeticMean(mu_1, beta, theta)
b1 = prestrain_b1(rho_1, beta, alpha,theta)
b2 = prestrain_b2(rho_1, beta, alpha,theta)
return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output)
# TEST
# def classifyMin_geo(alpha,beta,theta,q3,q12,mu_1, rho_1, print_Cases=False, print_Output=False):
# mu_1 = 1.0
# rho_1 = 1.0
# q12 = 0.0
# q1 = harmonicMean(mu_1, beta, theta)
# q2 = arithmeticMean(mu_1, beta, theta)
# b1 = prestrain_b1(rho_1, beta, alpha,theta)
# b2 = prestrain_b2(rho_1, beta, alpha,theta)
# q3 = q1
#
#
# return classifyMin(q1, q2, q3, q12, b1, b2)
# Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types )
# where
# R1 : unique local (global) minimizer which is not axial
# R2 : continuum of local (global) minimizers which are not axial
# R3 : one or two local (global) minimizers which are axial
# Partition given by
# R1 = E1
# R2 = P1.2
# R3 = E2 U E3 U P1.1 U P2 U H
# -------------------------------------------------------------
def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False): #ClassifyMin_hom?
if print_Output: print("Run ClassifyMin...")
CaseCount = 0
epsilon = sys.float_info.epsilon
epsilon = sys.float_info.epsilon #Machine epsilon
#TEST
# A = np.array([[1, 2], [3, 4]])
# kappa = 1.0
# angle = 1.57
# type = 2
B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
# print('Matrix A:', A)
# print('Matrix B:', B)
# print('Matrix rank of A:', np.linalg.matrix_rank(A))
# print('shape of A:', A.shape)
# print('shape of B:', B.shape)
# print('Matrix [A,B]:', np.c_[A, B])
# print('Matrix rank of [A,B]:', np.linalg.matrix_rank(np.c_[A, B]))
# print('shape of [A,B]:', C.shape)
#
# x = np.linalg.solve(A, B) # works only if A is not singular!!
# print("Solution LGS:", x)
# # print("sym Matrix", sym.Matrix(([A],[B])))
#
# # Test
# C = np.array([[1, 0], [0, 0]])
# d = np.array([5, 0])
# y = np.linalg.lstsq(C, d)[0]
# print("Solution LGS:", y)
# T = np.c_[C, d]
# print('T:', T)
# Trref = sym.Matrix(T).rref()[0]
# Out = np.array(Trref, dtype=float)
# print('rref:', Out)
determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
print("determinant:", determinant)
if print_Cases: print("determinant:", determinant) # TODO ..add everywhere if print_Cases:
# Define values for axial-Solutions (b1*,0) & (0,b2*)
b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
# ------------------------------------ Parabolic Case -----------------------------------
if abs(determinant) < epsilon:
print('parabolic case (determinant equal zero)')
if print_Cases: print('P : parabolic case (determinant equal zero)')
# if print_Cases: print('P : parabolic case (determinant equal zero)')
# check if B is in range of A (TODO)
# check if rank(A) == rank([A,B])
#
# OK this way? (or use Sympy?)
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.c_[A, B]):
if print_Cases: print('P1 (B is in the range of A)')
if (q12+q3)/2.0 <= -1.0*epsilon:
print('should not happen(not compatible with det = 0)')
if (abs(B[0]) < epsilon and abs(B[1]) < epsilon) and (q12+q3)/2.0 >= epsilon:
if print_Cases: print('P1.1')
a1 = 0.0
a2 = 0.0
type = 3
CaseCount += 1
if (abs(B[0]) >= epsilon or abs(B[1]) >= epsilon) and (q12+q3)/2.0 >= epsilon:
# Continuum of minimizers located along a line of negative slope in Lambda
if print_Cases: print('P1.2 (Continuum of minimizers located along a line of negative slope in Lambda) ')
# TODO - what to set for a1, a2 ?
# Just solve Aa* = b (alternatively using SymPy ?)
# we know that A is singular (det A = 0 ) here..
# & we know that there are either infinitely many solutions or a unique solution ...
# ---- determine one via Least Squares
# "If A is square and of full rank, then x (but for round-off error)
# is the “exact” solution of the equation. Else, x minimizes the
# Euclidean 2-norm || b-Ax ||. If there are multiple minimizing solutions,
# the one with the smallest 2-norm is returned. ""
a = np.linalg.lstsq(A, B)[0] # TODO check is this Ok ?
print("Solution LGS: a =", a)
a1 = a[0]
a2 = a[1]
type = 2
CaseCount += 1
else:
if print_Cases: print('P2 (B is not in the range of A)')
# local Minimizers occur on the boundary of Lambda...
# There are at most two local minima and they are either
# (b1_star, 0) or (0, b2_star)
# could also outsource this to another function..
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = b1_star
a2 = 0.0
type = 3 # 1
CaseCount += 1
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = 0
a2 = b2_star
type = 3 # 2
CaseCount += 1
# TODO Problem? angle depends on how you choose?...
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
# Two Minimizers pick one
a1 = b1_star
a2 = 0.0
type = 3 # 4
CaseCount += 1
# ------------------------------------ Elliptic Case -----------------------------------
if determinant >= epsilon:
print('elliptic case (determinant greater zero)')
a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
if print_Cases: print('E : elliptic case (determinant greater zero)')
a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3)) / \
(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2)) / \
(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
prod = a1_star*a2_star
b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
if prod >= epsilon:
print('(E1) - inside Lambda ')
if print_Cases: print('(E1) - inside Lambda ')
a1 = a1_star
a2 = a2_star
type = 1 # non-axial Minimizer
CaseCount += 1
if abs(prod) < epsilon: # same as prod = 0 ? or better use <=epsilon ?
if print_Cases: print('(E2) - on the boundary of Lambda ')
a1 = a1_star
a2 = a2_star
type = 3 # non-axial Minimizer
type = 3 # could check which axis if a1_star or a2_star close to zero.. ?
CaseCount += 1
if abs(prod) < epsilon:
print('(E2) - on the boundary of Lambda ')
# if q2*b2**2 < q1*b1**2: # Needs to be updated to Mixed Term!!! just define f as a function and check value?!
#check
if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
# check
# if f(b1_star,0,q1,q2,q3,q12,b1,b2) < f(0, b2_star, q1,q2,q3,q12,b1,b2):
# a1 = b1_star
# a2 = 0.0
# type = 1
# CaseCount += 1
# if f(b1_star,0,q1,q2,q3,q12,b1,b2) > f(0, b2_star, q1,q2,q3,q12,b1,b2):
# a1 = 0
# a2 = b2_star
# type = 2
# CaseCount += 1
# if f(b1_star,0,q1,q2,q3,q12,b1,b2) = f(0, b2_star, q1,q2,q3,q12,b1,b2):
# # Two Minimizers pick one
# a1 = b1_star
# a2 = 0.0
# type = 4
# CaseCount += 1
if prod <= -1.0*epsilon:
if print_Cases: print('(E3) - Outside Lambda ')
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = b1_star
a2 = 0.0
type = 1
type = 3 # 1
CaseCount += 1
if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = 0
a2 = b2_star
type = 2
type = 3 # 2
CaseCount += 1
if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) = f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
# TODO ...does type4 happen here?
# TODO Problem? angle depends on how you choose?...
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
# Two Minimizers pick one
a1 = b1_star
a2 = 0.0
type = 4
CaseCount += 1
if prod < -1.0*epsilon:
print('(E3) - Outside Lambda ')
if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
a1 = b1_star
a2 = 0.0
type = 1
CaseCount += 1
if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
a1 = 0
a2 = b2_star
type = 2
type = 3 # 4
CaseCount += 1
#TODO ...does type4 happen here?
# ------------------------------------ Hyperbolic Case -----------------------------------
if determinant <= -1.0*epsilon:
if print_Cases: print('H : hyperbolic case (determinant smaller zero)')
# One or two minimizers wich are axial
type = 3
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = b1_star
a2 = 0.0
type = 3 # 1
CaseCount += 1
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
a1 = 0
a2 = b2_star
type = 3 # 2
CaseCount += 1
# TODO can add this case to first or second ..
if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
# Two Minimizers pick one
a1 = b1_star
a2 = 0.0
type = 3 # 4
CaseCount += 1
# ---------------------------------------------------------------------------------------
if (CaseCount > 1):
print('Error: More than one Case happened!')
# compute a3
a3 = math.sqrt(2.0*a1*a2) # ever needed?
return A, angle ,type ,kappa
# compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
e = [math.sqrt((a1/(a1+a2))), math.sqrt((a2/(a1+a2)))]
angle = math.atan2(e[1], e[0])
# compute kappa
kappa = (a1 + a2)
#Test
Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]],dtype=object)
# Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]])
if print_Output:
print('--- Output ClassifyMin ---')
print("Minimizing Matrix G:")
print(Minimizer)
print("angle = ", angle)
print("type: ", type)
print("kappa = ", kappa)
# ---------------------------------------------- Main ---------------------
print('Running Python Code')
return Minimizer, angle, type, kappa
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# --- Input Parameters ----
mu_1 = 1.0
rho_1 = 1.0
alpha = 14.0
beta = 20.0
theta = 1.0/4.0
# ---------------------------------------------- Main ---------------------
# define q1, q2 , mu_gamma, q12
# 1. read from Cell-Output
# 2. define values from analytic formulas (expect for mu_gamma)
q1 = harmonicMean(mu_1,beta,theta)
q2 = arithmeticMean(mu_1,beta,theta)
q12 = 0.0 # (analytical example)
# --- Input Parameters ----
# mu_1 = 1.0
# rho_1 = 1.0
# alpha = 9.0
# beta = 2.0
# theta = 1.0/8.0
# # define q1, q2 , mu_gamma, q12
# # 1. read from Cell-Output
# # 2. define values from analytic formulas (expect for mu_gamma)
# q1 = harmonicMean(mu_1, beta, theta)
# q2 = arithmeticMean(mu_1, beta, theta)
# # TEST
# q12 = 0.0 # (analytical example)
# # q12 = 12.0 # (analytical example)
# set mu_gamma to value or read from Cell-Output
mu_gamma = q1 # TODO read from Cell-Output
b1 = prestrain_b1(rho_1,beta,alpha)
b2 = prestrain_b2(rho_1,beta,alpha)
print("q1: ", q1)
print("q2: ", q2)
print("mu_gamma : ", mu_gamma)
print("b1: ", b1)
print("b2: ", b2)
print("machine epsilon", sys.float_info.epsilon)
## ------- options --------
# print_output = false;
A, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2)
print("Matrix A")
print(A)
print("angle",angle)
print("type",type)
print("kappa", kappa)
Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
print("Test", Test)
# Run several times / make Array :
# Gamma_Values = np.linspace(0.01, 2.5, num=12)
# print(Gamma_Values)
# Make Plot
# plt.figure()
# plt.title(r'$\mu_\gamma(\gamma)$-Plot') # USE MATHEMATICAL EXPRESSIONS IN TITLE
# plt.plot(Gamma_Values, mu_gamma)
# plt.scatter(Gamma_Values, mu_gamma)
# plt.xlabel("$\gamma$")
# plt.ylabel("$\mu_\gamma$")
# plt.legend()
# plt.show()
# mu_gamma = q1 # TODO read from Cell-Output
# b1 = prestrain_b1(rho_1, beta, alpha, theta)
# b2 = prestrain_b2(rho_1, beta, alpha, theta)
# print('---- Input parameters: -----')
# print('mu_1: ', mu_1)
# print('rho_1: ', rho_1)
# print('alpha: ', alpha)
# print('beta: ', beta)
# print('theta: ', theta)
# print("q1: ", q1)
# print("q2: ", q2)
# print("mu_gamma: ", mu_gamma)
# print("q12: ", q12)
# print("b1: ", b1)
# print("b2: ", b2)
# print('----------------------------')
# # print("machine epsilon", sys.float_info.epsilon)
#
#
# # ------- Options --------
# print_Cases = True
# print_Output = True
# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
#
# G, angle, type, kappa = classifyMin_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
#
# Out = classifyMin_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
# print('TEST:')
# Out = classifyMin_geo(alpha, beta, theta)
# print('Out[0]', Out[0])
# print('Out[1]', Out[1])
# print('Out[2]', Out[2])
# print('Out[3]', Out[3])
# #supress certain Outout..
# _,_,T,_ = classifyMin_geo(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
# print('Output only type..:', T)
# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
# print("Test", Test)
# -----------------------------------------------------------------------------------------------------------------------------------------------------------------
src/PhaseDiagram.py 0 → 100644
+ 141
− 0
View file @ 7f4fe101
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 mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
from vtk.util import numpy_support
from pyevtk.hl import gridToVTK
# print(sys.executable)
# from subprocess import Popen, PIPE
# Not Needed:
# def ndm(*args):
# return [x[(None,)*i+(slice(None),)+(None,)*(len(args)-i-1)] for i, x in enumerate(args)]
# ---------------------------------------------- Main ---------------------s
# --- Input Parameters ----
mu_1 = 1.0
rho_1 = 1.0
alpha = 9.0
beta = 2.0
theta = 1.0/8.0
# print('Type of theta', type(theta))
# ------- Options --------
# print_Cases = True
# print_Output = True
make_3D_plot = False
make_2D_plot = False
make_3D_PhaseDiagram = True
make_2D_PhaseDiagram = False
# --- Define q1, q2 , mu_gamma, q12 ---
# 1. read from Cell-Output
# 2. define values from analytic formulas (expect for mu_gamma)
q1 = harmonicMean(mu_1, beta, theta)
q2 = arithmeticMean(mu_1, beta, theta)
# TEST
q12 = 0.0 # (analytical example)
# --- Set mu_gamma to value or read from Cell-Output
mu_gamma = q1 # TODO read from Cell-Output
b1 = prestrain_b1(rho_1, beta, alpha, theta)
b2 = prestrain_b2(rho_1, beta, alpha, theta)
print('---- Input parameters: -----')
print('mu_1: ', mu_1)
print('rho_1: ', rho_1)
print('alpha: ', alpha)
print('beta: ', beta)
print('theta: ', theta)
print("q1: ", q1)
print("q2: ", q2)
print("mu_gamma: ", mu_gamma)
print("q12: ", q12)
print("b1: ", b1)
print("b2: ", b2)
print('----------------------------')
# print("machine epsilon", sys.float_info.epsilon)
# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12, b1, b2, print_Cases, print_Output)
# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
# print("Test", Test)
# ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
if make_3D_PhaseDiagram:
alphas_ = np.linspace(-20, 20, 80)
# betas_ = np.linspace(0.01,40.01,1)
#fix to one value:
betas_ = 2.0;
thetas_ = np.linspace(0.01,0.99,80)
# print('type of alphas', type(alphas_))
# print('Test:', type(np.array([mu_gamma])) )
alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
classifyMin_geoVec = np.vectorize(classifyMin_geo)
G, angles, Types, curvature = classifyMin_geoVec(alphas,betas,thetas,mu_gamma,q12, mu_1, rho_1)
# print('size of G:', G.shape)
# print('G:', G)
# Out = classifyMin_geoVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
VTKOutputName = "./PhaseDiagram3D"
gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
print('Written to VTK-File:', VTKOutputName )
if make_2D_PhaseDiagram:
alphas_ = np.linspace(-20, 20, 80)
# betas_ = np.linspace(0.01,40.01,1)
#fix to one value:
betas_ = 2.0;
beta = np.array([betas_])
thetas_ = np.linspace(0.01,0.99,80)
# print('type of alphas', type(alphas_))
# print('Test:', type(np.array([mu_gamma])) )
alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
classifyMin_geoVec = np.vectorize(classifyMin_geo)
G, angles, Types, curvature = classifyMin_geoVec(alphas,beta,thetas,mu_gamma,q12, mu_1, rho_1)
# print('size of G:', G.shape)
# print('G:', G)
# Out = classifyMin_geoVec(alphas,betas,thetas)
# T = Out[2]
# --- Write to VTK
# VTKOutputName = "./PhaseDiagram2D"
# gridToVTK(VTKOutputName , alphas, thetas, beta, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
# print('Written to VTK-File:', VTKOutputName )
# --- Make 3D Scatter plot
if(make_3D_plot or make_2D_plot):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# pnt3d=ax.scatter(x,y,z,c=z.flat)
colors = cm.plasma(Types)
# pnt3d=ax.scatter(x,y,z,c=T.flat)
if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
# pnt3d=ax.scatter(x,y,z,c=colors.flat)
cbar=plt.colorbar(pnt3d)
cbar.set_label("Values (units)")
ax.set_xlabel('alpha')
ax.set_ylabel('beta')
ax.set_zlabel('theta')
plt.show()
# plt.savefig('common_labels.png', dpi=300)
# print('TEST 3D Scatter')
# print('T:', T)
# print('Type 1 occured here:', np.where(T == 1))
# print('Type 2 occured here:', np.where(T == 2))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment