Итак, я работаю в кристаллографическом проекте c, и я столкнулся с проблемой. В этом коде я пытаюсь построить обратное пространство моего кристалла. Что еще более важно, некоторые конкретные c плоскости моего кристалла.
Я полагаю, что на этом форуме не так много химиков и физиков, поэтому я попытаюсь объяснить свою проблему без фона. У меня есть код, который строит трехмерный график моих «атомов» (точек) в моем кристалле. Каждый атом имеет индекс H, K и L в трехмерной матрице, и для каждого набора из 3 индексов у вас есть пространственная координата (x, y, z). Пример полного трехмерного графика: 1 Но что мне действительно нужно, так это двухмерный график набора плоскостей в этом пространстве. Пример плоскости с нормалью 1-11: 2 для окончательное представление выглядит так: 3
Итак, вот мой код:
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')
В последней части я попытался посмотреть, могу ли я поиграть с углами обзора 3D-сюжет, чтобы построить мою плоскость и, возможно, после этого исчезнуть ось. Но вы не можете повернуть трехмерный сюжет как хотите. Извините, я знаю, что это недостаточно хорошо объяснено, и код очень длинный, но я буду рад, если кто-то сможет мне помочь. Большое спасибо.
edit1 : вот код для минимального примера. Для входов у вас есть a, b и c, которые являются параметрами решетки моего кристалла (например, a = 5 b = 5 и c = 5 для куба), а также h, k и l, которые являются координатами вектора перпендикулярно плоскости, которую вы хотите (если вы поставите 0 0 0, у вас будет полный кристалл, попробуйте поставить 1 0 0 для простой плоскости и 1 1 0 для одной, я не знаю, как построить в 2d). Попробуйте изменить n = 1 или 2 (в начале кода) для лучшей видимости.
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: Думаю, я найду ответ, поищите алгоритм Грама-Шмидта *. 1021 *