Bug: RectangularPatches.getpatchgeometry returns strike off by 180°

[ecaballerol]

Summary

In RectangularPatches.py, getpatchgeometry computes:

strike = np.arctan2(vs[0], vs[1]) - np.pi
if strike < 0.:
    strike += 2*np.pi

The - np.pi term makes every patch's returned strike exactly 180° from its true value, regardless of dip direction or corner numbering. The previous version (no - np.pi) was correct.

Reproducer : using planarfault.buildPatches directly

Using planarfault's actual patch-construction code (discretize + buildPatches), across a full 3×2 patch grid:

import numpy as np

def discretize(strike, length, n_strike):
    strike_rad = strike*np.pi/180.
    xc, yc = 0., 0.
    half_length = 0.5*length
    x0 = xc - half_length * np.sin(strike_rad)
    y0 = yc - half_length * np.cos(strike_rad)
    dist_strike = np.linspace(0, length, n_strike+1)
    xi = x0 + dist_strike * np.sin(strike_rad)
    yi = y0 + dist_strike * np.cos(strike_rad)
    return xi, yi

def buildPatches(strike, dip, f_length, f_width, n_strike, n_dip, top=0.):
    p_width = f_width/float(n_dip)
    xi, yi = discretize(strike, f_length, n_strike)
    zi = np.ones(xi.shape)*top
    dip_rad = dip*np.pi/180.
    dipdirection_rad = ((strike + 90) % 360) * np.pi/180.
    patches = []
    for i in range(n_dip):
        xt, yt, zt = xi, yi, zi
        xb = xt + p_width * np.cos(dip_rad) * np.sin(dipdirection_rad)
        yb = yt + p_width * np.cos(dip_rad) * np.cos(dipdirection_rad)
        zb = zt + p_width*np.sin(dip_rad)
        for j in range(xt.shape[0]-1):
            x1,y1,z1 = xt[j], yt[j], zt[j]
            x2,y2,z2 = xt[j+1], yt[j+1], zt[j+1]
            x3,y3,z3 = xb[j+1], yb[j+1], zb[j+1]
            x4,y4,z4 = xb[j], yb[j], zb[j]
            if y1 > y2:
                p2 = [x1,y1,z1]; p1 = [x2,y2,z2]; p4 = [x3,y3,z3]; p3 = [x4,y4,z4]
            else:
                p1 = [x1,y1,z1]; p2 = [x2,y2,z2]; p3 = [x3,y3,z3]; p4 = [x4,y4,z4]
            patches.append(np.array([p1,p2,p3,p4]))
        xi, yi, zi = xb, yb, zb
    return patches

def getpatchgeometry_strike(patch, subtract_pi):
    p1, p2, p3, p4 = patch
    vs, vd = p2 - p1, p4 - p1
    if vs[1]*vd[0] - vs[0]*vd[1] < 0.:
        vs = vs * -1
    strike = np.arctan2(vs[0], vs[1]) - (np.pi if subtract_pi else 0)
    if strike < 0.:
        strike += 2*np.pi
    return np.degrees(strike)

for true_strike in [10, 90, 200, 300]:
    patches = buildPatches(true_strike, dip=40, f_length=10, f_width=10, n_strike=3, n_dip=2)
    for idx, p in enumerate(patches):
        old = getpatchgeometry_strike(p, subtract_pi=False)
        new = getpatchgeometry_strike(p, subtract_pi=True)
        print(true_strike, idx, old, new)

Every one of the 6 patches per fault reproduces the same result: old matches the input strike exactly, current is off by 180° — confirming this isn't an artifact of a hand-built patch, and that planarfault's own corner-assignment rule (if y1>y2: swap) doesn't change the outcome.

Impact

strike feeds directly into Green's function construction (e.g. EDKSmp.dropSourcesInPatches), so this flips the along-strike convention on every patch — not a scaling error, a fixed rotation per patch. Any class inheriting RectangularPatches without overriding getpatchgeometry (e.g. RectangularPatchesKin) is affected too.

Suggested fix

strike = np.arctan2(vs[0], vs[1])
if strike < 0.:
    strike += 2*np.pi

The vnz check already above it handles both corner-winding cases correctly, so nothing further is needed. (For what it's worth, the commented-out alternative np.arctan2(-normal[0], normal[1]) - np.pi left in the code isn't a drop-in fix either — it doesn't recover the true strike reliably.)