Matlab Script to decompose a moment tensor

mtensor_decomp.m

function mtensor_decomp
%
% function to decompose a moment-tensor
% 
% INPUT:
%        Moment tensor, M, in Aki & Richards format
%
%            Mnn Mne Mnz
%            Men Mee Mez
%            Mzn Mze Mzz
%
% M = [ Mnn  Mne Mzn;  Mne, Mee,  Mze; Mzn,  Mze, Mzz]
% M = [ MTT -MTF MRT; -MTF, MFF, -MRF; MRT, -MRF, MRR]
% replace the next line with your moment tensor (or make it an argument to the ftn)
M = [-0.1348, 0.5585,1.339;0.5585,-1.686,-3.165;1.339,-3.165,1.821]
%
% Test cases:
%
% M = [0,0,0;0,-1,0;0,0,1];    % North-striking 45 degree DS
% M = [0,1,0;1,0,0;0,0,0];     % North-Striking VSS
% M = [0,0,0;0,0,-1;0,-1,0];   % North-striking VDS
% M = [1,6,0;6,-2,-1;0,-1,4];  % Jost and Herrman, SRL Example

disp( '=======================================================' )
disp( '=======================================================' )
disp( '=======================================================' )
disp(sprintf('\nMOMENT-TENSOR Decomposition'));
%
disp(sprintf('\nThe moment tensor:\n'));
disp(M);
% compute the eigenvalues
%
[EF,DF] = eig(M);
disp(sprintf('\nEigenvalue Matrix:\n'));
disp(DF);
disp(sprintf('\nEigenvector Matrix:\n'));
disp(EF);
disp( '=======================================================' )
%--------------------------------------------------------------------------
l1 = DF(1:1,1:1);
l2 = DF(2:2,2:2);
l3 = DF(3:3,3:3);
a1 = EF(:,1);
a2 = EF(:,2);
a3 = EF(:,3);
%
moment = sqrt((l1*l1+l2*l2+l3*l3)/2.0);
disp( '=======================================================' )
disp(sprintf('Seismic Moment: %6.3f',moment));
disp( '=======================================================' )

%--------------------------------------------------------------------------
% remove the isotropic component
%--------------------------------------------------------------------------
iso = ( M(1:1,1:1) + M(2:2,2:2) + M(3:3,3:3)) / 3;
disp(sprintf('The Isotropic Component: %6.2f', iso));
disp( '=======================================================' )
l1 = l1 - iso;
l2 = l2 - iso;
l3 = l3 - iso;
%--------------------------------------------------------------------------
% sort the eigenvalues
%--------------------------------------------------------------------------
if(abs(l2) < abs(l1)) ,
	ltmp = l1;
	atmp = a1;
	l1 = l2;
	a1 = a2;
	l2 = ltmp;
	a2 = atmp;
end
 %-------------------------------------------------------------------------
 % check to see if #2 is less than #3
 %-------------------------------------------------------------------------
if(abs(l3) < abs(l2)) ,
	ltmp = l2;
	atmp = a2;
	l2 = l3;
	a2 = a3;
	l3 = ltmp;
	a3 = atmp;
 %-------------------------------------------------------------------------
 % check to see if the new # 2 is less than the new #1
 %-------------------------------------------------------------------------
	if(abs(l2) < abs(l1)) ,
		ltmp = l1;
		atmp = a1;
		l1 = l2;
		a1 = a2;
		l2 = ltmp;
		a2 = atmp;
	end
%--------------------------------------------------------------------------
end
disp(sprintf('The eigenvalues: %.3f %.3f %.3f',l1,l2,l3));
%--------------------------------------------------------------------------
% Write out epsilon (related to percent double couple
%--------------------------------------------------------------------------
disp( '=======================================================' )
disp(sprintf('Epsilon:') );
disp( '=======================================================' )
disp(sprintf('Herrmann & Jost Definition (positive only): %.2f\n',abs(l1 / l3)));
a0 = abs(DF(1:1,1:1));
if(abs(DF(3:3,3:3)) > a0) a0 = abs(DF(3:3,3:3));end
signed_epsilon = -DF(2:2,2:2)/a0;
disp(sprintf('Non-positive Definition (Giardini,1983): %.2f',signed_epsilon))
disp('(epsilon > 0 implies tension) ');
disp('(0 > epsilon implies compression)');
%--------------------------------------------------------------------------
% compute the major double couple
disp( '=======================================================' )
disp( 'MAJOR DOUBLE COUPLE' )
disp( '=======================================================' )
disp(sprintf( '\nMo: %.2f\n',abs(l3)));
%--------------------------------------------------------------------------
% plane one
%--------------------------------------------------------------------------
if(l3 > 0) ,
		p = a2;
		t = a3;
else
		p = a3;
		t = a2;
   end
disp(sprintf('P-axis: %6.3f %6.3f %6.3f',p));
disp(sprintf('T-axis: %6.3f %6.3f %6.3f\n',t));
%--------------------------------------------------------------------------
 u = 1 / sqrt(2) * (t - p);
 u = u / norm(u);
nu = 1 / sqrt(2) * (t + p);
nu = nu / norm(u);
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
disp( 'For First Plane:')
disp(sprintf('Slip Vector:   %6.3f %6.3f %6.3f',u));
disp(sprintf('Normal Vector: %6.3f %6.3f %6.3f',nu));
%
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
%
disp(sprintf('\nStrike     Dip      Rake'));
disp(sprintf('%5.1f     %4.1f     %5.1f',strike_deg, dip_deg,rake_deg));
%
%--------------------------------------------------------------------------
% plane two
%--------------------------------------------------------------------------
atmp = u;
 u = nu;
nu = atmp;
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
disp(sprintf('%5.1f     %4.1f     %5.1f\n',strike_deg, dip_deg,rake_deg));
%
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% compute the minor double couple
disp( '=======================================================' )
disp( 'MINOR DOUBLE COUPLE' )
disp( '=======================================================' )
mystring = sprintf( '\nMo: %.2f\n',abs(l1));
disp(mystring);
%
if (abs(l1) ~= 0 ) ,
%
%--------------------------------------------------------------------------
% plane one
%--------------------------------------------------------------------------
if(l1 > 0) ,
		p = a2;
		t = a1;
else
		p = a1;
		t = a2;
end
disp(sprintf('P-axis: %6.3f %6.3f %6.3f',p));
disp(sprintf('T-axis: %6.3f %6.3f %6.3f\n',t));
%--------------------------------------------------------------------------
 u = 1 / sqrt(2) * (t - p);
nu = 1 / sqrt(2) * (t + p);
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
disp( 'For First Plane:')
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
disp(sprintf('Slip Vector:   %6.3f %6.3f %6.3f',u));
disp(sprintf('Normal Vector: %6.3f %6.3f %6.3f',nu));
%
disp(sprintf('\nStrike     Dip      Rake'));
disp(sprintf('%5.1f     %4.1f     %5.1f',strike_deg, dip_deg,rake_deg));
%--------------------------------------------------------------------------
% plane two
%--------------------------------------------------------------------------
atmp = u;
 u = nu;
nu = atmp;
%
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
%
disp(sprintf('%5.1f     %4.1f     %5.1f\n',strike_deg, dip_deg,rake_deg));
%--------------------------------------------------------------------------
end
%
%
% compute the minor double couple preserving P axis
disp( '=======================================================' )
disp( 'MINOR DOUBLE COUPLE - Preserve P- or T-Axis' )
disp( '=======================================================' )
%
disp(sprintf( '\nMo: %.2f\n',abs(l1)));
%
if (abs(l1) ~= 0 ) ,
%
%--------------------------------------------------------------------------
% plane one
%--------------------------------------------------------------------------
if(l3 < 0) ,
		p = a3;
		t = a1;
else
		p = a1;
		t = a3;
   end
disp(sprintf('P-axis: %6.3f %6.3f %6.3f',p));
disp(sprintf('T-axis: %6.3f %6.3f %6.3f\n',t));
%--------------------------------------------------------------------------
 u = 1 / sqrt(2) * (t - p);
nu = 1 / sqrt(2) * (t + p);
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
%
disp( 'For First Plane:')
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
disp(sprintf('Slip Vector:   %6.3f %6.3f %6.3f',u));
disp(sprintf('Normal Vector: %6.3f %6.3f %6.3f',nu));
%
disp(sprintf('\nStrike     Dip      Rake'));
disp(sprintf('%5.1f     %4.1f     %5.1f',strike_deg, dip_deg,rake_deg));
%--------------------------------------------------------------------------
% plane two
%--------------------------------------------------------------------------
atmp = u;
 u   = nu;
nu   = atmp;
[strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu);
%
%  Make the values of strike and rake fall between -180 and +180
%
strike_rad = Nice_Angle(strike_rad);
   dip_rad = Nice_Angle(   dip_rad);
  rake_rad = Nice_Angle(  rake_rad);
%
%
strike_deg = 180 / pi * strike_rad;
dip_deg    = 180 / pi * dip_rad;
rake_deg   = 180 / pi * rake_rad;
disp(sprintf('%5.1f     %4.1f     %5.1f\n',strike_deg, dip_deg,rake_deg));
%--------------------------------------------------------------------------
end
%
%
%   SUB-FUNCTIONS CALLED BY THE FUNCTION ABOVE
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% make an angle lie between -pi and pi
%
function a = Nice_Angle(a)
%
if(a >= 2*pi) ,
	a = a - 2*pi;
end
%
if(a < -pi),
	a = a + 2 * pi;
end
if(a > pi),
	a = a - 2 * pi;
end
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%
% convert the slip and normal vectors to dislocation angles
%
function [strike_rad, dip_rad, rake_rad] = Axes_to_Angles(u,nu)
%
dip_rad    = acos(-nu(3:3));
strike_rad = atan2(nu(2:2),nu(1:1)) - pi / 2;
rake_rad   = asin(-u(3:3) / sin(dip_rad));
%
%  this is a check of the non-unique inversion equations
%
ltmp = cos(rake_rad)*sin(strike_rad)-cos(dip_rad)*cos(strike_rad)*sin(rake_rad);
%
dtmp = ltmp - u(2:2);
if(abs(dtmp) > 0.001) ,
  rake_rad = pi - rake_rad;
end
%
if ( dip_rad > pi / 2) ,
		strike_rad = pi + strike_rad;
		dip_rad = pi - dip_rad;
  rake_rad = - rake_rad;
end

Good morning

i want to ask you, what kind of seismic data format data that used to make moment tensor in those script? does it also relocate a microearthquake locations? and if it's allowed, may i have your email address for further informations?

thanks you,
Arif Setiawan