pygmt Swapped compression/extension fills in meca module with mt convention

Description of the problem

For certain combinations of strike ($\phi$), dip ($\delta$), and rake ($\lambda$) angles, the moment tensor-based expression (mt convention) of the beach ball by meca module 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 aki convention does not cause this problem.

swapped_mech

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.

Minimal Complete Verifiable Example

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()

Full error message

No error message appeared.

System information

PyGMT information:
  version: v0.13.0
System information:
  python: 3.12.5 | packaged by conda-forge | (main, Aug  8 2024, 18:36:51) [GCC 12.4.0]
  executable: /home/tktmyd/miniforge3/envs/seismo24/bin/python
  machine: Linux-6.8.0-45-generic-x86_64-with-glibc2.35
Dependency information:
  numpy: 1.26.4
  pandas: 2.2.2
  xarray: 2024.9.0
  netCDF4: 1.7.1
  packaging: 24.1
  contextily: None
  geopandas: None
  IPython: 8.27.0
  rioxarray: None
  gdal: 3.9.2
  ghostscript: 10.03.1
GMT library information:
  version: 6.5.0
  padding: 2
  share dir: /home/tktmyd/miniforge3/envs/seismo24/share/gmt
  plugin dir: /home/tktmyd/miniforge3/envs/seismo24/lib/gmt/plugins
  library path: /home/tktmyd/miniforge3/envs/seismo24/lib/libgmt.so
  cores: 36
  grid layout: rows
  image layout: 
  binary version: 6.5.0

Dec 12 '24 03:12 tktmyd

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

Dec 12 '24 03:12 welcome[bot]

@seisman, could you check if this is the same issue upstream as https://github.com/GenericMappingTools/gmt/issues/8009? I'm not a seismologist, but would guess that this has to do with some numerical precision error in GMT's code for meca.