You appear to be describing rotations of the Riemann Sphere: they are conformal and preserve the spherical metric. There are several simple ways to write them down (not involving any trigonometry!). A good one is to consider coordinates (x,y) in the plane as complex numbers z = x + yI, I^2 = -1. Given any four complex numbers a, b, c, and d, a fractional linear transformation of z is the value (az + b)/(cz + d). Thinking of z as the image of a point on the sphere via stereographic projection, it's easy to work out that for a fractional linear transformation to be a true rotation, a and d must be complex conjugates of each other and b and -c must also be mutual complex conjugates. (This defines the matrix group SU(2,C).)
Let Q be the point (in projected coordinates) you would like to base a new projection on. All we need to do is rotate Q to infinity, because infinity corresponds to the projection's base point (the north pole). This implies the denominator cz + d must equal zero, giving a solution
a = Conjugate of Q,
b = 1,
c = -1,
d = Q.
In other words, the change of coordinates you seek has the formula
z --> (1 + Conjugate(Q)*z) / (Q - z).
As an example, suppose you want to make Q = (5, 5) = 5 + 5I the new origin of projection. Using the rules of complex arithmetic we can work out the transformation in terms of the coordinates (x,y):
z --> (1 + (5 - 5I)*z) / (5 + 5I - z)
= (1 + (5 - 5I)*(x + yI)) / (5 + 5I - (x + yI))
= (1 + 5x + 5y + (5y - 5x)I) / (5 - x + (5-y)I)
= [(5 - x + 50y - 5x^2 - 5y^2) + (-5 - 50x + y + 5x^2 + 5y^2)I]
/ ((5-x)^2 + (5-y)^2).
That is, writing the new coordinates of (x,y) as (x',y'), we have
x' = (5 - x + 50y - 5x^2 - 5y^2) / ((5-x)^2 + (5-y)^2)
y' = (-5 - 50x + y + 5x^2 + 5y^2) / ((5-x)^2 + (5-y)^2).
You can follow this up by any rotation around the origin to reorient the new map. A good choice is to multiply the result by Q/Conjugate(Q) = Q^2 / |Q|^2. This is because the transformation has a nice interpretation:
(1 + Conjugate(Q)*z) / (Q - z) * Q / Conjugate(Q)
= (1/Conjugate(Q) + z) (1 + z/Q + (z/Q)^2 + ...).
The first factor merely translates all points by the (small) amount 1/Conjugate(Q). In particular, this keeps the map oriented correctly. The second factor can be ignored when z/Q is very small. Typically, for small rotations, Q is already "near infinity": that is, it is large compared with any point on the map, because it is the projection a point close to the original north pole. This justifies approximating the change in projection by means of a translation and, in addition, it tells us how to measure the error: the size of z/Q (the first, and largest, neglected term in the series on the right) is the ratio of the sizes of z and Q (that is, the ratio of their distances from the map origin). In other words, when the original projection of Q is way beyond the extent of your map, you will likely be ok with the approximate formula.