def METHOD_NAME(mesh_path, n, order):
    mesh = trimesh.load(mesh_path, process=False)
    vectors_orig, phi, theta = sampleSphericalDirections(n)
    SH_orig = getSHCoeffs(order, phi, theta)
    w = 4.0 * math.pi / (n * n)
    origins = mesh.vertices
    normals = mesh.vertex_normals
    n_v = origins.shape[0]
    origins = np.repeat(origins[:, None], n, axis=1).reshape(-1, 3)
    normals = np.repeat(normals[:, None], n, axis=1).reshape(-1, 3)
    PRT_all = None
    for i in tqdm(range(n)):
        SH = np.repeat(SH_orig[None, (i * n):((i + 1) * n)], n_v, axis=0).reshape(-1, SH_orig.shape[1])
        vectors = np.repeat(vectors_orig[None, (i * n):((i + 1) * n)], n_v, axis=0).reshape(-1, 3)
        dots = (vectors * normals).sum(1)
        front = (dots > 0.0)
        delta = 1e-3 * min(mesh.bounding_box.extents)
        hits = mesh.ray.intersects_any(origins + delta * normals, vectors)
        nohits = np.logical_and(front, np.logical_not(hits))
        PRT = (nohits.astype(np.float) * dots)[:, None] * SH
        if PRT_all is not None:
            PRT_all += (PRT.reshape(-1, n, SH.shape[1]).sum(1))
        else:
            PRT_all = (PRT.reshape(-1, n, SH.shape[1]).sum(1))
    PRT = w * PRT_all
    # NOTE: trimesh sometimes break the original vertex order, but topology will not change.
    # when loading PRT in other program, use the triangle list from trimesh.
    return PRT, mesh.faces