So i am working in a crystallographic project and I'm in front of a problem. What i'm trying to do in this code is to plot the reciprocal space of my cristal. More importantly some specific planes of my crystal.
I suppose that there are not a lot of chemists and physicists in this forum so I'm going to try to explain my problem without the background. I have a code that build a 3D plot of my "atoms" (points) in my cristal. Each atom has an index H,K and L in a 3D matrix and for each set of 3 indexes you have a space coordinate (x,y,z). Exemple of full 3d plot : 1 But what I really want is a 2D plot of a set of planes in this space exemple of a plane with the normal 1-11 : 2 for a final representation looking like that : 3
So here is my code :
from numpy import *
from math import *
from cmath import *
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.pyplot import *
print('###########################################')
print('Bienvenue sur CR7-Cristallix the 7th wonder')
print('###########################################')
#ensemble des variables
###################################################################################
dir=array([0,0,0])
ar,br,cr=0,0,0 #vecteur réseau réciproque
cif=array([0.,0.,0.])
f=array([1,0,0]) # vecteurs unitaire du réseau direct f,g,h
g=array([0,1,0])
h=array([0,0,1])
n=4 # taille du tenseur
motif=array([0,0,0]) # motif du cristal, prévu pour plus tard si amélioration avec un fichier cif
fd=1. # facteur de diffusion atomique non utilisé pour le moment
v=array([0,0,0])#sert à comparer grace au produit scalaire entre la direction choisie et l'ensemble des vecteurs
y,t=0,0
fig = plt.figure()#utilisés pour définir ax.scatter
ax = fig.add_subplot(111, projection='3d')
###################################################################################
#ensemble des inputs
###################################################################################
for i in range (len(dir)):
dir[i]=int(input('just enter 0 0 1 one by one (it's the normal to the plane:'))
res, a, b, c=input('type of lattice (type F) : '),float(input('a : ')),float(input('b : ')),float(input('c : '))#entre les paramètres de maille et le réseau
n=int(input('dimension of reciprocal space (go for 2or3or4 : '))#entre la taille du tenseur, de hkl max et min. (-n,n)=(h/k/lmin,h/k/lmax)
###################################################################################
#fonction pour créer les vecteurs du réseau réciproque
###################################################################################
def rep():
global a,b,c,ar,br,cr
#vecteurs de maille pour chaque réseau
if res=='P':
a,b,c=a*f,b*g,c*h
elif res =='F':
a,b,c=(b/2)*g+(c/2)*h,(a/2)*f+(c/2)*h,(a/2)*f+(b/2)*g
elif res =='I':
a,b,c=-a/2*f+b/2*g+c/2*h,a/2*f-b/2*g+c/2*h,a/2*f+b/2*g-c/2*h
#calcul du volume de la maille
V=vdot(a,cross(b,c))
#calcul des vecteurs du réseau réciproque
ar,br,cr=2*pi*(1/V)*cross(b,c),2*pi*(1/V)*cross(c,a),2*pi*(1/V)*cross(a,b)
return(a,b,c,ar,br,cr)
rep()
###################################################################################
#fonction qui calcule le facteur de structure et les extinctions
###################################################################################
def facteur(res,cif):#rajouter de nouveaux paramètres à terme
global ext
ext=zeros((2*n+1,2*n+1,2*n+1))
#calcul des raies d'extinction pour chaque réseau avec le facteur de structure
if res=='P':
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
ext[i][j][k]=real(fd*exp(2j*pi*(i*cif[0]+j*cif[1]+k*cif[2])))
elif res =='F':
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
ext[i][j][k]=real(fd*exp(2j*pi*(i*cif[0]+j*cif[1]+k*cif[2]))*(1+exp(1j*pi*(i+j))+exp(1j*pi*(i+k))+exp(1j*pi*(k+j))))/4
elif res =='I':
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
ext[i][j][k]=real(fd*exp(2j*pi*(i*cif[0]+j*cif[1]+k*cif[2]))*(1+exp(1j*pi*(i+j+k))))/2
else:
print('Choisissez entre ces réseaux P, F ou I (mon programmeur est encore trop limité)')
#renvoie un tenseur d'ordre 3 avec pour des indices i,j,k correspondant à h,k,l, une valeur 1 ou 0 selon raie ou extinction
facteur(res,cif)
###################################################################################
#création du réseau réciproque et affichage de celui-ci en fonction du plan indiqué
###################################################################################
coord=zeros((2*n+1,2*n+1,2*n+1,3))
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
coord[i][j][k]=ext[i][j][k]*(ar+br+cr)*array([i,j,k])
print(coord) # bugtest
###################################################################################
###################################################################################
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
v=vdot(coord[i][j][k],dir)
if v==0:#on ne conserve que les vecteurs appartenant au plan de la normal dir
#figure()
s=str((i,j,k))
ax.scatter(coord[i][j][k][0],coord[i][j][k][1],coord[i][j][k][2]) #plot du réseau, manque l'automatisation du
ax.text(coord[i][j][k][0],coord[i][j][k][1],coord[i][j][k][2],s,size=7)
z=array([0,0,1])
e=array([0,1,0])
w=[float(dir[0]),float(dir[1]),float(dir[2])]
calphi=vdot(z,w)/(linalg.norm(z)*linalg.norm(w))
calthe=vdot(z,w)/(linalg.norm(z)*linalg.norm(w))
angphi=real(acos(calphi)*180/pi)
angthe=real(acos(calthe)*180/pi)
# ax.view_init(90,0)
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
In the last part i tried to see if i could play with viewing angles of the 3D plot in order to plot my plane and maybe after that make disappear the axis. But you can't rotate the 3d plot as you want. Sorry I know that it is not well explained and the code is very long but I will be happy if someone can help me. Thank you very much.
edit1: here the code for a minimal exemple. For the inputs you have a,b and c which are the lattice parameters of my cristal (for exemple a=5 b=5 and c=5 for a cube) and h, k and l which are the coordinate of the vector normal to the plane you want (if you put 0 0 0 you'll have the full cristal, try to put 1 0 0 for a simple plane and 1 1 0 for one i don't know how to plot in 2d). Try to change n=1 or 2 (at the beginning of the code) to have a better visibility .
from numpy import *
from math import *
from cmath import *
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.pyplot import *
print('###########################################')
print('Bienvenue sur CR7-Cristallix the 7th wonder')
print('###########################################')
#########################parameters######################################
dir=array([0,0,0])
n=4 #dimension of cristal
v=array([0,0,0])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
coord=zeros((2*n+1,2*n+1,2*n+1,3))
##########################inputs#########################################
a, b, c=float(input('a : ')),float(input('b : ')),float(input('c : '))
for i in range (len(dir)):
dir[i]=int(input('enter the indexis of the normal to the plane you want one by one: '))
###########################create the cells of my cristal################
ext=ones((2*n+1,2*n+1,2*n+1))
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
coord[i][j][k]=ext[i][j][k]*(a+b+c)*array([i,j,k])
############################create the plane perpendicular to the direction i put (dir)####################################
for i in range (-n,n+1):
for j in range (-n,n+1):
for k in range (-n,n+1):
v=vdot(coord[i][j][k],dir)
if v==0:
s=str((i,j,k))
ax.scatter(coord[i][j][k][0],coord[i][j][k][1],coord[i][j][k][2])
ax.text(coord[i][j][k][0],coord[i][j][k][1],coord[i][j][k][2],s,size=7)
###################################################################################
edit2: I think i find the answer, look for the algorithm of Gram-Schmidt.
