That is great and very helpful.

What I did to get the axis perpendicular to two vectors is through solving two equations.

a1*x+b1*y+c1*z=0;
a2*x+b2*y+c2*z=0;

I assume z=1 and get y = (c1*a2-c2*a1)/(a1*b2-a2*b1); x = (c1*b2-c2*b1)/(b1*a2-b2*a1); and then handle the condition that z=0; That is awkward.

Now I get it and the calculation will be straght forward. I am glad I learn much from this thread.

One thing puzzles me.

The surface normal is calculated from vectorA and vectorB. If we swap vectorA and vectorB, the plane is the same, but the calculated surface normal will be a "negative" of the formal one.

I use that axis to rotate an angle so that vectorA may conform with vectorB. But I dont know whether I should rotate it Clockwise or counter-clockwise. (or which surface normal I should pick, the first one or the second one ?)