For magnesium EBSD data, compute Schmid factors and plot #mtexScript

SchmidR2.m

%% Function to calculate the schmid factors for all magnesium slip/twin systems
%code input; EBSD data (ebsd variable), and a vector3d defining the applied
%tensile force (default is in x).

%output a 6 column array containing from left to right the SF of basal
%slip, the SF of prismatic slip, the SF of pyramidal slip, the SF of
%extension twinning, the lowest SF of the four available types and
%the index (column) of the lowest SF of the four systems.  Information for plotting is at
%bottom.

%the code was developed for magnesium, and will require customization of
%all slip systems for any other SF calculations.

%Note that to simplify things, the schmid factor output will always be >0
%for slip but for twinning, schmid factor output will always be -90 to +90, 
%as it is directional.

%note that for magnesium a plane and the normal will not necessarily have
%the same index values (unlike cubic systems), so use the ,'hkl' suffix for
%plane normals and the 'uvw' suffix for directions.

%Rev2-changed the orientations to a variable ori entered at top to simplify
%dataset changes.  Switched to function.

function [TotalArray,BasalArray,PrismaticArray,PyrArray,ExtTwinArray] = SchmidR2(dataset,r)
%Limit to the phase of interest
dataset=dataset('Magnesium');
cs=dataset.CS;
ori=dataset.orientations;

%define tensile axis or for an arbitrary vector use forDir = vector3d(1,-1,0);
if nargin<2
r=xvector;
end

%define  CRSS of the basal, prismatic, pyramidal, extension twin so that the schmid
%factor can be scaled relatively between the systems
%Values below from Figure 11, Chapuis, A. & Driver, J. H. Temperature dependency of slip and twinning in plane strain compressed magnesium single crystals Acta Mater., 2011, 59, 1986-1994 in MPa.   

% set basal slip as base value(1)
BasalFactor=4/4;

%(prism 4.5x harder)
PrismaticFactor=18/4;

%type II pyramidal (10x harder)
PyramidalFactor=40/4;

%{10-12} extension twinning (2.75xharder)
ExtTwinFactor=11/4;
%% Basal slip calculations
%define slip plane and directions
% normal to the slip plane
m1 = Miller(0,0,0,1,cs,'hkl'); 
% slip direction in the slip plane
n1 = Miller(1,1,-2,0,cs,'uvw'); 
n2 = Miller(-2,1,1,0,cs,'uvw');
n3 = Miller(1,-2,1,0,cs,'uvw');

% calculate angle between force and plane normal for each point in the map
BasalPlaneNormal= ori*m1;
theta=cos(angle(r, BasalPlaneNormal,'antipodal'));
%to check the angles between BasalPlaneNormal and r, use max(angle(r, BasalPlaneNormal,'antipodal')/degree)
%also check min, values should range from 0-90 degrees  We use antipodal
%symmetry to enforce this

%calculate the orientation of the 3 slip directions
BasalSlip1=ori*n1;
BasalSlip2=ori*n2;
BasalSlip3=ori*n3;
%calculate the angle between the slip and the force
lmda1=cos(angle(r,BasalSlip1,'antipodal'));
lmda2=cos(angle(r,BasalSlip2,'antipodal'));
lmda3=cos(angle(r,BasalSlip3,'antipodal'));

% Now calculate the schmid factor
tau1=theta.*lmda1;
tau2=theta.*lmda2;
tau3=theta.*lmda3;
%combine the variables and calculate the minimum value and the index where
%it occurs. max(BasalArray) to test.  The first 4 columns should be 0 to
%0.5, and the final should be an integer between 1 and 3
BasalArray=[tau1,tau2,tau3];
[val,I]=max(BasalArray,[],2);
BasalArray=[tau1,tau2,tau3,val,I];

%% Prismatic slip calculations
%define prism slip plane and directions
% normal to the slip plane
m4 = Miller(1,0,-1,0,cs,'hkl'); 
m5 = Miller(0,1,-1,0,cs,'hkl');
m6 = Miller(-1,1,0,0,cs,'hkl');
% slip direction in the slip plane
n4 = Miller(-1,2,-1,0,cs,'uvw'); 
n5 = Miller(2,-1,-1,0,cs,'uvw');
n6 = Miller(-1,-1,2,0,cs,'uvw');
%to check the angle use test=angle(m5,n5,'noSymmetry')/degree

% calculate Schmid factors for each point in the map
%angle between force and normal
PrismaticPlaneNormal4= ori*m4;
PrismaticPlaneNormal5= ori*m5;
PrismaticPlaneNormal6= ori*m6;
theta4=cos(angle(r, PrismaticPlaneNormal4,'antipodal'));
theta5=cos(angle(r, PrismaticPlaneNormal5,'antipodal'));
theta6=cos(angle(r, PrismaticPlaneNormal6,'antipodal'));
%to check the angles between PrismaticPlaneNormal4 and r, use max(angle(r, PrismaticPlaneNormal4,'antipodal')/degree)
%also check min, values should range from 0-90 degrees, 

%calcualte the orientation of the 3 slip directions
PrismaticSlip4=ori*n4;
PrismaticSlip5=ori*n5;
PrismaticSlip6=ori*n6;
%calculate the angle between the slip and the force
lmda4=cos(angle(r,PrismaticSlip4,'antipodal'));
lmda5=cos(angle(r,PrismaticSlip5,'antipodal'));
lmda6=cos(angle(r,PrismaticSlip6,'antipodal'));
%to check the angles between PrismaticPlaneNormal4 and r, use max(angle(r, PrismaticSlip4,'antipodal')/degree)
%also check min, values should range from 0-90 degrees, 

% Now calculate the schmid factor
tau4=theta4.*lmda4;
tau5=theta5.*lmda5;
tau6=theta6.*lmda6;

%combine the variables and calculate the minimum value and the index where
%it occurs. The first 4 columns should be 0 to
%0.5, and the final should be an integer between 1 and 3
PrismaticArray=[tau4,tau5,tau6];
[val,I]=max(PrismaticArray,[],2);
PrismaticArray=[tau4,tau5,tau6,val,I];
%% Pyramidal slip calculations
%define pyramidal slip plane and directions
% normal to the slip plane is 
m7 = Miller(1,1,-2,2,cs,'hkl');  
m8 = Miller(1,-2,1,2,cs,'hkl');
m9 = Miller(-2,1,1,2,cs,'hkl');
m10 = Miller(-1,-1,2,2,cs,'hkl'); 
m11 = Miller(-1,2,-1,2,cs,'hkl');
m12 = Miller(2,-1,-1,2,cs,'hkl');

% slip direction in the slip plane
n7 = Miller(-1,-1,2,3,cs,'uvw');
n8 = Miller(-1,2,-1,3,cs,'uvw');
n9 = Miller(2,-1,-1,3,cs,'uvw');
n10 = Miller(1,1,-2,3,cs,'uvw'); 
n11 = Miller(1,-2,1,3,cs,'uvw');
n12 = Miller(-2,1,1,3,cs,'uvw');

% calculate Schmid factors for each point in the map
%angle between force and normal
PyrPlaneNormal7= ori*m7;
PyrPlaneNormal8= ori*m8;
PyrPlaneNormal9= ori*m9;
PyrPlaneNormal10= ori*m10;
PyrPlaneNormal11= ori*m11;
PyrPlaneNormal12= ori*m12;

theta7=cos(angle(r, PyrPlaneNormal7,'antipodal'));
theta8=cos(angle(r, PyrPlaneNormal8,'antipodal'));
theta9=cos(angle(r, PyrPlaneNormal9,'antipodal'));
theta10=cos(angle(r, PyrPlaneNormal10,'antipodal'));
theta11=cos(angle(r, PyrPlaneNormal11,'antipodal'));
theta12=cos(angle(r, PyrPlaneNormal12,'antipodal'));
%to check the angles between PyrPlaneNormal12 and r, use max(angle(r, PyrPlaneNormal12,'antipodal')/degree)
%also check min, values should range from 0-90 degrees, 

%calculate the orientation of the 3 slip directions
PyrSlip7=ori*n7;
PyrSlip8=ori*n8;
PyrSlip9=ori*n9;
PyrSlip10=ori*n10;
PyrSlip11=ori*n11;
PyrSlip12=ori*n12;

%calculate the angle between the slip and the force
lmda7=cos(angle(r,PyrSlip7,'antipodal'));
lmda8=cos(angle(r,PyrSlip8,'antipodal'));
lmda9=cos(angle(r,PyrSlip9,'antipodal'));
lmda10=cos(angle(r,PyrSlip10,'antipodal'));
lmda11=cos(angle(r,PyrSlip11,'antipodal'));
lmda12=cos(angle(r,PyrSlip12,'antipodal'));

%to check the angles between PyrSlip7 and r, use max(angle(r, PyrSlip7,'antipodal')/degree)
%also check min, values should range from 0-90 degrees, 

% Now calculate the schmid factor
tau7=theta7.*lmda7;
tau8=theta8.*lmda8;
tau9=theta9.*lmda9;
tau10=theta10.*lmda10;
tau11=theta11.*lmda11;
tau12=theta12.*lmda12;

%combine the variables and calculate the minimum value and the index where
%it occurs. The first 6 columns should be 0 to
%0.5, and the final should be an integer between 1 and 6
PyrArray=[tau7,tau8,tau9,tau10,tau11,tau12];
[val,I]=max(PyrArray,[],2);
PyrArray=[tau7,tau8,tau9,tau10,tau11,tau12,val,I];
%% Ext Twin calculations
%define Extension twin plane and directions
% normal to the twin plane
m13 = Miller(-1,0,1,2,cs,'hkl'); 
m14 = Miller(1,0,-1,2,cs,'hkl');
m15 = Miller(1,-1,0,2,cs,'hkl');
m16 = Miller(-1,1,0,2,cs,'hkl'); 
m17 = Miller(0,1,-1,2,cs,'hkl');
m18 = Miller(0,-1,1,2,cs,'hkl');

% slip direction in the slip plane
n13 = Miller(1,0,-1,1,cs,'uvw'); 
n14 = Miller(-1,0,1,1,cs,'uvw');
n15 = Miller(-1,1,0,1,cs,'uvw');
n16 = Miller(1,-1,0,1,cs,'uvw'); 
n17 = Miller(0,-1,1,1,cs,'uvw');
n18 = Miller(0,1,-1,1,cs,'uvw');

% calculate Schmid factors for each point in the map
%angle between force and normal
ExtTwinNormal13= ori*m13;
ExtTwinNormal14= ori*m14;
ExtTwinNormal15= ori*m15;
ExtTwinNormal16= ori*m16;
ExtTwinNormal17= ori*m17;
ExtTwinNormal18= ori*m18;

%note that the 'antipodal' command is not used because twinning is
%directional.  this will give negative values of cos (range 1 to -1) when angle >90 degrees
%calculate the angle between the twin shear plane normal and the force
theta13=cos(angle(r, ExtTwinNormal13));
theta14=cos(angle(r, ExtTwinNormal14));
theta15=cos(angle(r, ExtTwinNormal15));
theta16=cos(angle(r, ExtTwinNormal16));
theta17=cos(angle(r, ExtTwinNormal17));
theta18=cos(angle(r, ExtTwinNormal18));
%to check the angles between ExtTwinNormal16 and r, use max(angle(r, ExtTwinNormal16)/degree)
%also check min, values should range from 0 to +180 degrees,  because antipodal is not used

%calculate the orientation of the 6 slip directions
ExtTwinDir13=ori*n13;
ExtTwinDir14=ori*n14;
ExtTwinDir15=ori*n15;
ExtTwinDir16=ori*n16;
ExtTwinDir17=ori*n17;
ExtTwinDir18=ori*n18;

%note that the 'antipodal' command is not used because twinning is
%directional.  this will give negative values of cos (range 1 to -1) when angle >90 degrees
%calculate the angle between the twin shear direction and the force
lmda13=cos(angle(r,ExtTwinDir13));
lmda14=cos(angle(r,ExtTwinDir14));
lmda15=cos(angle(r,ExtTwinDir15));
lmda16=cos(angle(r,ExtTwinDir16));
lmda17=cos(angle(r,ExtTwinDir17));
lmda18=cos(angle(r,ExtTwinDir18));

%to check the angles between ExtTwinDir16 and r, use max(angle(r, ExtTwinDir16,'antipodal')/degree)
%also check min, values should range from 0-90 degrees, 

% Now calculate the schmid factor
tau13=theta13.*lmda13;
tau14=theta14.*lmda14;
tau15=theta15.*lmda15;
tau16=theta16.*lmda16;
tau17=theta17.*lmda17;
tau18=theta18.*lmda18;

%combine the variables and calculate the max value and the index where
%it occurs.
ExtTwinArray=[tau13,tau14,tau15,tau16,tau17,tau18];
[val,I]=max(ExtTwinArray,[],2);
ExtTwinArray=[tau13,tau14,tau15,tau16,tau17,tau18,val,I];
%% Scale outputs, assemble
%next assemble all the maxima
TotalArray=[(BasalArray(:,4)/BasalFactor),(PrismaticArray(:,4)/PrismaticFactor),(PyrArray(:,7)/PyramidalFactor),(ExtTwinArray(:,7)/ExtTwinFactor)];
%find the maxima and the column in which it exists
[val,I]=max(TotalArray,[],2);
TotalArray=[(BasalArray(:,4)/BasalFactor),(PrismaticArray(:,4)/PrismaticFactor),(PyrArray(:,7)/PyramidalFactor),(ExtTwinArray(:,7)/ExtTwinFactor),val,I];


 
%% To plot the array
%to stop generating these, comment this section out.

%to map the SF scaled by system type
figure;
plot(dataset,TotalArray(:, 5));
setColorRange([0, 0.5]) ;
str=sprintf('lowest SF for any system');
set(gcf,'name',str,'NumberTitle','off');


% to map the active slip or twinning system
%with colorcoding-Basal red, prismatic yellow, pyramidal green, ext twin
%Blue

figure;
colors = zeros(size(TotalArray(:, 6)));   
colors(TotalArray(:, 6) == 1) = 1;     % red (1)
colors(TotalArray(:, 6) == 2) = 40;  % yellow (2)
colors(TotalArray(:, 6) == 3) = 65;      % green (3)
colors(TotalArray(:, 6) == 4) = 100;    %blue(4)
plot(dataset,colors);
setColorRange([0, 100]) ;
cmap = [1 0 0; 1 1 0; 0 1 0; 0 0 1];  
colormap(cmap);
str=sprintf('Lowest SF system; red=basal yellow=prism, green=pyr, etwin=blue');
set(gcf,'name',str,'NumberTitle','off');


%% Additional plots; to get these, uncomment 
% 
% %to see the basal schmid factors
% figure;plot(dataset,BasalArray(:, 4))
% setColorRange([0, 0.5]) ;
% str=sprintf('Basal SF');
% set(gcf,'name',str,'NumberTitle','off');
% 
% %to see the prismatic schmid factors
% figure;plot(dataset,PrismaticArray(:, 4))
% setColorRange([0, 0.5]) ;
% str=sprintf('Prism SF');
% set(gcf,'name',str,'NumberTitle','off');
% 
% %to see the pyramidal schmid factors
% figure;plot(dataset,PyrArray(:, 7))
% setColorRange([0, 0.5]) ;
% str=sprintf('Pyr SF');
% set(gcf,'name',str,'NumberTitle','off');
% 
% % to see the extension twin schmid factors-note negative SF value twins
% %  should not appear in theory, but in practice are quite common due to strain
% % accommodation effects
% 
% figure;plot(dataset,ExtTwinArray(:, 7))
% setColorRange([-0.5, 0.5]) ;
% str=sprintf('E-Twin SF');
% set(gcf,'name',str,'NumberTitle','off');
end