% A Matlab script to calculate silicon Young'g modulus and Poisson's ration in
% any crystal direction. This script was originally downloaded from
% http://www.kaajakari.net/PracticalMEMS
%
% You may modify this script to your own use and you may distribute this script 
% unmodified. 
% 
% Author: Ville Kaajakari
% Date: 2.4.2009

function yangle3


%110
Cr=Crotated(0,0);
% By definitions
S=inv(Cr);
Y100=1/S(1,1)

%110
Cr=Crotated(pi/4,0);
% By definitions
S=inv(Cr);
Y110=1/S(1,1)

%111
theta=-pi/4;
gamma=-asin(1/sqrt(3));
Cr=Crotated(theta,gamma);
% By definitions
S=inv(Cr);
Y111=1/S(1,1)


N=101;
theta=linspace(0,pi/2,N);
v_angle1=zeros(1,N);
v_angle2=zeros(1,N);
for kk=1:N,
    
    
    Cr=Crotated(theta(kk),0);
    % By definitions
    S=inv(Cr);
    Y_angle(kk)=1/S(1,1);
    v_angle1(kk)=-S(1,2)/S(1,1); % Between two (100)-plane vectors
    v_angle2(kk)=-S(1,3)/S(1,1); % Between (100)-plane vector and [001]-vector
    G_angle1(kk)=1/S(5,5); % Between (100)-plane vector and [001]-vector
    G_angle2(kk)=1/S(6,6); % Between two (100)-plane vectors

    
end

figure(1);
h=polar(theta,Y_angle);
ha=axis;
axis([0 ha(2) 0 ha(4)]);

figure(2);
polar(theta,v_angle2)
hold on;
polar(theta,v_angle1);
hold off;
ha=axis;
axis([0 ha(2) 0 ha(4)]);

figure(3);
polar(theta,G_angle1)
hold on;
polar(theta,G_angle2);
hold off;
ha=axis;
axis([0 ha(2) 0 ha(4)]);


function Crot=Crotated(theta,gamma)
% theta is rotation aroung z-axis
% gamma is rotation around y-axis

c1111=1.66;
c1122=0.64;
c2323=0.80;

%Non-zero componenents of Si stiffness tensor in [100] coordinates
Cc(1,1,1,1)=c1111;
Cc(2,2,2,2)=c1111;
Cc(3,3,3,3)=c1111;

Cc(1,1,2,2)=c1122;
Cc(1,1,3,3)=c1122;
Cc(2,2,1,1)=c1122;
Cc(2,2,3,3)=c1122;
Cc(3,3,1,1)=c1122;
Cc(3,3,2,2)=c1122;

Cc(1,3,1,3)=c2323;
Cc(3,1,1,3)=c2323;
Cc(1,3,3,1)=c2323;
Cc(3,1,3,1)=c2323;

Cc(2,3,2,3)=c2323;
Cc(3,2,2,3)=c2323;
Cc(2,3,3,2)=c2323;
Cc(3,2,3,2)=c2323;

Cc(1,2,1,2)=c2323;
Cc(2,1,1,2)=c2323;
Cc(1,2,2,1)=c2323;
Cc(2,1,2,1)=c2323;
% Rotation matrix
Q1=[
    [cos(theta) sin(theta) 0];
    [-sin(theta) cos(theta) 0];
    [0 0 1];
    ];

l1=cos(gamma); m1=0; n1=sin(gamma);
l2=0;             m2=1;             n2=0;
l3=-sin(gamma); m3=0; n3=cos(gamma);

Q2=[
    [l1 m1 n1];
    [l2 m2 n2];
    [l3 m3 n3];
];
Q=Q1*Q2;
%Carculate stiffness tensor in rotated coordinates
Crot=zeros(3,3,3,3);
for i=1:3,
    for j=1:3,
        for k=1:3,
            for l=1:3,
                for p=1:3,
                    for q=1:3,
                        for r=1:3,
                            for s=1:3,
                                Crot(i,j,k,l)=Crot(i,j,k,l)+Q(p,i)*Q(q,j)*Q(r,k)*Q(s,l)*Cc(p,q,r,s);
                            end
                        end
                    end
                end
            end
        end
    end
end

Crot=[
    [Crot(1,1,1,1) Crot(1,1,2,2) Crot(1,1,3,3) Crot(1,1,2,3) Crot(1,1,1,3) Crot(1,1,1,2)];
    [Crot(2,2,1,1) Crot(2,2,2,2) Crot(2,2,3,3) Crot(2,2,2,3) Crot(2,2,1,3) Crot(2,2,1,2)];
    [Crot(3,3,1,1) Crot(3,3,2,2) Crot(3,3,3,3) Crot(3,3,2,3) Crot(3,3,1,3) Crot(3,3,1,2)];
    [Crot(3,2,1,1) Crot(3,2,2,2) Crot(3,2,3,3) Crot(2,3,2,3) Crot(2,3,1,2) Crot(2,3,1,2)];
    [Crot(3,1,1,1) Crot(3,1,2,2) Crot(3,1,3,3) Crot(3,1,2,3) Crot(1,3,1,3) Crot(3,1,1,2)];
    [Crot(2,1,1,1) Crot(2,1,2,2) Crot(2,1,3,3) Crot(2,1,2,3) Crot(2,1,1,3) Crot(1,2,1,2)];
    ];