Skip to content

Commit b34e652

Browse files
beaudubeaudu
committed
improves pCDM source plotting
1 parent 1b6d563 commit b34e652

File tree

2 files changed

+171
-55
lines changed

2 files changed

+171
-55
lines changed

‎nikkhoo/ellipsoid3.m

+71
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
function [xx,yy,zz]=ellipsoid3(varargin)
2+
%ELLIPSOID3 Generate 3D rotated ellipsoid.
3+
% [X,Y,Z]=ELLIPSOID3(AX,AY,AZ,RX,RY,RZ,N) generates three
4+
% (N+1)-by-(N+1) matrices so that SURF(X,Y,Z) produces an
5+
% ellipsoid centered at (0,0,0) with semi-axis radii AX, AY, AZ,
6+
% and rotation angles RX, RY, RZ around each axis (in degree).
7+
%
8+
% [X,Y,Z]=ELLIPSOID3(AX,AY,AZ,RX,RY,RZ) uses N = 20.
9+
%
10+
% ELLIPSOID3(...) and ELLIPSOID3(...,N) with no output arguments
11+
% graph the ellipsoid as a SURFACE and do not return anything.
12+
%
13+
% ELLIPSOID3(AX,...) plots into AX instead of GCA.
14+
%
15+
% The ellipsoidal data is generated using the equation:
16+
%
17+
% X^2 Y^2 Z^2
18+
% -------- + -------- + -------- = 1
19+
% AX^2 AY^2 AY^2
20+
%
21+
% thus applies a 3D rotation using RX, RY and RZ angles.
22+
%
23+
% See also SPHERE, ELLIPSOID, CYLINDER.
24+
%
25+
%
26+
% Author: François Beauducel <[email protected]>
27+
% Created: 2020-06-07 in Yogyakarta, Indonesia
28+
%
29+
% Based on the original script by Laurens Schalekamp and Damian T. Packer
30+
% Copyright 1984-2002 The MathWorks, Inc.
31+
32+
% Parse possible Axes input
33+
narginchk(6,8);
34+
[cax,args,nargs] = axescheck(varargin{:});
35+
36+
[ax,ay,az,rx,ry,rz] = deal(args{1:6});
37+
n = 20;
38+
39+
if nargs > 6
40+
n = args{7};
41+
end
42+
43+
% generates a unitary sphere
44+
[x,y,z] = sphere(n);
45+
46+
% extends semi-axis radii
47+
x = ax*x;
48+
y = ay*y;
49+
z = az*z;
50+
51+
% rotation matrix
52+
Rx = [1 0 0;0 cosd(rx) sind(rx);0 -sind(rx) cosd(rx)];
53+
Ry = [cosd(ry) 0 -sind(ry);0 1 0;sind(ry) 0 cosd(ry)];
54+
Rz = [cosd(rz) sind(rz) 0;-sind(rz) cosd(rz) 0;0 0 1];
55+
R = Rz*Ry*Rx;
56+
57+
xyz = R*[x(:),y(:),z(:)]';
58+
59+
sz = size(x);
60+
x = reshape(xyz(1,:),sz);
61+
y = reshape(xyz(2,:),sz);
62+
z = reshape(xyz(3,:),sz);
63+
64+
if(nargout == 0)
65+
cax = newplot(cax);
66+
surf(x,y,z,'parent',cax)
67+
else
68+
xx = x;
69+
yy = y;
70+
zz = z;
71+
end

‎nikkhoo/plotpcdm.m

+100-55
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function varargout = plotpcdm(p,w,pview,varargin)
1+
function varargout = plotpcdm(p,w,varargin)
22
%PLOTPCDM Plot 3-D pCDM source
33
% PLOTPCDM([X0,Y0,Z0,OX,OY,OZ,A,B],W) adds to current axe a 3-D
44
% representation of pCDM source, with following parameters:
@@ -8,22 +8,57 @@
88
% B = vertical volume variation ratio
99
% W = source maximum width
1010
%
11+
% PLOTPCDM(...,'ellipsoid') plots a 3D ellipsoid instead of 3 rectangular
12+
% plans. Needs ELLIPSOID3 function.
13+
%
1114
% PLOTPCDM(...,'3d') plots the source in 3-D (default).
1215
% PLOTPCDM(...,'xy') plots a projection of the source in XY plan.
1316
% PLOTPCDM(...,'xz') plots a projection of the source in XZ plan.
1417
% PLOTPCDM(...,'zy') plots a projection of the source in ZY plan.
1518
%
1619
%
1720
% Authors: François Beauducel and Antoine Villié
21+
% inspired by the plotCDM function by Mehdi Nikkhoo (2016)
22+
1823
% Created: 2018-07-13, in Yogyakarta (Indonesia)
19-
% Updated: 2019-07-31
24+
% Updated: 2020-06-08
2025

26+
pview = '3d';
27+
pplan = true;
28+
hflag = ishold;
29+
if ~hflag
30+
hold on
31+
end
2132

33+
if nargin > 2
34+
kk = cellfun(@ischar,varargin);
35+
36+
% view type input argument
37+
k = ismember(varargin(kk),{'3d','xy','xz','zy'});
38+
if any(k)
39+
pview = varargin{kk(k)};
40+
varargin(kk(k)) = [];
41+
end
42+
43+
% source type input argument
44+
k = strcmp(varargin(kk),'ellipsoid');
45+
if any(k)
46+
pplan = false;
47+
varargin(kk(k)) = [];
48+
end
49+
end
2250

23-
if nargin < 3
24-
pview = '3d';
51+
if isscalar(w)
52+
w = repmat(w,size(p,1),1);
2553
end
26-
opt = {'LineWidth',2,'EdgeColor','black','FaceColor','none'};
54+
55+
%opt = {'LineWidth',2,'EdgeColor','black','FaceColor','none'};
56+
opt = {'LineWidth',.1,'EdgeColor','black'};
57+
58+
% default dislocation colors by Nikkhoo
59+
dcolx = [255 192 0]/255;
60+
dcoly = [0 176 80]/255;
61+
dcolz = [0 112 192]/255;
2762

2863
for n = 1:size(p,1)
2964
x0 = p(n,1);
@@ -35,64 +70,74 @@
3570
a = p(n,7);
3671
b = p(n,8);
3772

38-
% converts A and B aspect ratio to DVX/DVY/DVZ
39-
dvz = a;
40-
dvy = (1-dvz)*b;
41-
dvx = (1-dvz)*(1-b);
73+
% converts A and B aspect ratio to semi-axis DVX/DVY/DVZ with unitary
74+
% volume of the source
75+
dvx = 1/((1-a)*(1-b));
76+
dvy = 1/((1-a)*b);
77+
dvz = 1/a;
4278

4379
dvmax = max([dvx,dvy,dvz]);
44-
dvz = w*dvz/dvmax;
45-
dvx = w*dvx/dvmax;
46-
dvy = w*dvy/dvmax;
47-
48-
% unitary patch Z
49-
pz = [-1,-1, 1, 1;
50-
-1, 1, 1,-1;
51-
0, 0, 0, 0]/2;
80+
dvz = w(n)*dvz/dvmax;
81+
dvx = w(n)*dvx/dvmax;
82+
dvy = w(n)*dvy/dvmax;
5283

53-
% unitary patch X
54-
px = [ 0, 0, 0, 0;
55-
-1, 1, 1,-1;
56-
1, 1,-1,-1]/2;
84+
if pplan
85+
% unitary patches (4 full squares)
86+
i4 = [0,-1,-1,0,0,0,-1,-1,0,1,1,0,0,0,1,1,0]/2;
87+
j4 = [0,0,-1,-1,0,1,1,0,0,0,1,1,0,-1,-1,0,0]/2;
88+
z4 = zeros(size(i4));
5789

58-
% unitary patch Y
59-
py = [-1, 1, 1,-1;
60-
0, 0, 0, 0;
61-
1, 1,-1,-1]/2;
90+
pz = [i4*dvx;j4*dvy;z4];
91+
px = [z4;i4*dvy;j4*dvz];
92+
py = [i4*dvx;z4;j4*dvz];
6293

63-
% applies rotations
64-
Rx = [1 0 0;0 cosd(ox) sind(ox);0 -sind(ox) cosd(ox)];
65-
Ry = [cosd(oy) 0 -sind(oy);0 1 0;sind(oy) 0 cosd(oy)];
66-
Rz = [cosd(oz) sind(oz) 0;-sind(oz) cosd(oz) 0;0 0 1];
67-
R = Rz*Ry*Rx;
94+
% applies rotations
95+
Rx = [1 0 0;0 cosd(ox) sind(ox);0 -sind(ox) cosd(ox)];
96+
Ry = [cosd(oy) 0 -sind(oy);0 1 0;sind(oy) 0 cosd(oy)];
97+
Rz = [cosd(oz) sind(oz) 0;-sind(oz) cosd(oz) 0;0 0 1];
98+
R = Rz*Ry*Rx;
6899

69-
px = R*(px*dvx);
70-
py = R*(py*dvy);
71-
pz = R*(pz*dvz);
100+
px = R*px;
101+
py = R*py;
102+
pz = R*pz;
72103

73-
% plots result
74-
switch lower(pview)
75-
case 'xy'
76-
h(1) = patch(pz(1,:) + x0,pz(2,:) + y0,'k',opt{:},varargin{:});
77-
h(2) = patch(px(1,:) + x0,px(2,:) + y0,'k',opt{:},varargin{:});
78-
h(3) = patch(py(1,:) + x0,py(2,:) + y0,'k',opt{:},varargin{:});
79-
case 'xz'
80-
h(1) = patch(pz(1,:) + x0,pz(3,:) + z0,'k',opt{:},varargin{:});
81-
h(2) = patch(px(1,:) + x0,px(3,:) + z0,'k',opt{:},varargin{:});
82-
h(3) = patch(py(1,:) + x0,py(3,:) + z0,'k',opt{:},varargin{:});
83-
case 'zy'
84-
h(1) = patch(pz(3,:) + z0,pz(2,:) + y0,'k',opt{:},varargin{:});
85-
h(2) = patch(px(3,:) + z0,px(2,:) + y0,'k',opt{:},varargin{:});
86-
h(3) = patch(py(3,:) + z0,py(2,:) + y0,'k',opt{:},varargin{:});
87-
88-
otherwise
89-
h(1) = patch(pz(1,:) + x0,pz(2,:) + y0, pz(3,:) + z0,'k',opt{:},varargin{:});
90-
h(2) = patch(px(1,:) + x0,px(2,:) + y0, px(3,:) + z0,'k',opt{:},varargin{:});
91-
h(3) = patch(py(1,:) + x0,py(2,:) + y0, py(3,:) + z0,'k',opt{:},varargin{:});
104+
% plots result
105+
switch lower(pview)
106+
case 'xy'
107+
h(1) = patch(pz(1,:) + x0,pz(2,:) + y0,dcolz,opt{:},varargin{:});
108+
h(2) = patch(px(1,:) + x0,px(2,:) + y0,dcolx,opt{:},varargin{:});
109+
h(3) = patch(py(1,:) + x0,py(2,:) + y0,dcoly,opt{:},varargin{:});
110+
case 'xz'
111+
h(1) = patch(pz(1,:) + x0,pz(3,:) + z0,dcolz,opt{:},varargin{:});
112+
h(2) = patch(px(1,:) + x0,px(3,:) + z0,dcolx,opt{:},varargin{:});
113+
h(3) = patch(py(1,:) + x0,py(3,:) + z0,dcoly,opt{:},varargin{:});
114+
case 'zy'
115+
h(1) = patch(pz(3,:) + z0,pz(2,:) + y0,dcolz,opt{:},varargin{:});
116+
h(2) = patch(px(3,:) + z0,px(2,:) + y0,dcolx,opt{:},varargin{:});
117+
h(3) = patch(py(3,:) + z0,py(2,:) + y0,dcoly,opt{:},varargin{:});
118+
otherwise
119+
h(1) = patch(pz(1,:) + x0,pz(2,:) + y0, pz(3,:) + z0,dcolz,opt{:},varargin{:});
120+
h(2) = patch(px(1,:) + x0,px(2,:) + y0, px(3,:) + z0,dcolx,opt{:},varargin{:});
121+
h(3) = patch(py(1,:) + x0,py(2,:) + y0, py(3,:) + z0,dcoly,opt{:},varargin{:});
122+
end
123+
124+
else
125+
[x,y,z] = ellipsoid3(dvx/2,dvy/2,dvz/2,ox,oy,oz,50);
126+
if nargout == 3
127+
varargout{1} = x;
128+
varargout{2} = y;
129+
varargout{3} = z;
130+
else
131+
h = surf(x + x0,y + y0,z + z0,varargin{:});
132+
end
92133
end
93-
94-
if nargout > 0
95-
varargout = h;
134+
135+
if nargout == 1
136+
varargout{1} = h;
96137
end
138+
139+
if ~hflag
140+
hold off
97141
end
142+
98143
end

0 commit comments

Comments
 (0)