Set matrix rotation value (not add to it) 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?
Similarly to the answer to Convert mouse movement into rotation value (not in bge), customized to make the all selected objects rotate about the active object.

Need to subtract the pivot point of the rotation, which gives the "pivot arm" vector for the rotation.  In the previous answer, where just the active object was being rotated it had been hard coded to (0, 0, 0), since that is the result of subtracting an objects global location from itself.

import bpy, math
from bpy_extras.view3d_utils import (
        location_3d_to_region_2d,
        region_2d_to_vector_3d)

from mathutils import Matrix, Vector

class CustomRotate(bpy.types.Operator):
    bl_idname = "transform.custom_rotate"
    bl_label = "Custom Rotate"
    pivot_point = Vector((0, 0))
    v2d = Vector((0, 0))
    axis = Vector((1, 0, 0))
    mw = Matrix()

    def modal(self, context, event):
        ob = context.object
        if event.type in ['LEFTMOUSE','ENTER']:
            print('finished')
            return{'FINISHED'}

        if event.type in ['ESC','RIGHTMOUSE']: 
            ob.matrix_world = self.mw
            return {'CANCELLED'}

        if event.value == 'PRESS' and event.type == self.axis:
            # two times make it local
            axis = Vector([int(c == event.type) for c in 'XYZ'])
            self.axis = ob.matrix_world * axis  - ob.matrix_world * Vector()
            return {'RUNNING_MODAL'}                
        elif event.value == 'PRESS' and event.type in "XYZ":
            # to rotate about global axes
            self.axis = event.type
            # local
            return {'RUNNING_MODAL'}

        v2d = Vector((event.mouse_region_x, event.mouse_region_y)) - self.pivot_point
        if not v2d.length:
            return {'PASS_THROUGH'}
        #Get rotation angle and update object rotation
        a = self.v2d.angle_signed(v2d)
        rot_mat = Matrix.Rotation(a, 4, self.axis)
        t = ob.matrix_world.translation.copy()
        #Print angle in degrees to help with testing...
        #print('alpha angle:',math.degrees(a))
        for o in context.selected_objects:            
            mw = o.matrix_world
            mw.translation -=  t            
            mw = rot_mat * mw 
            mw.translation  += t
            o.matrix_world = mw

        self.v2d = v2d        

        return {'RUNNING_MODAL'}

    def invoke(self, context, event):        
        region = context.region
        rv3d = context.region_data

        ob = context.object
        self.pivot_point = location_3d_to_region_2d(region, rv3d, ob.matrix_world.translation)
        self.v2d = Vector((event.mouse_region_x, event.mouse_region_y)) - self.pivot_point
        z_axis = Vector((0, 0, -1))

        self.axis = rv3d.view_matrix.transposed().to_3x3() * z_axis
        self.mw = ob.matrix_world.copy()
        context.window_manager.modal_handler_add(self)

        return {'RUNNING_MODAL'}

def register():
    bpy.utils.register_class(CustomRotate)             
def unregister():
    bpy.utils.unregister_class(CustomRotate) 

if __name__ == "__main__":
    register()