I would use drivers for this. For a simple rotation around a static axis, it's quite straightforward, and the formulae can be entered directly in the driver window. No need for elaborate Python scripts.
A rotation quaternion is composed of four trigonometric values:
W = cos(a/2) for clockwise rotation or
W = -cos(a/2) for counter-clockwise rotation
X = dirX*sin(a/2)
Y = dirY*sin(a/2)
Z = dirZ*sin(a/2)
a is the rotation angle and
dirZ are the
Z values of a vector that describes the direction of the axis of rotation. This vector should have a magnitude of 1, i.e.
sqrt(dirX**2+dirY**2+dirZ**2) = 1
For linear rotation where 360° takes a known number of frames,
a/2 = (frame-1)*pi/(number of frames)
(Blender normally starts at frame 1, and the rotation at the first frame should be 0, hence the subtraction of 1 from the frame number).
In your example file the rotation is around the global Z axis, whose directional vector is
(0, 0, 1), and 360° takes 40 frames. This yields the following driver functions. Use the Scripted expression driver type.
W = cos((frame-1)*pi/40)
X = 0
Y = 0
Z = sin((frame-1)*pi/40)
In this image, I'm actually rotating around the Y axis in front view, but the basic idea is the same.
There are several methods to calculate the values of
dirZ. I will try to explain two of them.
You need to know the angles between the axis of rotation and the global X, Y and Z axes. Let's call these angles
aZ. Then it's as simple as
dirX = cos(aX)
dirY = cos(aY)
dirZ = cos(aZ)
These are the direction cosines of the axis, and if the angles are composed correctly, the vector
(dirX, dirY, dirZ) will have a magnitude of 1. If it doesn't, the angles define an axis that can't exist in 3D space.
Find the coordinates in the object's local space where the axis of rotation intersects the surface of the object. Let's call them
z. The vector
(x, y, z) will have a magnitude that may be 1, but quite possibly isn't. To normalize it, simply divide it by its own magnitude. Like this (dividing a vector by a scalar is the same as dividing its individual values by that scalar):
dirX = x/sqrt(x**2+y**2+z**2)
dirY = y/sqrt(x**2+y**2+z**2)
dirZ = z/sqrt(x**2+y**2+z**2)
A more in-depth explanation of the math behind all this, is most likely off topic here, but it probably would be right at home at math.SE.
I've found the following links quite useful:
Using Quaternion to Perform 3D rotations
Quaternions and spatial rotation