Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
function [hfig, save_density, valid_poles_weight] = gtMakePoleFigure(varargin)
% GTMAKEPOLEFIGURE Make a pole figure given a list of poles
% [hfig, save_density, valid_poles_weight] = gtMakePoleFigure({varargin})
% -----------------------------------------------------------------------
%
% the poles should cover a larger angular area than required, to avoid
% depleted pole density at the edges. Any un-needed poles will be crops to
% improve speed.
%
% can pass in either poles or a density matrix. If the density matrix, YOU
% have to make sure the other variables are consistant with it's size.
% if passing in density also need to pass valid_poles_weight for
% normalising
%
% Loads parameters.mat:cryst(phaseid):spacegroup,latticepar.
%
% INPUT: parse parameter pairs format.
%
% 'phaseid' = phase number <int> {1}
%
% 'poles' = poles from gtReadBoundariesStructure
% [poles, weigths] = gtReadBoundariesStructure(boundaries_structure, vThr, DoSym)
%
% 'range' = {'sst'}/'phi90'/'phi360'
% standard stereographic triangle, or quarter or whole hemisphere.
%
% 'label_poles' = 'all'/'bcc'/'main'/{'none'}
% all hkl poles or just the corners of the sst.
%
% 'mark_poles' = {'all'}/'bcc'/'main'/'none'
%
% 'weights' = can supply a vector the same length as poles to weight the
% pole data {[]}
%
% 'searchR' = smoothing radius in pole figure {0.1 rad}
%
% 'increment' = data point density in pole figure {0.01 rad}
%
% 'plot_figure' = can choose not to plot figure, and use the returned
% density values (non normalised) for something
% {true}
%
% Version 002 26-06-2012 by LNervo
% Add phaseid as input argument for multiple phases materials
%
% Version 001 Sep-2009 by AKing
% could be improved to pass in a list of which hkl poles to plot and to
% label.
% read in the specs
appdefaults.phaseid = 1;
appdefaults.poles = [];
appdefaults.density = [];
appdefaults.valid_poles_weight = [];
appdefaults.weights = [];
appdefaults.range = 'sst'; %or phi90, phi360
appdefaults.label_poles = 'none';
appdefaults.mark_poles = 'all';
appdefaults.increment = 0.01;
appdefaults.searchR = 0.1;
appdefaults.new_figure = true;
appdefaults.plot_figure = true;
app = parse_pv_pairs(appdefaults,varargin);
phaseid = app.phaseid;
poles = app.poles;
density = app.density;
searchR = app.searchR;
inc = app.increment;
weights = app.weights;
valid_poles_weight = app.valid_poles_weight;
% get the parameters file
parameters = [];
load('parameters.mat');
spacegroup = parameters.cryst(phaseid).spacegroup;
if ~isempty(poles)
% if no weighting data passed in
if isempty(weights)
weights=ones(size(poles, 1),1);
end
% before building figure, take only those poles that fall in the range 0-phimax,
% 0-psimax (allowing for seach radius)
if strcmp(app.range, 'sst')
% for standard stereographic triangle
phimax=pi/4;
psimax=atan(sqrt(2));
dum=find(poles(:,2)>-sin(searchR) & poles(:,2)<poles(:,1)+(sqrt(2)*sin(searchR)) & poles(:,3)>poles(:,1)-(sqrt(2)*sin(searchR)));
poles=poles(dum, :);
weights=weights(dum);
valid_poles_weight=sum(weights(find(poles(:,2)>=0 & poles(:,2)<=poles(:,1) & poles(:,3)>=poles(:,1))));
total_solid_angle=pi/12;
elseif strcmp(app.range, 'phi90')
% for a 90 degree section of the pole figure
phimax=pi/2;
psimax=pi/2;
dum=find(poles(:,1)>-sin(searchR) & poles(:,2)>-sin(searchR) & poles(:,3)>cos(psimax)-sin(searchR));
poles=poles(dum, :);
weights=weights(dum);
valid_poles_weight=sum(weights(find(poles(:,1)>=0 & poles(:,2)>=0 & poles(:,3)>=0)));
total_solid_angle=pi/2;
elseif strcmp(app.range, 'phi360')
% for 360 pole figure - reduce only using psi
phimax=2*pi;
psimax=pi/2;
dum=find(poles(:,3)>cos(psimax)-sin(searchR));
poles=poles(dum,:);
weights=weights(dum);
valid_poles_weight=sum(weights(find(poles(:,3)>=0)));
total_solid_angle=2*pi;
else
disp('unrecognised plot range :-(')
return
end
% analyse the pole data
i=1;j=1;
for phi=[0:inc:phimax phimax] %around z-axis
j=1;
for psi=[0:inc:psimax psimax] % out from z-axis
z=cos(psi);
x=sin(psi)*cos(phi);
y=sin(psi)*sin(phi);
test_normal=[x y z];
angle_dev = acos((dot(poles, repmat(test_normal, size(poles,1), 1), 2)));
density(i,j)=sum(weights(find(angle_dev<searchR)));
j=j+1;
end
if mod(i, round(length([0:inc:phimax phimax])/10))==0
fprintf('%d percent done\n', round(100*i/length([0:inc:phimax phimax])))
end
i=i+1;
end
% before normalising, save raw density values
save_density=density;
% normalise density into something meaningful
% total solid angle covered 4pi for a whole sphere
% we have reduced the search area
search_solid_angle=pi*(searchR^2);
% random density is
random_density=valid_poles_weight*(search_solid_angle/total_solid_angle);
density=density./random_density;
else
% we need the range of phi and phi
if strcmp(app.range, 'sst')
% for standard stereographic triangle
phimax=pi/4;
psimax=atan(sqrt(2));
total_solid_angle=pi/12;
elseif strcmp(app.range, 'phi90')
% for a 90 degree section of the pole figure
phimax=pi/2;
psimax=pi/2;
total_solid_angle=pi/2;
elseif strcmp(app.range, 'phi360')
phimax=2*pi;
psimax=pi/2;
total_solid_angle=2*pi;
else
disp('unrecognised plot range :-(')
return
end
% normalise density passed in using valid poles_weight
search_solid_angle=pi*(searchR^2);
random_density=valid_poles_weight*(search_solid_angle/total_solid_angle);
density=density/random_density;
end
if app.plot_figure
% build the figure
% have: density as a function of phi and psi.
% Can convert to x,y,z, and use surf.
% modify for steorographic projection
phi=[0:inc:phimax phimax]; %around z-axis
psi=[0:inc:psimax psimax]; %away from z axis
%radius distorted for stereographic proj
r=tan(psi/2);
[phi,r]=meshgrid(phi,r);
[x,y]=pol2cart(phi,r);
if app.new_figure
hfig=figure();
end
surf(x,y,density')
hold on
shading('interp')
axis equal
if ~strcmp(app.mark_poles, 'none') %draw poles
Yoann Guilhem
committed
% Get symmetry operators
symm = gtCrystGetSymmetryOperators([], spacegroup);
% cubic materials
if spacegroup==225 || spacegroup==229
if strcmp(app.mark_poles, 'main')
uvw_poles=[0 0 1; 1 0 1; 1 1 1];
elseif strcmp(app.mark_poles, 'bcc')
uvw_poles=[0 0 1; 1 0 1; 1 1 1; 1 1 2; 2 1 3];
elseif strcmp(app.mark_poles, 'all')
uvw_poles=[0 0 1; 1 0 1; 1 1 1;...
2 0 1; 2 1 1; 2 2 1;...
3 0 1; 3 1 1; 3 2 0; 3 2 1; 3 2 2; 3 3 1; 3 3 2;...
4 0 1; 4 1 1; 4 2 1; 4 3 0; 4 3 1; 4 3 2; 4 3 3; 4 4 1; 4 4 3];
else
disp('unrecognised "mark_poles" option - should be all/main/none')
end
uvw_poles2=[];
Yoann Guilhem
committed
for i=1:length(symm)
uvw_poles2=[uvw_poles2; uvw_poles*symm(i).g];
end
uvw_poles=[uvw_poles2; -uvw_poles2];
% hexagonal materials
elseif spacegroup==194 || spacegroup==167 || spacegroup==663
uvw_poles_hex=[...
1 0 -1 0;...
0 0 0 1;...
0 1 -1 2;...
0 1 -1 1;...
1 1 -2 3;...
1 1 -2 0];
% permute
uvw_poles_hex2=[];
for i = 1:size(uvw_poles_hex,1)
Yoann Guilhem
committed
uvw_poles_hex2 = [uvw_poles_hex2 ; gtGetReflections(uvw_poles_hex(i,:), spacegroup)];
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
end
uvw_poles_hex=uvw_poles_hex2;
% convert to cartesian
uvw_poles(:,1)= uvw_poles_hex(:,1) + 0.5 * uvw_poles_hex(:,2);
uvw_poles(:,2)= 3^0.5/2 *uvw_poles_hex(:,2);
uvw_poles(:,3)= uvw_poles_hex(:,4);
% account for cell parameters
uvw_poles(:,1)=uvw_poles(:,1)*2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1));
uvw_poles(:,2)=uvw_poles(:,2)*2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1));
uvw_poles(:,3)=uvw_poles(:,3)/parameters.cryst(phaseid).latticepar(3);
% add invert through origin
uvw_poles=[uvw_poles; -uvw_poles];
else
disp('spacegroup not recognised...')
end
% get the relevent uvw_poles
if strcmp(app.range, 'sst')
%for standard stereographic triangle
dum=find(uvw_poles(:,2)>=0 & uvw_poles(:,2)<=uvw_poles(:,1) & uvw_poles(:,3)>=uvw_poles(:,1));
uvw_poles=uvw_poles(dum, :);
elseif strcmp(app.range, 'phi90')
%for a 90 degree section of the pole figure
dum=find(uvw_poles(:,1)>=0 & uvw_poles(:,2)>=0 & uvw_poles(:,3)>=0);
uvw_poles=uvw_poles(dum, :);
elseif strcmp(app.range, 'phi360')
%for 360 pole figure - reduce only using psi
uvw_poles=uvw_poles(find(uvw_poles(:,3)>=0),:);
end
% plot these uvw/hkl poles
phi_uvw=atan2(uvw_poles(:,2), uvw_poles(:,1));
psi_uvw=acos(uvw_poles(:,3)./(sqrt(sum(uvw_poles.*uvw_poles, 2))));
r_uvw=tan(psi_uvw/2);
dummy=find(r_uvw>1);
r_uvw(dummy)=[];
phi_uvw(dummy)=[];
[x_uvw,y_uvw]=pol2cart(phi_uvw, r_uvw);
dummy=find(x_uvw<-1 | x_uvw>1 | y_uvw<-1 | y_uvw>1);
x_uvw(dummy)=[];
y_uvw(dummy)=[];
plot3(x_uvw, y_uvw, ones(size(x_uvw))*(max(density(:))+1), 'k*')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot rings on phi360 figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(app.range, 'phi360')
for ring=15:15:75
ring_phi=(pi/180)*[0:3:360 0];
ring_psi=(pi/180)*ring*ones(size(ring_phi));
% convert to cartesian
ring_r=tan(ring_psi/2);
[ring_x,ring_y]=pol2cart(ring_phi,ring_r);
ring_z=ones(size(ring_phi))*(max(density(:))+1);
plot3(ring_x, ring_y, ring_z, 'w-');
if 0 % label rings
h=text(ring_x(1)+0.01, 0, (max(density(:))+1), num2str(ring));
set(h, 'color', 'w');
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tidy the figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(gca, 'XTickLabel','');
set(gca, 'YTickLabel','');
set(gca, 'GridLineStyle','none');
set(gca, 'Ycolor', 'w');
set(gca, 'Xcolor', 'w');
colorbar();
%change to x-y view
set(gca, 'CameraPosition', [(min(x(:))+max(x(:)))/2 (min(y(:))+max(y(:)))/2 5*(max(density(:))+1)]);
set(gca, 'CameraUpVector', [0 1 0]);
if strcmp(app.range, 'sst')
%for sst - apply a mask to the figure
%line from 101 to 111 is z=x
line_phi=[0:0.01:pi/4 pi/4];
linex=cos(line_phi);
linez=linex;
liney=sin(line_phi);
line_points=[linex' liney' linez'];
%normalise
line_points=line_points./repmat(sqrt(sum(line_points.*line_points, 2)), 1, 3);
%convert to stereographic proj.
line_psi=acos(line_points(:,3));
line_r=tan(line_psi/2);
[line_x,line_y]=pol2cart(line_phi', line_r);
%white patch over unwanted bit of figure
mask_x=[line_x ; sqrt(2)-1+0.045 ; sqrt(2)-1+0.045 ; line_x(1)];
mask_y=[line_y ; line_y(end) ; 0 ; 0];
mask_z=ones(size(mask_y))*(max(density(:))+0.5); %between the pole points and the surface
p=patch(mask_x, mask_y, mask_z, 'w');
set(p, 'edgecolor', 'w')
%black outline on polefigure
mask_x=[line_x ; 0 ; line_x(1)];
mask_y=[line_y ; 0 ; 0];
mask_z=ones(size(mask_y))*(max(density(:))+1);
plot3(mask_x, mask_y, mask_z, 'k')
%tweak axis limits and therefore adjust camera position
set(gca, 'xlim', [-0.02 sqrt(2)-1+0.045])
set(gca, 'ylim', [-0.025 0.3769])
set(gca, 'CameraPosition', [(sqrt(2)-1+0.025)/2 (-0.025+0.3769)/2 (max(density(:))+2)]);
set(gca, 'CameraUpVector', [0 1 0]);
%redo the colorbar, and adjust position
%warning('nice sst colorbar switched off')
c=colorbar;
set(c, 'Position', [0.18 0.3 0.04 0.6])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% label poles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~strcmp(app.label_poles, 'none')
label_z=max(density(:))+1;
if strcmp(app.label_poles, 'all') || strcmp(app.label_poles, 'main') || strcmp(app.label_poles, 'bcc')
%label the poles son the corners of the sst
text(0, 0-0.015, label_z, '(001)')
text(0.366+0.01, 0.366-0.005, label_z, '(111)')
text(0.4142-0.015, -0.015, label_z, '(101)')
end
if strcmp(app.label_poles, 'all')
%label the long list of poles
text(0.1213-0.04, 0.1213, label_z, '(114)')
text(0.1583-0.04, 0.1583, label_z, '(113)')
text(0.2247-0.04, 0.2247, label_z, '(112)')
text(0.2808-0.04, 0.2808, label_z, '(223)')
text(0.3052-0.045, 0.3052, label_z, '(334)')
text(0.3845+0.01, 0.2884, label_z, '(434)')
text(0.3901+0.01, 0.2601, label_z, '(323)')
text(0.4000+0.01, 0.2000, label_z, '(212)')
text(0.4077+0.01, 0.1359, label_z, '(313)')
text(0.4105+0.01, 0.1026, label_z, '(414)')
text(0.1231-0.02, -0.015, label_z, '(104)')
text(0.1623-0.01, -0.015, label_z, '(103)')
text(0.2361-0.015, -0.015, label_z, '(102)')
text(0.3028-0.02, -0.015, label_z, '(203)')
text(0.3333-0.01, -0.015, label_z, '(304)')
text(0.2967+0.01, 0.1483, label_z, '(213)')
text(0.2330+0.01, 0.1165, label_z, '(214)')
text(0.3297+0.01, 0.1099, label_z, '(314)')
text(0.3197+0.01, 0.2131, label_z, '(324)')
elseif strcmp(app.label_poles, 'bcc')
%label bcc slip planes
text(0.2247-0.04, 0.2247, label_z, '(112)')
text(0.3, 0.153, label_z, '(213)')
end
end
else
hfig=[];
end
Yoann Guilhem
committed
end % end of function