Moebius transformations in 3d

Moebius3d

//Moebius transformations in 3d, by reverse stereographic projection to the 3-sphere,
//rotation in 4d space, and projection back.
//by Daniel Piker 09/08/20
//Feel free to use, adapt and reshare. I'd appreciate a mention if you post something using this.
//You can also now find this transformation as a component in Grasshopper/Rhino
//I first wrote about these transformations here:
//https://spacesymmetrystructure.wordpress.com/2008/12/11/4-dimensional-rotations/

//If you want to transform about a given circle. Points on the circle and its axis stay on those curves.
//You can skip these 2 lines if you want to always use the origin centred unit circle.
Pt.Transform(Transform.PlaneToPlane(circle.Plane, Plane.WorldXY));
Pt /= circle.Radius;

double xa = Pt.X;
double ya = Pt.Y;
double za = Pt.Z;

double p = 0; //set to the same as q for isoclinic rotations
double q = 1; 

//reverse stereographic projection to hypersphere
double xb = 2 * xa / (1 + xa * xa + ya * ya + za * za);
double yb = 2 * ya / (1 + xa * xa + ya * ya + za * za);
double zb = 2 * za / (1 + xa * xa + ya * ya + za * za);
double wb = (-1 + xa * xa + ya * ya + za * za) / (1 + xa * xa + ya * ya + za * za);

//rotate hypersphere by amount t
double xc = +(xb) * (Math.Cos((p) * (t))) + (yb) * (Math.Sin((p) * (t)));
double yc = -(xb) * (Math.Sin((p) * (t))) + (yb) * (Math.Cos((p) * (t)));
double zc = +(zb) * (Math.Cos((q) * (t))) - (wb) * (Math.Sin((q) * (t)));
double wc = +(zb) * (Math.Sin((q) * (t))) + (wb) * (Math.Cos((q) * (t)));

//project stereographically back to flat 3D
double xd = xc / (1 - wc);
double yd = yc / (1 - wc);
double zd = zc / (1 - wc);

//the transformed point
Point3d P = new Point3d(xd, yd, zd);

//transform back to our input circle (you can skip this part if always using the origin centred unit circle)
P *= circle.Radius;
P.Transform(Transform.PlaneToPlane(Plane.WorldXY, circle.Plane));