Schmid Factor-plot on pole figure or ipf

[grandr...@gmail.com]

Hello All, I just put together a code to plot values of schmid factor for various situations and slip systems on pole figures and inverse pole figures.  This was developed for magnesium so it would have to be customized to other situations.  The top half is pole figures, the bottom half is IPF, and no variables are shared between the two.  Feel free to use it and let me know of any issues.



%% Plotting of SF on PF and IPF

% Rev 0-Created Jan 21, 2017 by J. Hiscocks 
%rev 1 fixed IPF point generation issue
%rev 2 added more graph outputs to IPF and titles etc.

%on a pole figure(PF), we can plot the Schmid factor(SF) of a specific
%orientation that has forces applied in a variety of orientations.  For a
%pole figure, it must be ONE orientation and MULTIPLE force directions.  If
%you had multiple orientations, in many cases you could have two orientations
%that result in the same point on the pole figure (ie different Euler angles
%but the same location on a given pole figure).

%note that what you are seeing on the pole figure is actually vector3d
%representing the applied force direction.
 
%% Define some variables
%alternatively if you have no ebsd variable
cs= crystalSymmetry('6/mmm',[3.2,3.2,5.2],'X||a*','Y||b','Z||c');



%% Designate the orientation

%input by Euler angles
ori=orientation('euler',0*degree, 0*degree, 0 *degree,'ZXZ',cs);

%or input by cliking on an existing EBSD map
% disp('select the orientation to use')
%  [x,y]=ginput(1);
%  ori=ebsd(x,y).orientations;

%or input by taking the mode of the pole figure for the whole map (ie the
%orientation of greatest intensity)
%first calculate the odf to generate continuous data.
% odf=calcODF(ebsd('Magnesium').orientations,'halfwidth',5*degree);
% %select the maxima from the odf, and load the 3d vector into the
% % variable 'modes', and skip return of 'values' by using ~ in that place
% [modes, ~] = calcModes(odf,1);
% %transforms the mode output into the appropriate crystallographic vector
% ori = modes(1);



%%  Calculate the force directions

%pick a value every 5 degrees, and label this variable x
% define a grid of tension directions
r = plotS2Grid('resolution',5*degree);
%reshape this into a linear array
r=r(:);
% %this can then be visualised on the pole figure
% h=[Miller(0,0,1,cs), Miller(1,0,0,cs)]
% figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
% hold on
% %use plot instead of plotPDF because we are plotting vector3D here instead
% %of orientations
% plot(r)
% hold off


%% Calculate the schmid factors for all possible slip/twinning systems

%define Basal 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 the orientation in
% question
BasalPlaneNormal= ori*m1;
%calculate the angles between each force (r) and the normal
theta=cos(angle(r, BasalPlaneNormal,'antipodal'));

%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];

%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];

%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];

%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];

%% Plot either the SF for one type of slip  OR plot the lowest value


%plot for basal slip the schmid factors
%this can then be visualised on the pole figure
colour=BasalArray(:,4);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
set(gcf,'name','SF values for Basal Slip','NumberTitle','off');
setColorRange([0, 0.5])


%plot for prism slip the schmid factors
%this can then be visualised on the pole figure
colour=PrismaticArray(:,4);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
set(gcf,'name','SF values for Prismatic Slip','NumberTitle','off');
setColorRange([0, 0.5]);

%plot for pyramidal slip the schmid factors
%this can then be visualised on the pole figure
colour=PyrArray(:,7);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
set(gcf,'name','SF values for Pyramidal Slip','NumberTitle','off');
setColorRange([0, 0.5]);

%plot for Extension twinning the schmid factors
%this can then be visualised on the pole figure
colour=ExtTwinArray(:,7);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
set(gcf,'name','SF values for Extension Twinning','NumberTitle','off');
setColorRange([0, 0.5])

%plot for minimum CRSS
%first define the relative CRSS of the different types of slip, scaled such
%that basal factor =1
BasalFactor=4/4;
PrismaticFactor=18/4;
%type II pyramidal
PyramidalFactor=40/4;
%{10-12} extension twinning
ExtTwinFactor=11/4;

%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];

%plot the max schmid factors for any system scaled relative to each other
%this can then be visualised on the pole figure
colour=TotalArray(:,5);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
set(gcf,'name','max SF values for all systems','NumberTitle','off');
setColorRange([0, 0.5]);

%plot the type of system that is active
%plot the max schmid factors for any system scaled relative to each other
%this can then be visualised on the pole figure
colour=TotalArray(:,6);
h=[Miller(0,0,1,cs), Miller(1,0,0,cs)];
figure;plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',1);
hold on
%use plot instead of plotPDF because we are plotting vector3D here instead
%of orientations
plot(r,colour)
plotPDF(ori,h,'upper','projection','eangle','points','all', 'MarkerSize',10)
hold off
%1 is basal, 2 is prism, 3 is pyramidal, 4 is extension twinning
set(gcf,'name','Active system: dk blue=basal, lt blue=prism, green=pyramidal, yellow=Etwin','NumberTitle','off');

%% for IPF
% In this case we're examining the relationship between a single applied force
% and all the orientations within a crystal

%define the symmetry of magnesium
cs= crystalSymmetry('6/mmm',[3.2,3.2,5.2],'X||a*','Y||b','Z||c');

%create a grid of orientations
ori=equispacedSO3Grid(cs,'resolution',5*degree);

%make this variable linear for ease of tracking max values
ori=ori(:);

%define the force direction
r=xvector;

%% calculate the schmid factors

%define Basal 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 the orientation in
% question
BasalPlaneNormal= ori*m1;
%calculate the angles between each force (r) and the normal
theta=cos(angle(r, BasalPlaneNormal,'antipodal'));

%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];

%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];

%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];

%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];

%% plot the results


%plot for basal slip the schmid factors
%this can then be visualised on the pole figure
colour=BasalArray(:,4);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for Basal Slip','NumberTitle','off');
setColorRange([0, 0.5]);
title('SF values for Basal Slip')

%plot for prism slip the schmid factors
%this can then be visualised on the pole figure
colour=PrismaticArray(:,4);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for Prism Slip','NumberTitle','off');
setColorRange([0, 0.5]);
title('SF values for Prismatic Slip')

%plot for pyramidal slip the schmid factors
%this can then be visualised on the pole figure
colour=PyrArray(:,7);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for Pyramidal Slip','NumberTitle','off');
setColorRange([0, 0.5]);
title('SF values for Pyramidal Slip')

%plot for Extension twinning the schmid factors
%this can then be visualised on the pole figure
colour=ExtTwinArray(:,7);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for Extension Twinning','NumberTitle','off');
setColorRange([0, 0.5]);
title('SF values for Extension Twinning')

%plot for minimum CRSS
%first define the relative CRSS of the different types of slip, scaled such
%that basal factor =1
BasalFactor=4/4;
PrismaticFactor=18/4;
%type II pyramidal
PyramidalFactor=40/4;
%{10-12} extension twinning
ExtTwinFactor=11/4;

%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];

%plot the max schmid factors for any system scaled relative to each other
%this can then be visualised on the pole figure
colour=TotalArray(:,5);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for All systems, scaled','NumberTitle','off');
setColorRange([0, 0.5]);
title('Scaled SF for all systems')

%plot the type of system that is active
%plot the max schmid factors for any system scaled relative to each other
%this can then be visualised on the pole figure
colour=TotalArray(:,6);
figure;plotIPDF(ori,r,'upper','points','all', 'MarkerSize',10,'property',colour);
set(gcf,'name','SF values for All systems, scaled','NumberTitle','off');
%1 is basal, 2 is prism, 3 is pyramidal, 4 is extension twinning
set(gcf,'name','Active system: dk blue=basal, lt blue=prism, green=pyramidal, yellow=Etwin','NumberTitle','off');
title('Active system')

[Thomas Simm]

Hi Jessica,

Have you tried the (new) slip systems and Schmid factor functionality in MTEX? I think this can make your code easier? The directonality of twins is also incoporated in the code.

https://mtex-toolbox.github.io/files/doc/PlasticDeformation.html

Thomas

[Thomas Simm]

Actually not the best link as it doesn't include slip systems, but should be in Ralf's presentation from MTEX2017

[grandr...@gmail.com]

I did a comparison of two codes, mine to one written by Ralf using this functionality a while back (see the thread labeled Schmid Factor Code)
There were some issues with the slip systems generated at the time, so I wanted to stick to the explicit version of the code because I was certain exactly what calculations were happening.  There is no doubt that Ralf's version will be much much shorter input code, but this also let me link my code (SchmidR2 on gist) to the twinning code and the schmid pole figure code I have there, so that the variants were matched up (eg twin variant #1 from each code referred to the same axis). 

Jessica

[Thomas Simm]

Aah I see I missed that post,

As Ralf said you can create your own slip system like this:

if you type >> edit slipSystem.m

and then below 'methods' add a new function e.g. 

 

function sS = myslipsystem(cs,varargin)

      % <11-20>{0001}

      crss = ones(12,1);

      crss(1:6) = 1e99;

      sS = slipSystem(Miller(1,1,-2,0,cs,'UVTW'),Miller(0,0,0,1,cs,'HKIL'),crss);

      sS = symmetrise(sS);

      sS.CRSS = [1 2 3 1 2 3 1 2 3 1 2 3];

      

    end

you can then call it as slipSystem.myslipsystem(CS) 

the line  sS.CRSS = allows you to take account of differing CRSSs of a slip system e.g. for a twin or if you want to stack different slip system types

[Thomas Simm]

Sorry should be:

function sS = myslipsystem(cs,varargin)

      % <11-20>{0001}

      sS = slipSystem(Miller(1,1,-2,0,cs,'UVTW'),Miller(0,0,0,1,cs,'HKIL'));

      sS = symmetrise(sS);

      sS.CRSS = [1 2 3 1 2 3 1 2 3 1 2 3];

      

    end

[Thomas Simm]

Very nice scripts, 

If you added your calculated properties to grains (e.g. grains.prop.daughterArea) this can help keep things together and makes some analysis easier to follow. You could also reconstruct the grain using the parent (and its GOS and orientation), have a parentID for each twin and a list of daughters for each parent. This then makes it easy to plot things like orientation of parents against twin number or area.

We discussed different differentiation methods of twins in an earlier post (i.e. it is not obvious what is the daughter and the parent when using merge), so would be interesting to see how successful your differentiation is (I think you use length of boundaries but also have a checking facility (incorporating ginput)). 

Taylor will probably not be that helpful. Also another idea is to use a stress matrix (instead of just [1 0 0; 0 0 0; 0 0 0] try [1 0 0; 0 -.5 0; 0 0 -.5]) that includes compression in the y and z directions (if pulling in the x-direction) 

Thomas

"Hello Thomas, I updated the document to say the scripts have been uploaded to gist with the tag #mtexScript instead of in an appendix.  I did this so that I can make frequent updates easily.  The schmid factor function I uploaded there (SchmidR2) is capable of calculating the SF for all 6 variants (set up for magnesium) and returning that as well as the lowest SF on any system.  I have recently done a fair bit of work on SF and twin variants, and agree there is some, but not much influence of SF on twin variants present.

My larger scripts on that site (Parent_Twin_AreaR8.m and Parent_Twin_InteractR18.m) are designed to tabulate the statistics for twin variant formation, grain size etc so that you can graph any grain information against SF of the parent, of the variants formed, of the SF of the variants formed...etc. Any factors that could influence twinning which I could think of, but I got my statistics that way.  No Taylor factor information though.

Regards, Jessica"

[grandr...@gmail.com]

If you added your calculated properties to grains (e.g. grains.prop.daughterArea) this can help keep things together and makes some analysis easier to follow.

The output of the script is two arrays, one which could be combined with mergedGrains (and has one line per mergedGrain) and one which could be combined with grains (and has one line per grain).  At the time I didn't consider combining them, but there are some advantages to keeping them separate.  By outputting them as arrays, I can import them into excel which allows me to easily sort and filter them (eg look only at grains where one twin variant is present vs grains where 2-6 are).  Of course you can do this in Matlab as well, but I am more familiar with Excel AND this allows me to combine results from multiple datasets into one larger dataset.  Since the ID for each mergedGrain or grain is stored in the output array (along w/ parentID), it would be easy to copy the information over (and duplicate it) to allow the user to deal with it in whatever fashion they prefer.

I think you use length of boundaries but also have a checking facility (incorporating ginput)

Not quite; for the script with no user input, the largest grain of a given mergedGrain is automatically selected as the parent.  The script with user input starts with that as the default, and in cases where the mergedGrain is large enough to be significant and the default parent not 95% of the grain area, the user is asked for verification.  I find that I manually altered about 10% of these verification cases as a rough estimate.

The boundary checking at the end of the scripts is actually for identification of grains that have undergone secondary twinning, which is done as a last stage (after twins are identified).

On the whole I found the script satisfactory, but there are occasional issues stemming from the wide tolerances that must be used due to grain misorientations across large grains, or boundaries where some segments meet the twin criteria but not most of them.  I suspect it could use some fine tuning, but I found it to work well in almost all of the cases, which can be checked by plotting a grain and having it label them as a twin or parent or secondary twin.

I also found it handy to overlay the location of all twins over the map so I can check the spatial distribution of them.