You might want to take a look at the CelestialCoord class.  It might already do some of the spherical trig you're looking for here.  Or at least make it easier to do the calculation.  I'll also take a look and see if I agree with your math.

For your last point, I think the answer is yes.  That's the convention I adopted in the wcs module.  Namely that north is up and east is left in world coordinates.  If the world profile is defined in that way, then CelestialWCS classes (such as those normally used in fits files) will transform that correctly into image coordinates.  So we should keep that convention here for consistency.  And we should document it here as well to be clear.