PF and IPF of Schmid factor (for Mg) #mtexScript

PFandIPFbySchmid.txt

%% Creates pole figure and inverse pole figure maps of schmid factor

% 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.

%note that IPF section is the second half, and variables are not shared.

%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')