%The function calculates the matrix for rotating the elastic stiffness matrix around z-axis.
%Inputs:
    %alpha- rotation angle.
%Outputs:
    %M- the matrix used to rotate elastic stiffness matrix around z-axis.
function M=BondZ(alpha)
M=[cos(alpha)^2 sin(alpha)^2 0 0 0 sin(2*alpha);
   sin(alpha)^2 cos(alpha)^2 0 0 0 -sin(2*alpha);
   0 0 1 0 0 0;
   0 0 0 cos(alpha) -sin(alpha) 0;
   0 0 0 sin(alpha) cos(alpha) 0;
   -sin(2*alpha)/2 sin(2*alpha)/2 0 0 0 cos(2*alpha)];