The moment tensor in the above example has an equal plunges between P-and T-axis, so this issue of GMT might be relevant. But note that the problem also occurs in the mt convention, not in the principal axis representation (-Sy option in GMT or principal_axis convention in PyGMT).
I also noticed that we cound avoid the problem by adding a very small value (such as 0.0001 degree) to the rake angle.
def sdr2moment(strike, dip, rake, M0):
"""calculate moment tensor component from strike, dip, and rake
Reference: Aki and Richards (2002) Eq (1) of Box 4.4
"""
sinδ = np.sin( np.deg2rad(dip ))
cosδ = np.cos( np.deg2rad(dip ))
sin2δ = np.sin(2 * np.deg2rad(dip ))
cos2δ = np.cos(2 * np.deg2rad(dip ))
sinλ = np.sin( np.deg2rad(rake ))
cosλ = np.cos( np.deg2rad(rake ))
sinφ = np.sin( np.deg2rad(strike))
cosφ = np.cos( np.deg2rad(strike))
sin2φ = np.sin(2 * np.deg2rad(strike))
cos2φ = np.cos(2 * np.deg2rad(strike))
Mxx = - M0 * (sinδ * cosλ * sin2φ + sin2δ * sinλ * sinφ * sinφ )
Mxy = M0 * (sinδ * cosλ * cos2φ + sin2δ * sinλ * sin2φ / 2 )
Mxz = - M0 * (cosδ * cosλ * cosφ + cos2δ * sinλ * sinφ )
Myy = M0 * (sinδ * cosλ * sin2φ - sin2δ * sinλ * cosφ * cosφ )
Myz = - M0 * (cosδ * cosλ * sinφ - cos2δ * sinλ * cosφ )
Mzz = M0 * ( sin2δ * sinλ )
return Mxx, Myy, Mzz, Myz, Mxz, Mxy
def meca_aki2mt(strike, dip, rake, Mw):
M0 = 10**( 1.5 * (Mw + 10.7) ) # in dyn-cm unit
exponent = int(np.floor(np.log10(M0)))
# strike, dip, rake -> moment tensor in the cartesian coordinate
Mxx, Myy, Mzz, Myz, Mxz, Mxy = sdr2moment(strike, dip, rake, M0/10**exponent)
# from cartesian to polar coordinates
Mrr, Mtt, Mff, Mtf, Mrf, Mrt = Mzz, Mxx, Myy, -Mxy, -Myz, Mxz
spec = {"mrr": Mrr, "mtt": Mtt, "mff": Mff,
"mrf": Mrf, "mrt": Mrt, "mtf": Mtf,
"exponent": exponent}
return spec
## main script
fig = pygmt.Figure()
strike = 75
dip = 30
rake = 180
fig.basemap(projection='X10c/10c', region=[0, 4, 0, 4], frame = 'wsen')
fig.text(x=2, y=3.8, font="12p,Helvetica,Red", text="aki")
fig.text(x=3, y=3.8, font="12p,Helvetica,Blue", text="mt")
fig.text(x=1, y=3.8, font="10p,Helvetica,Black", text="(@~f@~,@~d@~,@~l@~)", justify='RM')
for i in range(3):
rake = 179 + i
fig.text(x=1, y=3-i, text=f"({strike}, {dip}, {rake})", justify = 'RM')
fig.meca(spec={"strike": strike, "dip": dip, "rake": rake, 'magnitude': 5},
longitude= 2, latitude=3-i, scale="2c", compressionfill='red')
fig.meca(spec = meca_aki2mt(strike, dip, rake, 5),
longitude= 3, latitude=3-i, scale="2c", compressionfill='blue')
fig.plot(x=[2.5, 3.5, 3.5, 2.5, 2.5], y=[1.5, 1.5, 2.5, 2.5, 1.5], pen='thicker,orange')
fig.show()No error message appeared.
Description of the problem
For certain combinations of strike ($\phi$ ), dip ($\delta$ ), and rake ($\lambda$ ) angles, the moment tensor-based expression (
mtconvention) of the beach ball bymecamodule gives wrongly swapped compression & expression fill. The following figure is a typical example of this problem, showing the mechanism solution with the rake angle changing by 1 degree.The
akiconvention does not cause this problem.The moment tensor in the above example has an equal plunges between P-and T-axis, so this issue of GMT might be relevant. But note that the problem also occurs in the
mtconvention, not in the principal axis representation (-Syoption in GMT orprincipal_axisconvention in PyGMT).I also noticed that we cound avoid the problem by adding a very small value (such as
0.0001degree) to the rake angle.Minimal Complete Verifiable Example
Full error message
System information