I know this question has been asked and answered before but the behavior I want is what Shift+Numpad 7 exactly does in edit mode which is to align the viewport to "face" the element selected (in my case a face or polygon) and "rotate" the viewport on its axis to match the orientation of the normal.

My code takes care of the first part which is to position the viewport to "face" the polygon and have no problem orienting itself on their normals except for faces that are rotated in multiple axis. e.g. the active face in the picture below, for the other faces which are angled only on a single axis the script has no problem orienting itself on their normals.

In my original code, the variables are delivered by raycast so I'd like to achieve this effect using only hit, normal plus the viewport and object properties. Face index is avoided since the target objects are using unapplied modifiers.

import bpy
import bmesh
from mathutils import *; from math import *

for area in bpy.context.screen.areas:
    if area.type == 'VIEW_3D':
       rv3d = area.spaces[0].region_3d
       
obj = bpy.context.edit_object
mat = obj.matrix_world
me = obj.data

bm = bmesh.from_edit_mesh(me)
f = bm.select_history.active

start = f.calc_center_median()
normal = f.normal

stop  = start + normal
norm = start - stop
norm.normalize()

q = norm.to_track_quat('-Z', 'Y')
rv3d.view_rotation = q
rv3d.view_location = start

bmesh.update_edit_mesh(me, True)

rotateview_blend_file