How to rotate the unit cell in crystal?
3 次查看(过去 30 天)
显示 更早的评论
Hi,
I have a structure (attached in the question, it is in xyz file), i wnanto to rote it the crystal. Here is my code for rotae
% around the z axis [0 0 1]
x = x1*cosd(ang) - y1*sind(ang);
y = x1*sind(ang) + y1*cosd(ang);
z = z1 ;
What want to make the rotate in 1 1 0 or 0 1 1 or 1 1 1, I know already how to 1 0 0 or 0 1 0
0 个评论
采纳的回答
Alex Hanes
2022-10-24
编辑:Alex Hanes
2022-10-26
You can use quaternions to rotate about an arbitrary axis in 3-dimensions. The quatrotate is capable of performing a rotation of all of your atoms in the unit cell by a quaternion (you will have to specify the quaternion using the x, y, z coordinate of the vector to perform the rotation about and the angle to rotate by). If you want to rotate about an axis K = (kx, ky, kz) by angle theta, then your quaternion for "quatrotate" is:
clear; clc;
% Rotate about the (1,4,3) axis by theta = 60 deg:
theta = pi/3; % (rad), angle to rotate by
K = [1, 4, 3]; % [kx,ky,kz], axis to rotate about
k = K./vecnorm(K); % unit vector
r = [1, 1, 1]; % vector to rotate
% Define Quaternion:
q = [cos(theta/2), sin(theta/2).*k];
vOut = quatRotate(q,r); % call quatRotate()
% Example from quatrotate() page:
q = [1 0 1 0];
r = [1 1 1];
vOut = quatRotate(q,r);
% Use Home-Written quatRotate following quatrotate documentation:
% (https://www.mathworks.com/help/aerotbx/ug/quatrotate.html)
function vOut = quatRotate(q,vIn)
q = q./vecnorm(q); % ensure quaternion is unit vector
q0 = q(1); q1 = q(2); q2 = q(3); q3 = q(4); % sort quaternion components
Q = [ (1 - 2*q2^2 - 2*q3^2), ( 2*q1*q2 + 2*q0*q3), ( 2*q1*q3 - 2*q0*q2); ...
( 2*q1*q2 - 2*q0*q3), (1 - 2*q1^2 - 2*q3^2), ( 2*q2*q3 + 2*q0*q1); ...
( 2*q1*q3 - 2*q0*q2), ( 2*q2*q3 - 2*q0*q1), (1 - 2*q1^2 - 2*q2^2) ];
vOut = Q*vIn';
end
You could also use Euler angles, but the quaternion method is compact and numerically more stable.
EDIT: I wrote a simple script that uses "quatRotate" if you don't have the Aerospace Toolbox. You can modify this to take array inputs for "vIn" as necessary.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!