% crystal shape demo
% needs @crystalShape avaiable in mtex 4.5.3
% script plots a map and the user can select points where a little crystal shape
% will be dispalyed -crystal will be cc according to the map

mtexdata titanium
[grains,ebsd.grainId]=calcGrains(ebsd,'angle',5*degree); 
% plot something
cs=ebsd('t').CS;
f=newMtexFigure;
om=ipdfHSVOrientationMapping(cs);
om.inversePoleFigureDirection=xvector;
plot(ebsd('t'),om.orientation2color(ebsd('t').orientations))
hold on
plot(grains.boundary)
%% select positions
% click on some positions where you'd like to draw some points
% hit enter/return when done
xy=ginput();
ebsdsel = ebsd(ebsd.findByLocation(xy),'t');

%%
% define crystal shape
m = Miller({1,0,-1,0},cs);  % hexagonal prism
r = Miller({1,0,-1,1},cs);  % positive rhomboedron, usally bigger then z
z = Miller({0,1,-1,1},cs);  % negative rhomboedron
cS = crystalShape([m,r,z],1.2,[1,1.2,1]);

%% plot cs at position 
%elevation
elev=100;
% make sure to plot the crystals above the plane
zdir=getMTEXpref('zAxisDirection');
if strcmp(zdir,'intoPlane'); elev=elev*-1; end

%positon vectors for crystals
posv = vector3d(ebsdsel('t').x,ebsdsel('t').y,elev)
% size of crystal shapes, probably one needs to
% find some good value depending on the maps size
scaling = 100;

% color of crystal shapes at position
colpos = om.orientation2color(ebsdsel.orientations);

% plot at place
mP=getappdata(f.currentAxes,'mapPlot');
for i=1:length(posv)
plot(posv(i) + scaling.*(cS*ebsdsel(i).orientations) ,'FaceColor',colpos(i,:),'parent',mP.ax)
hold on
end
hold off