/* I. Create a FCC lattice II. Determine the projection plane using Miller indices notation III. Calculate the angles between a projection of the normal to the projection plane on 3 different planes (XY, XZ and YZ) and axes X [1;0;0], Y [0;1;0] and Z [0;0;1] respectively IV. Rotate the grid (do the projection) */ #include <iostream>#include <math.h> using std::cout; // Using standart namespaceusing std::cin;using std::endl; class Atom{ public: float x, y, z; // Coordinates of the atoms}; typedef class Atom AtomType; AtomType Atom1 [10][10][10]; // Create a 3 dimensional array to fill the space with atoms of type 1 AtomType Atom2 [10][10][10]; // Create a 3 dimensional array to fill the space with atoms of type 2 AtomType Atom3 [10][10][10]; // Create a 3 dimensional array to fill the space with atoms of type 3 AtomType Atom4 [10][10][10]; // Create a 3 dimensional array to fill the space with atoms of type 4 AtomType Atom11 [10][10][10]; // Create a 3 dimensional array of atoms of type 1 that belong to projection plane AtomType Atom22 [10][10][10]; // Create a 3 dimensional array of atoms of type 2 that belong to projection plane AtomType Atom33 [10][10][10]; // Create a 3 dimensional array of atoms of type 3 that belong to projection plane AtomType Atom44 [10][10][10]; // Create a 3 dimensional array of atoms of type 4 that belong to projection plane //float x, y, z; int i, j, k; // Variables for array Atom [i][j][k] int a, b, c; // Integers, which represent reciprocal values of interceptions of the projection plane // with axes X, Y and Z respectively & also represent a normal to the projection plane void projection(); int main (){ float angleX, angleY, angleZ; // Variables to represent angles between a projection of the normal to the projection plane // on 3 different planes (XY, XZ and YZ) and axes X [1;0;0], Y [0;1;0] and Z [0;0;1] respectively // Do item I // 3 nested 'for' loops create the grid of atoms of type 1, 2, 3 and 4 cout << "Atom1 Atom2 Atom3 Atom4" << endl; cout << "-----------------------------------------------------------------------" << endl; for (i=0; i<10; i++) { for (j=0; j<10; j++) { for (k=0; k<10; k++) { Atom1 [i][j][k].x=i-4; // We get '-4' because we want our origin to be in the centre of array Atom1 [i][j][k].y=j-4; Atom1 [i][j][k].z=k-4; Atom2 [i][j][k].x=i-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along X axis Atom2 [i][j][k].y=j-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along Y axis Atom2 [i][j][k].z=k-4; Atom3 [i][j][k].x=i-4; Atom3 [i][j][k].y=j-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along Y axis Atom3 [i][j][k].z=k-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along Z axis Atom4 [i][j][k].x=i-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along X axis Atom4 [i][j][k].y=j-4; Atom4 [i][j][k].z=k-4+0.5; // We get '+0.5' because the atoms of type 2 are shifted on 0.5 along Z axis cout << Atom1 [i][j][k].x << " " << Atom1 [i][j][k].y << " " << Atom1 [i][j][k].z << " "; cout << Atom2 [i][j][k].x << " " << Atom2 [i][j][k].y << " " << Atom2 [i][j][k].z << " "; cout << Atom3 [i][j][k].x << " " << Atom3 [i][j][k].y << " " << Atom3 [i][j][k].z << " "; cout << Atom4 [i][j][k].x << " " << Atom4 [i][j][k].y << " " << Atom4 [i][j][k].z << " " << endl; } } } // Do item II cout << "Please enter the projection plane using Miller indices notation (e.g. 111),\npressing Enter between integers:"; cin >> a >> b >> c; // Do item III angleX=acos(fabs(a*1+b*0+0*0)/(sqrt(a*a+b*b+0*0)*sqrt(1*1+0*0+0*0))); // The angle is in radians, because any functions from angleY=acos(fabs(0*0+b*1+c*0)/(sqrt(0*0+b*b+c*c)*sqrt(0*0+1*1+0*0))); // 'math.h' that operate on angles use radians as the angleZ=acos(fabs(a*0+0*0+c*1)/(sqrt(a*a+0*0+c*c)*sqrt(0*0+0*0+1*1))); // unit of angle // Do item IV /* To obtain a new position of atoms we need to multiply the current position by rotation matrix. Current position: [AtomN [i][j][k].x, AtomN [i][j][k].y, AtomN [i][j][k].z] // N - any number from 1 to 4 Rotation about X axis matrix: 1 0 0 0 cos(angleX) sin(angleX) 0 -sin(angleX) cos(angleX) Rotation about Y axis matrix: cos(angleY) 0 -sin(angleY) 0 1 0 sin(angleY) 0 cos(angleY) Rotation about Z axis matrix: cos(angleZ) sin(angleZ) 0 -sin(angleZ) cos(angleZ) 0 0 0 1 */ cout << "Atom1 Atom2 Atom3 Atom4" << endl; cout << "-----------------------------------------------------------------------" << endl; // Do rotation about X axis for (i=0; i<10; i++) { for (j=0; j<10; j++) { for (k=0; k<10; k++) { Atom1 [i][j][k].x=Atom1 [i][j][k].x*1+Atom1 [i][j][k].y*0+Atom1 [i][j][k].z*0; Atom1 [i][j][k].y=Atom1 [i][j][k].x*0+Atom1 [i][j][k].y*cos(angleX)+Atom1 [i][j][k].z*(-sin(angleX)); Atom1 [i][j][k].z=Atom1 [i][j][k].x*0+Atom1 [i][j][k].y*sin(angleX)+Atom1 [i][j][k].z*cos(angleX); Atom2 [i][j][k].x=Atom2 [i][j][k].x*1+Atom2 [i][j][k].y*0+Atom2 [i][j][k].z*0; Atom2 [i][j][k].y=Atom2 [i][j][k].x*0+Atom2 [i][j][k].y*cos(angleX)+Atom2 [i][j][k].z*(-sin(angleX)); Atom2 [i][j][k].z=Atom2 [i][j][k].x*0+Atom2 [i][j][k].y*sin(angleX)+Atom2 [i][j][k].z*cos(angleX); Atom3 [i][j][k].x=Atom3 [i][j][k].x*1+Atom3 [i][j][k].y*0+Atom3 [i][j][k].z*0; Atom3 [i][j][k].y=Atom3 [i][j][k].x*0+Atom3 [i][j][k].y*cos(angleX)+Atom3 [i][j][k].z*(-sin(angleX)); Atom3 [i][j][k].z=Atom3 [i][j][k].x*0+Atom3 [i][j][k].y*sin(angleX)+Atom3 [i][j][k].z*cos(angleX); Atom4 [i][j][k].x=Atom4 [i][j][k].x*1+Atom4 [i][j][k].y*0+Atom4 [i][j][k].z*0; Atom4 [i][j][k].y=Atom4 [i][j][k].x*0+Atom4 [i][j][k].y*cos(angleX)+Atom4 [i][j][k].z*(-sin(angleX)); Atom4 [i][j][k].z=Atom4 [i][j][k].x*0+Atom4 [i][j][k].y*sin(angleX)+Atom4 [i][j][k].z*cos(angleX); } } } // Do rotation about Y axis for (i=0; i<10; i++) { for (j=0; j<10; j++) { for (k=0; k<10; k++) { Atom1 [i][j][k].x=Atom1 [i][j][k].x*cos(angleY)+Atom1 [i][j][k].y*0+Atom1 [i][j][k].z*sin(angleY); Atom1 [i][j][k].y=Atom1 [i][j][k].x*0+Atom1 [i][j][k].y*1+Atom1 [i][j][k].z*0; Atom1 [i][j][k].z=Atom1 [i][j][k].x*(-sin(angleY))+Atom1 [i][j][k].y*0+Atom1 [i][j][k].z*cos(angleY); Atom2 [i][j][k].x=Atom2 [i][j][k].x*cos(angleY)+Atom2 [i][j][k].y*0+Atom2 [i][j][k].z*sin(angleY); Atom2 [i][j][k].y=Atom2 [i][j][k].x*0+Atom2 [i][j][k].y*1+Atom2 [i][j][k].z*0; Atom2 [i][j][k].z=Atom2 [i][j][k].x*(-sin(angleY))+Atom2 [i][j][k].y*0+Atom2 [i][j][k].z*cos(angleY); Atom3 [i][j][k].x=Atom3 [i][j][k].x*cos(angleY)+Atom3 [i][j][k].y*0+Atom3 [i][j][k].z*sin(angleY); Atom3 [i][j][k].y=Atom3 [i][j][k].x*0+Atom3 [i][j][k].y*1+Atom3 [i][j][k].z*0; Atom3 [i][j][k].z=Atom3 [i][j][k].x*(-sin(angleY))+Atom3 [i][j][k].y*0+Atom3 [i][j][k].z*cos(angleY); Atom4 [i][j][k].x=Atom4 [i][j][k].x*cos(angleY)+Atom4 [i][j][k].y*0+Atom4 [i][j][k].z*sin(angleY); Atom4 [i][j][k].y=Atom4 [i][j][k].x*0+Atom4 [i][j][k].y*1+Atom4 [i][j][k].z*0; Atom4 [i][j][k].z=Atom4 [i][j][k].x*(-sin(angleY))+Atom4 [i][j][k].y*0+Atom4 [i][j][k].z*cos(angleY); } } } // Do rotation about Z axis for (i=0; i<10; i++) { for (j=0; j<10; j++) { for (k=0; k<10; k++) { Atom1 [i][j][k].x=Atom1 [i][j][k].x*cos(angleZ)+Atom1 [i][j][k].y*(-sin(angleZ))+Atom1 [i][j][k].z*0; Atom1 [i][j][k].y=Atom1 [i][j][k].x*sin(angleZ)+Atom1 [i][j][k].y*cos(angleZ)+Atom1 [i][j][k].z*0; Atom1 [i][j][k].z=Atom1 [i][j][k].x*0+Atom1 [i][j][k].y*0+Atom1 [i][j][k].z*1; Atom2 [i][j][k].x=Atom2 [i][j][k].x*cos(angleZ)+Atom2 [i][j][k].y*(-sin(angleZ))+Atom2 [i][j][k].z*0; Atom2 [i][j][k].y=Atom2 [i][j][k].x*sin(angleZ)+Atom2 [i][j][k].y*cos(angleZ)+Atom2 [i][j][k].z*0; Atom2 [i][j][k].z=Atom2 [i][j][k].x*0+Atom2 [i][j][k].y*0+Atom2 [i][j][k].z*1; Atom3 [i][j][k].x=Atom3 [i][j][k].x*cos(angleZ)+Atom3 [i][j][k].y*(-sin(angleZ))+Atom3 [i][j][k].z*0; Atom3 [i][j][k].y=Atom3 [i][j][k].x*sin(angleZ)+Atom3 [i][j][k].y*cos(angleZ)+Atom3 [i][j][k].z*0; Atom3 [i][j][k].z=Atom3 [i][j][k].x*0+Atom3 [i][j][k].y*0+Atom3 [i][j][k].z*1; Atom4 [i][j][k].x=Atom4 [i][j][k].x*cos(angleZ)+Atom4 [i][j][k].y*(-sin(angleZ))+Atom4 [i][j][k].z*0; Atom4 [i][j][k].y=Atom4 [i][j][k].x*sin(angleZ)+Atom4 [i][j][k].y*cos(angleZ)+Atom4 [i][j][k].z*0; Atom4 [i][j][k].z=Atom4 [i][j][k].x*0+Atom4 [i][j][k].y*0+Atom4 [i][j][k].z*1; cout << Atom1 [i][j][k].x << " " << Atom1 [i][j][k].y << " " << Atom1 [i][j][k].z << " "; cout << Atom2 [i][j][k].x << " " << Atom2 [i][j][k].y << " " << Atom2 [i][j][k].z << " "; cout << Atom3 [i][j][k].x << " " << Atom3 [i][j][k].y << " " << Atom3 [i][j][k].z << " "; cout << Atom4 [i][j][k].x << " " << Atom4 [i][j][k].y << " " << Atom4 [i][j][k].z << " " << endl; } } } projection (); // function projection evaluates whether atoms belong to the projection plane return 0;} void projection (){ for (i=0; i<10; i++) { for (j=0; j<10; j++) { for (k=0; k<10; k++) { // Condition for Atom1 if (Atom1 [i][j][k].x*a + Atom1 [i][j][k].y*b + Atom1 [i][j][k].z*c==1) // Intercept equation of the projection plane is x/m+y/n+z/p=1, but m=1/a; n=1/b; p=1/c; { Atom11 [i][j][k].x=Atom1 [i][j][k].x; Atom11 [i][j][k].y=Atom1 [i][j][k].y; Atom11 [i][j][k].z=Atom1 [i][j][k].z; } else break; // Condition for Atom2 if (Atom2 [i][j][k].x*a + Atom2 [i][j][k].y*b + Atom2 [i][j][k].z*c==1) { Atom22 [i][j][k].x=Atom2 [i][j][k].x; Atom22 [i][j][k].y=Atom2 [i][j][k].y; Atom22 [i][j][k].z=Atom2 [i][j][k].z; } else break; // Condition for Atom3 if (Atom3 [i][j][k].x*a + Atom3 [i][j][k].y*b + Atom3 [i][j][k].z*c==1) { Atom33 [i][j][k].x=Atom3 [i][j][k].x; Atom33 [i][j][k].y=Atom3 [i][j][k].y; Atom33 [i][j][k].z=Atom3 [i][j][k].z; } else break; // Condition for Atom4 if (Atom4 [i][j][k].x*a + Atom4 [i][j][k].y*b + Atom4 [i][j][k].z*c==1) { Atom44 [i][j][k].x=Atom4 [i][j][k].x; Atom44 [i][j][k].y=Atom4 [i][j][k].y; Atom44 [i][j][k].z=Atom4 [i][j][k].z; } else break; } } }}
C++ (Qt)bool OnPlane( const Atom & atom, float a, float b, float c, float d, float epsilon ){ float dist = atom.x * a + atom.y * b + atom.z * c + d; return fabs(dist) < epsilon;}
C++ (Qt)QPointF Atom2Screen( const Atom & atom, int w, int h, // ширина и высота экрана (области вывода) float sizeX, float sizeY ) // размер решетки атомов{ float scale = qMin(w / sizeX, h / sizeY); float x = (atom.x - sizeX / 2) * scale + w / 2; float y = h / 2 - (atom.y - sizeY / 2) * scale; return QPointF(x, y);}
C++ (Qt)// на основании Вашего текстаQVector3D axisX(10, 0, 0);QVector3D axisY(0, 10, 0);QVector3D axisZ(0, 0, 10);QVector3D cntr(-4, -4, -4); void Miller2Plane( int a, int b, int с, QVector <QVector4D> & plane ){ for (int i = 0; i < qMax(a, 1); ++i) { for (int j = 0; j < qMax(b, 1); ++j) { for (int k = 0; k < qMax(c, 1); ++k) { // находим 3 точки плоскости QVector <QVector3D> vec; if (a) vec.push_back(cntr + axisX / (i + 1)); if (b) vec.push_back(cntr + axisY / (j + 1)); if (c) vec.push_back(cntr + axisZ / (k + 1)); assert(vec.size() > 0); if (!a) vec.push_back(vec.back() + axisX); if (!b) vec.push_back(vec.back() + axisY); if (!c) vec.push_back(vec.back() + axisZ); // записываем уравнение плоскости (QVector4D) QVector3D nor = QVector3D::normal(vec[0], vec[1], vec[2]).normalize(); plane.push_back(QVector4D(nor, -QVector3D::dotProduct(nor, vec[0]))); } } }}
plane.push_back(QVector4D(nor, -QVector3D::dotProduct(nor, vec[0])));