Ah.  Moller / Trumbore.  Right, the planes aren't pre-calculated for that.  Moller / Trumbore.  Hmm. 

This was actually faster than the... Moller / Trumbore.

def bary(rayOrg,rayDelta,n,p0,p1,p2,d): # Neil Hayes' awesome code - inputs rayOrg, rayDelta (@axis), trinormal, p0, p1, p2, polydistancefromplane
    #d is a plane, as in Geom.polyplanes
    """
    Adapted from jj_wrap.pl  version 2.2b (shrink-wrapping perlscript for Modo) by:
    julian johnson: julian@exch.demon.co.uk
    http://www.exch.demon.co.uk/jj_modoscripts.htm
    """
    ### Is ray pointing to front side of tri and not parallel
    ndot = n[0]*rayDelta[0] + n[1]*rayDelta[1] + n[2]*rayDelta[2] #dot(n,rayDelta)

    if not ndot < 0.0: return (-1, [0.0,0.0,0.0], ndot)

    ### compute parametric point of intersection
    t = d - (n[0]*rayOrg[0] + n[1]*rayOrg[1] + n[2]*rayOrg[2]) #dot(n,rayOrg)

    ### Is ray origin on backside of poly?
    if not t <= 0.0: return (-2, [0.0,0.0,0.0], t)

    ### Ray intersects plane so compute intersection point
    t /= ndot
    p = [0.0,0.0,0.0]
    p[0] = rayOrg[0] + rayDelta[0] * t
    p[1] = rayOrg[1] + rayDelta[1] * t
    p[2] = rayOrg[2] + rayDelta[2] * t

    ### Find dominant axis to project onto
    if abs(n[0]) > abs(n[1]):

        if abs(n[0]) > abs(n[2]):
            ##remove X
            u0 = p[1] - p0[1]
            u1 = p1[1] - p0[1]
            u2 = p2[1] - p0 [1]

            v0 = p[2] - p0[2]
            v1 = p1[2] - p0[2]
            v2 = p2[2] - p0[2]

        else:
            ##remove Z
            u0 = p[0] - p0[0]
            u1 = p1[0] - p0[0]
            u2 = p2[0] - p0[0]

            v0 = p[1] - p0[1]
            v1 = p1[1] - p0[1]
            v2 = p2[1] - p0[1]

    else:
        if abs(n[1]) > abs(n[2]):
            ## remove Y
            u0 = p[0] - p0[0]
            u1 = p1[0] - p0[0]
            u2 = p2[0] - p0[0]

            v0 = p[2] - p0[2]
            v1 = p1[2] - p0[2]
            v2 = p2[2] - p0[2]

        else:
            ## remove Z
            u0 = p[0] - p0[0]
            u1 = p1[0] - p0[0]
            u2 = p2[0] - p0 [0]

            v0 = p[1] - p0[1]
            v1 = p1[1] - p0[1]
            v2 = p2[1] - p0[1]

    ### computer denominator for 2D projection
    temp = (u1v2) - (v1u2)
    if not (temp != 0.0):
        return (-3, [0.0,0.0,0.0], temp) #try next poly

    temp = 1.0/temp

    ### compute each barycentric (check for out of range)
    alpha = (u0v2 - v0u2) * temp
    if not (alpha >= 0.0):
        return (-4, [0.0,0.0,0.0], alpha) #try next poly

    beta = (u1v0 - v1u0) * temp
    if not (beta >= 0.0):
        return (-5, [0.0,0.0,0.0], beta) #try next poly

    gamma = 1.0 - alpha - beta
    if not (gamma >= 0.0):
        return (-6, [0.0,0.0,0.0], gamma) #try next poly

    #if made here then ray hit a target tri
    return (1, p, t)

===========================sigline======================================================

Cage can be an opinionated jerk who posts without thinking.  He apologizes for this.  He's honestly not trying to be a turkeyhead.

Cage had some freebies, compatible with Poser 11 and below.  His Python scripts were saved at archive.org, along with the rest of the Morphography site, where they were hosted.