Matlab画震源球

2018-04-03

今天有人问我会不会用matlab画震源球，他画出来的是扁的。给他推荐GMT或者Python画完后用AI修图软件给加到底图上去，但需要把成千上百个震源球投影到matlab画的一张底图上，不太好弄。

Matlab画震源球代码

搜索到matlab函数bb，该函数可以画double-couple 和全矩张量,代码如下

function bb(fm, centerX, centerY, diam, ta, color)
%function bb(fm, centerX, centerY, diam, ta, color)
%
%	Draws beachball diagram of earthquake double-couple focal mechanism(s). S1, D1, and
%		R1, the strike, dip and rake of one of the focal planes, can be vectors of
%		multiple focal mechanisms.
%	fm - focal mechanism that is either number of mechnisms (NM) by 3 (strike, dip, and rake)
%		or NM x 6 (mxx, myy, mzz, mxy, mxz, myz - the six independent components of
%		the moment tensor). The strike is of the first plane, clockwise relative to north.
%		The dip is of the first plane, defined clockwise and perpedicular to strike, 
%		relative to horizontal such that 0 is horizontal and 90 is vertical. The rake is 
%		of the first focal plane solution. 90 moves the hanging wall up-dip (thrust),
%		0 moves it in the strike direction (left-lateral), -90 moves it down-dip
%		(normal), and 180 moves it opposite to strike (right-lateral).
%	centerX - place beachball(s) at position centerX
%	centerY - place beachball(s) at position centerY
%	diam - draw with this diameter.  If diam is zero, beachball
%		is drawn on stereonet.
%	ta - type of axis. If 0, this is a normal axis, if 1 it is a map axis. In case
%		of the latter, centerX and centerY are Lons and Lats, respectively.
%	color - color to use for quadrants of tension; can be a string, e.g. 'r'
%		'b' or three component color vector, [R G B].
%

% eg: figure(1);clf;bb([20,30,40],0,0,10,0,'r');axis equal;set(gca,'Visible','off');

[ne,n] = size(fm);
if n == 6
	for j = 1:ne
		[s1(j),d1(j),r1(j)] = mij2sdr(fm(j,1),fm(j,2),fm(j,3),fm(j,4),fm(j,5),fm(j,6));
	end
else
	s1 = fm(:,1);
	d1 = fm(:,2);
	r1 = fm(:,3);
end

r2d = 180/pi;
d2r = pi/180;
ampy = cos(mean(centerY)*d2r);

if ne > 1
	[ds,i] = sort(diam,1,'descend');
	diam = diam(i);
	s1 = s1(i);
	d1 = d1(i);
	r1 = r1(i);
	centerX = centerX(i);
	centerY = centerY(i);
end

mech = zeros(ne,1);
j = find(r1 > 180);
r1(j) = r1(j) - 180;
mech(j) = 1;
j = find(r1 < 0);
r1(j) = r1(j) + 180;
mech(j) = 1;

% Get azimuth and dip of second plane
[s2,d2,r2] = AuxPlane(s1,d1,r1);

if diam(1) > 0
	hold on
end
for ev = 1:ne
	S1 = s1(ev);
	D1 = d1(ev);
	S2 = s2(ev);
	D2 = d2(ev);
	P = r1(ev);
	CX = centerX(ev);
	CY = centerY(ev);
	D = diam(ev);
	M = mech(ev);

if M > 0
   P = 2;
else
   P = 1;
end

if D1 >= 90
   D1 = 89.9999;
end
if D2 >= 90
   D2 = 89.9999;
end

phi = 0:.01:pi;
d = 90 - D1;
m = 90;
l1 = sqrt(d^2./(sin(phi).^2 + cos(phi).^2 * d^2/m^2));

d = 90 - D2;
m = 90;
l2 = sqrt(d^2./(sin(phi).^2 + cos(phi).^2 * d^2/m^2));

if D == 0
   stereo(phi+S1*d2r,l1,'k')
   hold on
   stereo(phi+S2*d2r,l2,'k')
end

inc = 1;
[X1,Y1] = pol2cart(phi+S1*d2r,l1);
if P == 1
   lo = S1 - 180;
   hi = S2;
   if lo > hi
      inc = -inc;
   end
   th1 = S1-180:inc:S2;
   [Xs1,Ys1] = pol2cart(th1*d2r,90*ones(1,length(th1)));
   [X2,Y2] = pol2cart(phi+S2*d2r,l2);
   th2 = S2+180:-inc:S1;
else
   hi = S1 - 180;
   lo = S2 - 180;
   if lo > hi
      inc = -inc;
   end
   th1 = hi:-inc:lo;
   [Xs1,Ys1] = pol2cart(th1*d2r,90*ones(1,length(th1)));
   [X2,Y2] = pol2cart(phi+S2*d2r,l2);
   X2 = fliplr(X2);
   Y2 = fliplr(Y2);
   th2 = S2:inc:S1;
end
[Xs2,Ys2] = pol2cart(th2*d2r,90*ones(1,length(th2)));

X = cat(2,X1,Xs1,X2,Xs2);
Y = cat(2,Y1,Ys1,Y2,Ys2);

if D > 0
   X = ampy*X * D/90 + CY;
   Y = Y * D/90 + CX;
   phid = 0:.01:2*pi;
   [x,y] = pol2cart(phid,90);
   xx = x*D/90 + CX;
   yy = ampy*y*D/90 + CY;
   if ta == 0
      fill(xx,yy,'w')
      fill(Y,X,color)
      line(xx,yy,'color','k','linewidth',0.2);
   else
      fillm(yy,xx,'w')
      fillm(X,Y,color)
      linem(yy,xx,'color','k','linewidth',0.2);
   end
else
   if ta == 0
      fill(X,Y,color)
   else
      fillm(Y,X,color)
   end
   view(90,-90)
end

end

function [strike, dip, rake] = AuxPlane(s1,d1,r1);
%function [strike, dip, rake] = AuxPlane(s1,d1,r1);
% Get Strike and dip of second plane, adapted from Andy Michael bothplanes.c
r2d = 180/pi;

z = (s1+90)/r2d;
z2 = d1/r2d;
z3 = r1/r2d;
%/* slick vector in plane 1 */
sl1 = -cos(z3).*cos(z)-sin(z3).*sin(z).*cos(z2);
sl2 = cos(z3).*sin(z)-sin(z3).*cos(z).*cos(z2);
sl3 = sin(z3).*sin(z2);
[strike, dip] = strikedip(sl2,sl1,sl3);

n1 = sin(z).*sin(z2);  %/* normal vector to plane 1 */
n2 = cos(z).*sin(z2);
n3 = cos(z2);
h1 = -sl2; %/* strike vector of plane 2 */
h2 = sl1;
%/* note h3=0 always so we leave it out */

z = h1.*n1 + h2.*n2;
z = z./sqrt(h1.*h1 + h2.*h2);
z = acos(z);

rake = zeros(size(strike));
j = find(sl3 > 0);
rake(j) = z(j)*r2d;
j = find(sl3 <= 0);
rake(j) = -z(j)*r2d;

function [strike, dip] = strikedip(n, e, u)
%function [strike, dip] = strikedip(n, e, u)
%       Finds strike and dip of plane given normal vector having components n, e, and u
%

% Adapted from Andy Michaels stridip.c

r2d = 180/pi;

j = find(u < 0);
n(j) = -n(j);
e(j) = -e(j);
u(j) = -u(j);

strike = atan2(e,n)*r2d;
strike = strike - 90;
while strike >= 360
        strike = strike - 360;
end
while strike < 0
        strike = strike + 360;
end

x = sqrt(n.^2 + e.^2);
dip = atan2(x,u)*r2d;

function hpol = stereo(theta,rho,line_style)
%function hpol = stereo(theta,rho,line_style)

if nargin < 1
    error('Requires 2 or 3 input arguments.')
elseif nargin == 2 
    if isstr(rho)
        line_style = rho;
        rho = theta;
        [mr,nr] = size(rho);
        if mr == 1
            theta = 1:nr;
        else
            th = (1:mr)';
            theta = th(:,ones(1,nr));
        end
    else
        line_style = 'auto';
    end
elseif nargin == 1
    line_style = 'auto';
    rho = theta;
    [mr,nr] = size(rho);
    if mr == 1
        theta = 1:nr;
    else
        th = (1:mr)';
        theta = th(:,ones(1,nr));
    end
end
if isstr(theta) | isstr(rho)
    error('Input arguments must be numeric.');
end
if ~isequal(size(theta),size(rho))
    error('THETA and RHO must be the same size.');
end

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');
set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state

% make a radial grid
    hold on;
    maxrho = max(abs(rho(:)));
    hhh=plot([-maxrho -maxrho maxrho maxrho],[-maxrho maxrho maxrho -maxrho]);
    set(gca,'dataaspectratio',[1 1 1],'plotboxaspectratiomode','auto')
    set(gca,'xlim',[-90 90])
    set(gca,'ylim',[-90 90])
    v = [get(cax,'xlim') get(cax,'ylim')];
    ticks = sum(get(cax,'ytick')>=0);
    delete(hhh);
% check radial limits and ticks
    rmin = 0; rmax = v(4); rticks = max(ticks-1,2);
%    if rticks > 5   % see if we can reduce the number
%        if rem(rticks,2) == 0
%            rticks = rticks/2;
%        elseif rem(rticks,3) == 0
%            rticks = rticks/3;
%        end
%    end
	 rticks = 1;

% define a circle
    th = 0:pi/50:2*pi;
    xunit = cos(th);
    yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
    inds = 1:(length(th)-1)/4:length(th);
    xunit(inds(2:2:4)) = zeros(2,1);
    yunit(inds(1:2:5)) = zeros(3,1);
% plot background if necessary
    if ~isstr(get(cax,'color')),
       patch('xdata',xunit*rmax,'ydata',yunit*rmax, ...
             'edgecolor',tc,'facecolor',get(gca,'color'),...
             'handlevisibility','off');
    end

% draw radial circles
    c82 = cos(82*pi/180);
    s82 = sin(82*pi/180);
    rinc = (rmax-rmin)/rticks;
    for i=(rmin+rinc):rinc:rmax
        hhh = plot(xunit*i,yunit*i,ls,'color',tc,'linewidth',1,...
                   'handlevisibility','off');
%        text((i+rinc/20)*c82,(i+rinc/20)*s82, ...
%            ['  ' num2str(i)],'verticalalignment','bottom',...
%            'handlevisibility','off')
    end
    set(hhh,'linestyle','-') % Make outer circle solid

% plot spokes
    th = (1:6)*2*pi/12;
    cst = cos(th); snt = sin(th);
%    cs = [-cst; cst];
%    sn = [-snt; snt];
%    plot(rmax*cs,rmax*sn,ls,'color',tc,'linewidth',1,...
%         'handlevisibility','off')

% annotate spokes in degrees
    rt = 1.1*rmax;
    for i = 1:length(th)
        text(rt*cst(i),rt*snt(i),int2str(i*30),...
             'horizontalalignment','center',...
             'handlevisibility','off');
        if i == length(th)
            loc = int2str(0);
        else
            loc = int2str(180+i*30);
        end
        text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center',...
             'handlevisibility','off')
    end

% set view to 2-D
    view(2);
% set axis limits
    axis(rmax*[-1 1 -1.15 1.15]);
end

% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits',fUnits );

% transform data to Cartesian coordinates.
xx = rho.*cos(theta);
yy = rho.*sin(theta);

% plot data on top of grid
if strcmp(line_style,'auto')
    q = plot(xx,yy);
else
    q = plot(xx,yy,line_style);
end
if nargout > 0
    hpol = q;
end
if ~hold_state
    set(gca,'dataaspectratio',[1 1 1]), axis off; set(cax,'NextPlot',next);
end
set(get(gca,'xlabel'),'visible','on')
set(get(gca,'ylabel'),'visible','on')


function [str,dip,rake] = mij2sdr(mxx,myy,mzz,mxy,mxz,myz)
%function [str,dip,rake] = mij2sdr(mxx,myy,mzz,mxy,mxz,myz)
%
%   INPUT
%	mij - siz independent components of the moment tensor
%
%   OUTPUT
%	str - strike of first focal plane (degrees)
%	dip - dip of first focal plane (degrees)
%	rake - rake of first focal plane (degrees)
%
%Adapted from code, mij2d.f, created by Chen Ji and given to me by Gaven Hayes.
%

a = [mxx mxy mxz; mxy myy myz; mxz myz mzz];
[V,d] = eig(a);

D = [d(3,3) d(1,1) d(2,2)];
V(2:3,1:3) = -V(2:3,1:3);
V = [V(2,3) V(2,1) V(2,2); V(3,3) V(3,1) V(3,2); V(1,3) V(1,1) V(1,2)];

IMAX = find(D == max(D));
IMIN = find(D == min(D));
AE = (V(:,IMAX)+V(:,IMIN))/sqrt(2.0);
AN = (V(:,IMAX)-V(:,IMIN))/sqrt(2.0);
AER = sqrt(AE(1)^2+AE(2)^2+AE(3)^2);
ANR = sqrt(AN(1)^2+AN(2)^2+AN(3)^2);
AE = AE/AER;
AN = AN/ANR;
if (AN(3) <= 0.)
	AN1 = AN;
	AE1 = AE;
else
	AN1 = -AN;
	AE1 = -AE;
end 
[ft,fd,fl] = TDL(AN1,AE1);
str = 360 - ft;
dip = fd;
rake = 180 - fl;

function [FT,FD,FL] = TDL(AN,BN)
XN=AN(1);
YN=AN(2);
ZN=AN(3);
XE=BN(1);
YE=BN(2);
ZE=BN(3);
AAA=1.0E-06;
CON=57.2957795;
if (abs(ZN) < AAA)
	FD=90.;
	AXN=abs(XN);
	if (AXN > 1.0) 
		AXN=1.0;
	end
	FT=asin(AXN)*CON;
	ST=-XN;
	CT=YN;
	if (ST >= 0. & CT < 0) 
		FT=180.-FT;
	end
	if (ST < 0. & CT <= 0) 
		FT=180.+FT;
	end
	if (ST < 0. & CT > 0) 
		FT=360.-FT;
	end
	FL=asin(abs(ZE))*CON;
	SL=-ZE;
	if (abs(XN) < AAA) THEN
		CL=XE/YN;
	else
		CL=-YE/XN;
	end 
	if (SL >= 0. & CL < 0) 
		FL=180.-FL;
	end
	if (SL < 0. & CL <= 0) 
		FL=FL-180.;
	end
	if (SL < 0. & CL > 0) 
		FL=-FL;
	end
else
	if (-ZN > 1.0) 
		ZN=-1.0;
	end
	FDH=acos(-ZN);
	FD=FDH*CON;
	SD=sin(FDH);
	if  (SD == 0)
		return;
	end 
	ST=-XN/SD;
	CT=YN/SD;
	SX=abs(ST);
	if (SX > 1.0) 
		SX=1.0;
	end
	FT=asin(SX)*CON;
	if (ST >= 0. & CT < 0) 
		FT=180.-FT;
	end
	if (ST < 0. & CT <= 0) 
		FT=180.+FT;
	end
	if (ST < 0. & CT > 0) 
		FT=360.-FT;
	end
	SL=-ZE/SD;
	SX=abs(SL);
	if (SX > 1.0) 
		SX=1.0;
	end
	FL=asin(SX)*CON;
	if (ST == 0) THEN
		CL=XE/CT;
	else
		XXX=YN*ZN*ZE/SD/SD+YE;
		CL=-SD*XXX/XN;
		if (CT == 0) 
			CL=YE/ST;
		end
	end 
	if (SL >= 0. & CL < 0) 
		FL=180.-FL;
	end
	if (SL < 0. & CL <= 0) 
		FL=FL-180.;
	end
	if (SL < 0. & CL > 0) 
		FL=-FL;
	end
end

在matlab下执行下列命令
figure(1);
clf;
bb([20,30,40],0,0,10,0,'r');
axis equal; % 少了 axis equal 球就是扁的了。
set(gca,'Visible','off'); %去掉坐标轴

其他画震源球的方法(不全)

FocalMech,但只能画全矩张量。
Obspy的函数以及命令行式的obspy-mopad，该命令可以做分量分解以及画图(推荐).

补充

double-couple 转 moment tensor

建议用 sdr2mt转成mt后，调用FocalMech画图，如果叠加到底图上，注释掉FocalMech文件中的colormap

感谢秀文的技术支持


function [mt]=sdr2mt(strike,dip,rake)
%SDR2MT    Convert strike-dip-rake to moment tensor
%
%    Usage:    mt=sdr2mt(sdr)
%              mt=sdr2mt(strike,dip,rake)
%
%    Description:
%     MT=SDR2MT(SDR) converts a double-couple (aka focal mechanism) given
%     in strike-dip-rake form (Nx3) as [strike dip rake] to Harvard moment
%     tensor form (Nx6) as [Mrr Mtt Mpp Mrt Mrp Mtp] which uses the Up,
%     South, East coordinate system.  N allows for multiple focal
%     mechanisms to be converted simultaneously.  The main purpose of this
%     function is to allow plotting of the focal mechanism with PLOTMT.
%     Strike is positive clockwise from North, dip is positive downward
%     from the horizontal and rake is positive counter-clockwise in the
%     fault plane from the strike direction.  Note that the strike must be
%     such that when you look along the direction of the strike the fault
%     dips to your right.
%
%     MT=SDR2MT(STRIKE,DIP,RAKE) allows inputing the strike, dip & rake
%     separately.  Note that the inputs should all be Nx1 column vectors or
%     scalars.
%
%    Notes:
%     - Tested OK against George Helffrich's website:
%        http://www1.gly.bris.ac.uk/~george/focmec.html
%     - See Aki & Richards (2002) Figure 4.20 & Box 4.4 for details.
%
%    Examples:
%     % Plot a dip-slip fault with EW strike and 45 deg dip:
%     plotmt(0,0,sdr2mt(90,45,-90))
%     axis tight equal off
%
%     % Validate SDR2MT & MT2SDR using GlobalCMT catalog:
%     cmts=findcmts;
%     sdr1=[cmts.strike1 cmts.dip1 cmts.rake1];
%     sdr2=[cmts.strike2 cmts.dip2 cmts.rake2];
%     sdr3=mt2sdr(sdr2mt(sdr1));
%     max(min(abs(azdiff(sdr3,sdr1)),abs(azdiff(sdr3,sdr2))))
%     sdr4=mt2sdr(sdr2mt(sdr2));
%     max(min(abs(azdiff(sdr4,sdr1)),abs(azdiff(sdr4,sdr2))))
%
%    See also: MT2SDR, AUXPLANE, PLOTMT, MT_DECOMP, MT_DIAG, MT_UNDIAG

%     Version History:
%        Mar.  8, 2010 - initial version
%        June  1, 2011 - harvard output, now with docs
%        Mar. 21, 2013 - doc update, clean up code, no mt cmp output
%        Mar. 25, 2013 - minor doc update
%
%     Written by Garrett Euler (ggeuler at wustl dot edu)
%     Last Updated Mar. 25, 2013 at 13:50 GMT

% todo:

% check nargin

% convert strike-dip-rake to HRV mt
% f = strike, d = dip, l = rake
%1 Mrr =  Mzz =  Mo sin2d sinl
%2 Mtt =  Mxx = -Mo(sind cosl sin2f +     sin2d sinl (sinf)^2 )
%3 Mpp =  Myy =  Mo(sind cosl sin2f -     sin2d sinl (cosf)^2 )
%4 Mrt =  Mxz = -Mo(cosd cosl cosf  +     cos2d sinl sinf )
%5 Mrp = -Myz =  Mo(cosd cosl sinf  -     cos2d sinl cosf )
%6 Mtp = -Mxy = -Mo(sind cosl cos2f + 0.5 sin2d sinl sin2f )

mt=nan(numel(strike),6);
mt(:,1)=sind(2*dip).*sind(rake);
mt(:,2)=sind(dip).*cosd(rake).*sind(2*strike) ...
    + sind(2*dip).*sind(rake).*sind(strike).^2;
mt(:,3)=sind(dip).*cosd(rake).*sind(2*strike) ...
    - sind(2*dip).*sind(rake).*cosd(strike).^2;
mt(:,4)=cosd(dip).*cosd(rake).*cosd(strike) ...
    + cosd(2*dip).*sind(rake).*sind(strike);
mt(:,5)=cosd(dip).*cosd(rake).*sind(strike) ...
    - cosd(2*dip).*sind(rake).*cosd(strike);
mt(:,6)=sind(dip).*cosd(rake).*cosd(2*strike) ...
    + 0.5.*sind(2*dip).*sind(rake).*sind(2*strike);
mt(:,2:2:6)=-mt(:,2:2:6);

end

function focalmech(fm, centerX, centerY, diam, varargin)
%	>> focalmech(fm, centerX, centerY, diam, varargin)
%   Draws full moment tensor beachball diagram of earthquake focal mechanism.
%   Required inputs
%	fm: 3x3 or 1x6 vector of the six independent components of the 
%           moment tensor (mrr, mtt, mff, mrt, mrf, mtf).
%	centerX, centerY: position to place beachball at position
%	diam: diameter of drawn beachball.
%
%	Optional inputs
%	color for tensional region changed by 'r', 'g', 'b', 'y', 'm', 'c' or 
%       3 element color vector. default is black.
%   scalar argument is assumed to be an aspect ratio for non-equal axes.
%       Value is stretch in E-W direction.
%       'map' or 'Mercator' will do automatic Mercator aspect ratio
%       assuming centerY = lat
%   'dc' will force to best fitting double-couple solution
%   'nofill' will leave transparent (nodal lines only)
%   'text' will add following argument as a title to the beachball
%   'FontSize' with following argument will change fontsize for text
%
%   Written to give a MatLab equivalent of psmeca -Sm from GMT.
%   As such, this function requires a moment tensor. For planes (e.g.,
%   strikes and dips), use bb.m. bb is a nice beachball plotting function
%   from Oliver Boyd at CERI, but aimed strictly at double-couple
%   mechanisms.
%
%   Written by James Conder, Southern Illinois University, Jan 10, 2017
%   Please cite:
%   Conder, J.A. and C.A. Arciniegas, Conjugate Faulting in the Wabash
%       Valley Fault Zone exhibited by the Nov 20, 2012 M3.6 earthquake, 
%       a Mt Carmel Late Aftershock, Seismological Research Letters, 88,
%       1203-1209, doi:10.1785/0220170021, 2017
%
%   2/20/2017 realized that a vertical vector may not actually be within
%   a single contour if there is a large isotropic component.
%   Algorithm changed to using sign of vector explicitly inside contour.
%   3/12/2017 Vectorized main loop to improve efficiency.
%   (Speed gains of 4-10x).

% defaults
fillcolor = [0 0 0];  % black for tensional region(s)
unitratio = 1;      % axis equal
DC = false;         % show non double-couple components
colorit = true;     % fill in tensional regions
ctitle = [];        % text above beachball
fsize = 10;         % fontsize for text

% set parameters based on varargin inputs
if nargin > 4
    icheck(1:nargin-4) = true;
    for i = 1:nargin-4
        if icheck(i)
        if isnumeric(varargin{i})
            if isscalar(varargin{i})
                unitratio = varargin{i};
            end
            if length(varargin{i}) == 3
                fillcolor = varargin{i};
            end
        elseif ischar(varargin{i})
            switch varargin{i}
                case {'r'}
                    fillcolor = [ 1 0 0 ];
                case {'g'}
                    fillcolor = [ 0 1 0 ];
                case {'b'}
                    fillcolor = [ 0 0 1 ];
                case {'y'}
                    fillcolor = [ 1 1 0 ];
                case {'m'}
                    fillcolor = [ 1 0 1 ];
                case {'c'}
                    fillcolor = [ 0 1 1 ];
                
                case {'map','Map','mercator','Mercator'}
                    unitratio = 1/cos(centerY*pi/180);
                    
                case {'dc','DC','double-couple','doublecouple'}
                    DC = true;

                case {'nocolor','nofill','transparent'}
                    colorit = false;

                case {'text','title'}
                    ctitle = varargin{i+1};
                    icheck(i+1) = false;

                case {'fontsize','FontSize','Fontsize','fsize'}
                    fsize = varargin{i+1};
                    icheck(i+1) = false;
            end 
        end
        end
    end     
end

%colormap([1 1 1; fillcolor])    % 2 color beachball fill   % 注释掉，默认为黑白球
sE = unitratio*0.5*diam; sN = 0.5*diam; % scaling for plotting
u = sE*cos(0:0.02:2*pi); w = sN*sin(0:0.02:2*pi);   % reference circle


%%% put fm (A&R convention, rtf) into 3x3 cartesian (xyz)
%%% x = north, y = East, Z = down (Aki & Richards pg 113)
M = eye(3);
if length(fm) == 6
   M(1,1) = fm(2); M(2,2) = fm(3); M(3,3) = fm(1);
   M(2,1) = -fm(6); M(1,2) = M(2,1);
   M(3,1) = fm(4); M(1,3) = M(3,1);
   M(3,2) = -fm(5); M(2,3) = M(3,2);
else
   M(1,1) = fm(2,2); M(2,2) = fm(3,3); M(3,3) = fm(1,1);
   M(2,1) = -fm(3,2); M(1,2) = M(2,1);
   M(3,1) = fm(2,1); M(1,3) = M(3,1);
   M(3,2) = -fm(3,1); M(2,3) = M(3,2);
end
    
% Tensional and compressional regions are determined by sign of 
%   local 'radius'.
% Negative length = compression, and positive length = tension
% Length is found by linear combination of three eigen vectors.

% find eigenvalues and eigenvectors of system and put in descending order
[V,D] = eig(M);
D = diag(D);
eig1 = max(D);      % largest eigenvalue
eig3 = min(D);      % smallest eigenvalue
if eig1 == eig3
    eig2 = eig1;
    i1 = 1; i2 = 2; i3 = 3;
else
    i1 = find(D == eig1,1);
    i3 = find(D == eig3,1);
    i2 = 6 - i1 - i3;
    eig2 = D(i2);
end

vT = V(:,i1);       % eigenvectors: TBP
vB = V(:,i2);
vP = V(:,i3);

% check for explosions or implosions
if eig1 <= 0     % implosion
    fill(centerX+u,centerY+w,'w')   % fill white circle
    return
end
if eig3 >= 0     % explosion
    fill(centerX+u,centerY+w,fillcolor)   % filled circle
    return
end

%%% if forced to be double couple
if DC
    eigB = 0.5*(eig1 - eig3);
    eig1 = eigB; eig2 = 0; eig3 = -eigB;
end

% Sweep out hemisphere to find negative & positive regions
dx = 0.02; dy = dx;     % grid loop
x = -1:dx:1; y = -1:dy:1;
nx = length(x); ny = length(y); % vectorization of previous code begins here
vij = zeros(3,nx,ny);

x2 = repmat(x',1,ny);
y2 = repmat(y,nx,1);
r2 = x2.*x2 + y2.*y2;        
trend = atan2(y2,x2);
plunge = pi/2 - 2*asin(sqrt(r2/2));  % equal area projection
ir = r2 > 1;
vij(1,:,:) = cos(trend).*cos(plunge);   % set up local vector grids
vij(2,:,:) = sin(trend).*cos(plunge);
vij(3,:,:) = sin(plunge);

% project eigenvectors onto local vectors
uT = 0*vij; uB = 0*vij; uP = 0*vij; 
uT(1,:,:) = vij(1,:,:)*vT(1) + vij(2,:,:)*vT(2) + vij(3,:,:)*vT(3);
uT(3,:,:) = uT(1,:,:)*vT(3);
uT(2,:,:) = uT(1,:,:)*vT(2);
uT(1,:,:) = uT(1,:,:)*vT(1);
uB(1,:,:) = vij(1,:,:)*vB(1) + vij(2,:,:)*vB(2) + vij(3,:,:)*vB(3);
uB(3,:,:) = uB(1,:,:)*vB(3);
uB(2,:,:) = uB(1,:,:)*vB(2);
uB(1,:,:) = uB(1,:,:)*vB(1);
uP(1,:,:) = vij(1,:,:)*vP(1) + vij(2,:,:)*vP(2) + vij(3,:,:)*vP(3);
uP(3,:,:) = uP(1,:,:)*vP(3);
uP(2,:,:) = uP(1,:,:)*vP(2);
uP(1,:,:) = uP(1,:,:)*vP(1);

% get weights for each eigenvector across hemisphere space
wT = uT(1,:,:).*uT(1,:,:) + uT(2,:,:).*uT(2,:,:) + uT(3,:,:).*uT(3,:,:);
wB = uB(1,:,:).*uB(1,:,:) + uB(2,:,:).*uB(2,:,:) + uB(3,:,:).*uB(3,:,:);
wP = uP(1,:,:).*uP(1,:,:) + uP(2,:,:).*uP(2,:,:) + uP(3,:,:).*uP(3,:,:);

uz = wT*eig1 + wB*eig2 + wP*eig3;   % lengths grid
uz = squeeze(uz); uz(ir) = nan;
clear vij uT uB uP wT wB wP

%%% plot
% Note centerX & centerY are in conventional cartesian coords unlike x,y
% above where x is north.
hold on
if ~colorit
    plot(centerX+u,centerY+w,'k')   % plot circle boundary
    c = contourc(centerX+sE*y,centerY+sN*x,uz,[0 0]);       % no fill
    i = c(2,:) >= 1;
    c(:,i) = nan;
    plot(centerX+sE*c(1,:),centerY+sN*c(2,:),'color',fillcolor)   % plot circle boundary    
else
    fill(centerX+u,centerY+w,'w')   % fill white circle
    plot(centerX+u,centerY+w,'k')   % plot circle boundary


% fill in tensional regions
%contourf(centerX+sE*y,centerY+sN*x,uz,[0 0]);   % fill in tensional regions
% Need to use 'fill' as contourf will change pallette of entire figure

% get contours of zero level. *Should* be one or two segments
c = contourc(y,x,uz,[0 0]); c = c';

% find contour breaks
nc = 1;
[~,I] = sort(c(:,2),'descend');
if c(I(2),2) > 7
    nc = 2;
end
I = I(1:nc);

%reset c with just contours of interest (maximum of 2)
if nc == 1
    c = [ c(I(1),:); c(I(1)+1:I(1)+c(I(1),2),:)];
else
    c = [ c(I(1),:); c(I(1)+1:I(1)+c(I(1),2),:); ...
          c(I(2),:); c(I(2)+1:I(2)+c(I(2),2),:)];
end
i = find(c(:,2) > 2);  % indices of contour breaks

% fill tensional regions
if nc == 1       % 1 contour (closed unless single vertical plane)
    c(1,:) = [];        % remove contour header 
    % closed contour must have vertical plunge in it - check for sign
    % not necessarily true if large isotropic component - fix
    
    % check for single vertical plane
    d2 = (c(1,1)-c(end,1))^2 + (c(1,2)-c(end,2))^2;
    
    if d2 > 2 % single vertical plane
        ang = atan2(c(end,2),c(end,1));     % location of contour end 
        dang = 0.1;
        vij = [ sin(ang+dang) cos(ang+dang) 0];
        uT = dot(vij,vT)*vT; wT = norm(uT)^2;
        uB = dot(vij,vB)*vB; wB = norm(uB)^2;  
        uP = dot(vij,vP)*vP; wP = norm(uP)^2;
        uvc = wT*eig1 + wB*eig2 + wP*eig3;
        if uvc > 0       
            a1 = (ang:0.02:ang+pi)';    % ccw along edge
        else
            a1 = flipud((ang-pi:0.02:ang)');    % cw along edge
        end   
        ac = [cos(a1) sin(a1)];
        c1 = [c; ac];
        fill(centerX+sE*c1(:,1),centerY+sN*c1(:,2),fillcolor)  
    else        % single closed contour
        xi = mean(c(:,1)); yi = mean(c(:,2)); % spot inside contour
        trend = atan2(xi,yi);   % x = north...     
        r2 = xi*xi + yi*yi;        
        plunge = pi/2 - 2*asin(sqrt(r2/2));  % equal area projection        
        vij(1) = cos(trend)*cos(plunge);
        vij(2) = sin(trend)*cos(plunge);
        vij(3) = sin(plunge);

        uT = dot(vij,vT)*vT; wT = norm(uT)^2;
        uB = dot(vij,vB)*vB; wB = norm(uB)^2;  
        uP = dot(vij,vP)*vP; wP = norm(uP)^2;
        uv = wT*eig1 + wB*eig2 + wP*eig3;
    
        if uv > 0   % fill inside contour
            fill(centerX+sE*c(:,1),centerY+sN*c(:,2),fillcolor)
        elseif uv < 0   % fill outside contour
            fill(centerX+u,centerY+w,fillcolor)
            fill(centerX+sE*c(:,1),centerY+sN*c(:,2),'w') 
        end
    end
else                    % 2 contours
    % find 4 spots where contours touch edge
    ia = ([ 2 i(2)-1 i(2)+1 length(c)])';
    dang = 0.02*pi;
    
    % determine whether contours bound positive or negative regions
    % get mid contour vector
    p1 = c(3,:);
    p2 = c(ia(2)-1,:);
    pmid = 0.5*(p1 + p2);
    ang = atan2(pmid(2),pmid(1));
    vij = [ sin(ang) cos(ang) 0]; 
    uT = dot(vij,vT)*vT; wT = norm(uT)^2;
    uB = dot(vij,vB)*vB; wB = norm(uB)^2;  
    uP = dot(vij,vP)*vP; wP = norm(uP)^2;
    uvc = wT*eig1 + wB*eig2 + wP*eig3;
        
    % first contour ending spot
    p2 = c(ia(2),:);
    ang2 = atan2(p2(2),p2(1));

    % angle between contour ends
    p1 = c(2,:);
    p2 = c(ia(2),:);
    ang = acos(dot(p1,p2)/(norm(p1)*norm(p2)));
    
    % determine direction to trace along edge
    d1 = acos(dot(vij,[sin(ang2+dang) cos(ang2+dang) 0]));
    d2 = acos(dot(vij,[sin(ang2-dang) cos(ang2-dang) 0]));
    if d1 < d2
        a1 = (ang2+0.01:0.02:ang2+ang-0.01)';    % ccw along edge
    else
        a1 = flipud((ang2-ang+0.01:0.02:ang2-0.01)');    % cw along edge
    end
        
    ac = [cos(a1) sin(a1)];
    c1 = [c(2:ia(2),:); ac; c(2,:)];
    if uvc < 0
        fill(centerX+u,centerY+w,fillcolor)   % fill colored circle
        fill(centerX+sE*c1(:,1),centerY+sN*c1(:,2),[ 1 1 1]) 
    else
        fill(centerX+sE*c1(:,1),centerY+sN*c1(:,2),fillcolor)
    end
     
    % second contour ending spot
    p2 = c(end-1,:);
    ang2 = atan2(p2(2),p2(1));
    
    % mid contour vector
    p1 = c(ia(3)+1,:);
    p2 = c(ia(end)-1,:);
    pmid = 0.5*(p1 + p2);
    ang = atan2(pmid(2),pmid(1));
    vij = [ sin(ang) cos(ang) 0]; 

    % angle between contour ends
    p1 = c(ia(3),:);
    p2 = c(end,:);
    ang = acos(dot(p1,p2)/(norm(p1)*norm(p2)));
    
    % determine direction to trace along edge
    d1 = acos(dot(vij,[sin(ang2+dang) cos(ang2+dang) 0]));
    d2 = acos(dot(vij,[sin(ang2-dang) cos(ang2-dang) 0]));
    if d1 < d2
        a1 = (ang2+0.01:0.02:ang2+ang-0.01)';    % ccw along edge
    else
        a1 = flipud((ang2-ang+0.01:0.02:ang2-0.01)');    % cw along edge
    end    

    ac = [cos(a1) sin(a1)];
    c2 = [c(ia(3):end,:); ac; c(ia(3),:)];
    if uvc < 0
        fill(centerX+sE*c2(:,1),centerY+sN*c2(:,2),'w')        
    else
        fill(centerX+sE*c2(:,1),centerY+sN*c2(:,2),fillcolor)
    end        

end
end

if ~isempty(ctitle)
    text(centerX,centerY+sN*1.1,ctitle,'HorizontalAlignment','center','FontSize',fsize)
end

end