The below code will rotate a selection of objects around an arbitrary location (using the active object's location here) by the chosen angle (45 in this example). I cobbled this code together from a few places but I don't fully understand how to manipulate matrices very well yet.

If you create a few objects and press run script you'll see 45 degrees is added to the objects while transforming them around the active object. This is almost exactly what I need minus one issue. I need this code to set the rotation value to 45 (or whatever arbitrary value) when the script is run not add 45 degrees to the selection of object's matrices.

import bpy
from math import radians
from mathutils import Matrix

ob = bpy.context.object
rot_point = ob.location
rot_angle = 45
rot_axis = 'Y'

mat = (Matrix.Translation(rot_point) *
       Matrix.Rotation(radians(rot_angle), 4, rot_axis) *
       Matrix.Translation(-rot_point))


for obj in bpy.context.selected_objects:
    obj.matrix_world = mat * obj.matrix_world

EDIT: If you create a few random objects, select them and then run the script, each time you run the script 45 degrees (45 degrees is an arbitrary value, it means nothing, just a random input value for this example) will be added to the object rotation values (as displayed in the clip below). This is good, and I need the objects to move like they're moving. Ultimately the value controlled by rot_angle will be fed into from mouse coordinates in a modal though, so I don't need the code to add the rot_angle value to these object's values, I need to be able to set them. Ultimately, the code should rotate the objects around whatever chosen rot_point in the manner displayed below, however once rotated, pressing run script shouldn't do anything. As the code stands, when fed into a modal it breaks everything because the code continually adds to the values with the mouse movement instead of setting them.

UPDATE: The below code sort of does what I want. The problem is it starts breaking down once the object gets rotated in to many ways. If I could subtract the rotation matrix values to get the same results for the rotation angle difference (ra_dif) instead of comparing against the Euler rotation value of the object, I'm guessing it would work better. How could I do that? Or does someone see a better (less dense) way of accomplishing what I'm trying to do here?

import bpy
from math import radians,degrees
from mathutils import Matrix

ob = bpy.context.object
rot_point = bpy.context.scene.cursor_location
rot_angle = 45
rot_axis = 'Y'

  
ra_dif = rot_angle - round(degrees(ob.rotation_euler.y))

mat = (Matrix.Translation(rot_point) *
       Matrix.Rotation(radians(ra_dif), 4, rot_axis) *
       Matrix.Translation(-rot_point))
   

for obj in bpy.context.selected_objects:
    obj.matrix_world = mat * obj.matrix_world

UPDATE 2: The code posted by batFINGER works great with one issue. Seems to be some strange sliding going on when objects have been rotated by a few odd degrees. If you open up the following file, press Run Script, then press space and access "Custom Rotate" then immediately press Z to enter global rotation, you'll notice the selected objects will slide if you pull the mouse in towards the pivot point of the rotation. What could be causing this? How could it be corrected?