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