%% General inputs
%Note that to simplify things, the schmid factor output will always be >0
%for slip and this code does not distinguish if the slip is in the negative
%direction. 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.
%%define system crystallography
cs=ebsd('Magnesium').CS;
%define tensile axis or for an arbitrary vector use r = vector3d(1,-1,0);
r=xvector;
%define CRSS of the basal, prismatic, pyramidal, extension twin so that the schmid
%factor can be scaled relatively between the systems
%based on 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.
BasalFactor=4/4;
PrismaticFactor=18/4;
%type II pyramidal
PyramidalFactor=40/4;
%{10-12} extension twinning
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= ebsd('Magnesium').orientations*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=ebsd('Magnesium').orientations*n1;
BasalSlip2=ebsd('Magnesium').orientations*n2;
BasalSlip3=ebsd('Magnesium').orientations*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'));
%to check the angles between BasalSlip1 and r, use max(angle(r, BasalSlip1,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,
% 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 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= ebsd('Magnesium').orientations*m4;
PrismaticPlaneNormal5= ebsd('Magnesium').orientations*m5;
PrismaticPlaneNormal6= ebsd('Magnesium').orientations*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=ebsd('Magnesium').orientations*n4;
PrismaticSlip5=ebsd('Magnesium').orientations*n5;
PrismaticSlip6=ebsd('Magnesium').orientations*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 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= ebsd('Magnesium').orientations*m7;
PyrPlaneNormal8= ebsd('Magnesium').orientations*m8;
PyrPlaneNormal9= ebsd('Magnesium').orientations*m9;
PyrPlaneNormal10= ebsd('Magnesium').orientations*m10;
PyrPlaneNormal11= ebsd('Magnesium').orientations*m11;
PyrPlaneNormal12= ebsd('Magnesium').orientations*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=ebsd('Magnesium').orientations*n7;
PyrSlip8=ebsd('Magnesium').orientations*n8;
PyrSlip9=ebsd('Magnesium').orientations*n9;
PyrSlip10=ebsd('Magnesium').orientations*n10;
PyrSlip11=ebsd('Magnesium').orientations*n11;
PyrSlip12=ebsd('Magnesium').orientations*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 slip plane and directions
% normal to the slip 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= ebsd('Magnesium').orientations*m13;
ExtTwinNormal14= ebsd('Magnesium').orientations*m14;
ExtTwinNormal15= ebsd('Magnesium').orientations*m15;
ExtTwinNormal16= ebsd('Magnesium').orientations*m16;
ExtTwinNormal17= ebsd('Magnesium').orientations*m17;
ExtTwinNormal18= ebsd('Magnesium').orientations*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=ebsd('Magnesium').orientations*n13;
ExtTwinDir14=ebsd('Magnesium').orientations*n14;
ExtTwinDir15=ebsd('Magnesium').orientations*n15;
ExtTwinDir16=ebsd('Magnesium').orientations*n16;
ExtTwinDir17=ebsd('Magnesium').orientations*n17;
ExtTwinDir18=ebsd('Magnesium').orientations*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 slip 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 minimum 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
TotalArray=[tau1/BasalFactor,tau2/BasalFactor,tau3/BasalFactor,tau4/PrismaticFactor,tau5/PrismaticFactor,tau6/PrismaticFactor,tau7/PyramidalFactor,tau8/PyramidalFactor,tau9/PyramidalFactor,tau10/PyramidalFactor,tau11/PyramidalFactor,tau12/PyramidalFactor,tau13/ExtTwinFactor,tau14/ExtTwinFactor,tau15/ExtTwinFactor,tau16/ExtTwinFactor,tau17/ExtTwinFactor,tau18/ExtTwinFactor];
[val,I]=max(TotalArray,[],2);
TotalArray=[tau1/BasalFactor,tau2/BasalFactor,tau3/BasalFactor,tau4/PrismaticFactor,tau5/PrismaticFactor,tau6/PrismaticFactor,tau7/PyramidalFactor,tau8/PyramidalFactor,tau9/PyramidalFactor,tau10/PyramidalFactor,tau11/PyramidalFactor,tau12/PyramidalFactor,tau13/ExtTwinFactor,tau14/ExtTwinFactor,tau15/ExtTwinFactor,tau16/ExtTwinFactor,tau17/ExtTwinFactor,tau18/ExtTwinFactor,val,I];
%% to plot these values
%to see the basal schmid factors
%to see the basal schmid factors
figure;plot(ebsd('Magnesium'),BasalArray(:, 4))
%or scaled
figure;plot(ebsd('Magnesium'),BasalArray(:, 4)/BasalFactor)
%to see the prismatic schmid factors
figure;plot(ebsd('Magnesium'),PrismaticArray(:, 4))
%or scaled
figure;plot(ebsd('Magnesium'),PrismaticArray(:, 4)/PrismaticFactor)
%to see the pyramidal schmid factors
figure;plot(ebsd('Magnesium'),PyrArray(:, 7))
%or scaled
figure;plot(ebsd('Magnesium'),PyrArray(:, 7)/PyramidalFactor)
%to see the extension twin schmid factors-note negative values will not produce
%twinning undeer any level of stress
figure;plot(ebsd('Magnesium'),ExtTwinArray(:, 6))
top=caxis;
top=top(1,2);
setColorRange([0, top]) ;
%or scaled
figure;plot(ebsd('Magnesium'),ExtTwinArray(:, 6)/ExtTwinFactor)
top=caxis;
top=top(1,2)
setColorRange([0, top]) ;
%to see the highest schmid factors of all systems, scaled relative to each
%other
figure;plot(ebsd('Magnesium'),TotalArray(:, 19))
%to map the active slip or twinning system
figure;
%with colorcoding-Basal red, prismatic yellow, pyramidal green, ext twin
%Blue
colors = zeros(size(TotalArray(:, 14)));
colors(TotalArray(:, 20) > 0 & TotalArray(:, 20) <= 3) = 1; % red (1)
colors(TotalArray(:, 20) > 3 & TotalArray(:, 20) <= 7) = 40; % yellow (2)
colors(TotalArray(:, 20) > 7 & TotalArray(:, 20) <= 12) = 65; % green (3)
colors(TotalArray(:, 20) > 12 & TotalArray(:, 20) <= 18) = 100; %blue(4)
plot(ebsd('Magnesium'),colors)
cmap = [1 0 0; 1 1 0; 0 1 0; 0 0 1];
colormap(cmap);
%plot an IPF in X
figure;oM = ipdfHSVOrientationMapping(ebsd('Magnesium')); oM.inversePoleFigureDirection =xvector;color = oM.orientation2color(ebsd('Magnesium').orientations); plot(ebsd('Magnesium'),color);
%plot a Pole figure- contour type
h=[Miller(0,0,1,ebsd('Magnesium').CS), Miller(1,0,0,ebsd('Magnesium').CS)];figure;o=ebsd('Magnesium').orientations; cs= ebsd('Magnesium').CS; odf = calcODF(ebsd('Magnesium').orientations,'halfwidth',5*degree); plotPDF(odf,h,'upper','projection','eangle','contourf','minmax');set(gcf,'name','ebsd','NumberTitle','off');
%plot a Pole figure- discrete type
h=[Miller(0,0,1,ebsd('Magnesium').CS), Miller(1,0,0,ebsd('Magnesium').CS)];figure; plotPDF(ebsd('Magnesium').orientations,h,'upper','projection','eangle','points','all', 'MarkerSize',1);set(gcf,'name','ebsd','NumberTitle','off');
save;