在管网空间信息系统的管线碰撞分析模块开发过程中，遇到计算异面直线的公垂线的问题。通过查找资料，得到了计算异面直线的公垂线距离的公式，但为了判断垂足是否在直线的确定的一段线段上，则必须计算出垂足的坐标。以下根据资料，演算了公垂线的垂足的坐标通用计算公式。
题目：已知两条异面直线AB、CD上两点坐标分别为A()、B()、C()、D()。推导直线AB与CD之间的公垂线在各自直线上的垂足坐标计算公式。
解题推导：
定义直线AB上的垂足M(Xm,Ym,Zm), 直线CD上的垂足N(Xn,Yn,Zn)。
定义 AM = t1*AB => (Xm-Xa,Y-Ya,Z-Za)=t1*(Xb-Xa,Yb-Ya,Zb-Za)
则 M坐标定义为
M( t1(Xb-Xa)+Xa , t1(Yb-Ya)+Ya , t1(Zb-Za)+Za );
同理 N坐标定义为
N( t2(Xd-Xc)+Xc , t2(Yd-Yc)+Yc , t2(Zd-Zc)+Zc );
所以 向量MN定义为 MN(t2(Xd-Xc)+Xc-t1(Xb-Xa)-Xa, t2(Yd-Yc)+Yc-t1(Yb-Ya)-Ya,
t2(Zd-Zc)+Zc-t1(Zb-Za)-Za);
因为 MN垂直于AB,MN垂直于CD，根据空间矢量点积特性，
(Xb-Xa)[t2(Xd-Xc)+Xc-t1(Xb-Xa)-Xa]+(Yb-Ya)[t2(Yd-Yc)+Yc-t1(Yb-Ya)-Ya]+(Zb-Za)[t2(Zd-Zc)+Zc-t1(Zb-Za)-Za]=0;
(Xd-Xc)[t2(Xd-Xc)+Xc-t1(Xb-Xa)-Xa]+(Yd-Yc)[t2(Yd-Yc)+Yc-t1(Yb-Ya)-Ya]+(Zd-Zc)[t2(Zd-Zc)+Zc-t1(Zb-Za)-Za]=0;
分别推导得到：
式1：
t2[(Xb-Xa)*(Xd-Xc)+(Yb-Ya)*(Yd-Yc)+(Zb-Za)*(Zd-Zc)]-t1[(Xb-Xa)* (Xb-Xa)+(Yb-Ya)*(Yb-Ya)+(Zb-Za)*(Zb-Za)]+[(Xb-Xa)*(Xc-Xa)+(Yb-Ya)*(Yc-Ya)+(Zb-Za)*(Zc-Za)]=0
式2：
t2[(Xd-Xc)*(Xd-Xc)+(Yd-Yc)*(Yd-Yc)+(Zd-Zc)*(Zd-Zc)]-t1[(Xb-Xa)*(Xd-Xc)+(Yb-Ya)*(Yd-Yc)+(Zb-Za)*(Zd-Zc)]+[(Xd-Xc)*(Xc-Xa)+(Yd-Yd)*(Yc-Ya)+(Zd-Zc)*(Zc-Za)]=0
令：
F1(a,b)=[(Xb-Xa)*(Xb-Xa)+(Yb-Ya)*(Yb-Ya)+(Zb-Za)*(Zb-Za)]
F1(c,d)= [(Xd-Xc)*(Xd-Xc)+(Yd-Yc)*(Yd-Yc)+(Zd-Zc)*(Zd-Zc)]
F2()=[(Xb-Xa)*(Xd-Xc)+(Yb-Ya)*(Yd-Yc)+(Zb-Za)*(Zd-Zc)]
F3(a,b)=[(Xb-Xa)*(Xc-Xa)+(Yb-Ya)*(Yc-Ya)+(Zb-Za)*(Zc-Za)]
F3(c,d)=[(Xd-Xc)*(Xc-Xa)+(Yd-Yc)*(Yc-Ya)+(Zd-Zc)*(Zc-Za)]
则：
式1=>式3：t2*F2()-t1*F1(a,b)+F3(a,b)=0;
式2=> 式4：t2*F1(c,d)-t1*F2()+F3(c,d)=0;
即：
t1=[F3(a,b)*F1(c,d)-F3(c,d)*F2()]/[F1(a,b)*F1(c,d)-F2()*F2()]
t2=[F3(c,d)*F1(a,b)-F2()*F3(a,b)]/[F2()*F2()-F1(a,b)*F1(c,d)]
由此得到两个垂足点的坐标：
M(Xm,Ym,Zm),
Xm=t1*(Xb-Xa)+Xa=(Xb-Xa)*[F3(a,b)*F1(c,d)-F3(c,d)*F2()]/[F1(a,b)*F1(c,d)-F2()*F2()]+Xa;
Ym=t1*(Yb-Ya)+Ya=(Yb-Ya)*[F3(a,b)*F1(c,d)-F3(c,d)*F2()]/[F1(a,b)*F1(c,d)-F2()*F2()]+Ya;
Zm=t1*(Zb-Za)+Za=(Zb-Za)*[F3(a,b)*F1(c,d)-F3(c,d)*F2()]/[F1(a,b)*F1(c,d)-F2()*F2()]+Za;
N(Xn,Yn,Zn),
Xn=t2*(Xd-Xc)+Xc=(Xd-Xc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Xc;
Yn=t2*(Yd-Yc)+Yc=(Yd-Yc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Yc;
Zn=t2*(Zd-Zc)+Zc=(Zd-Zc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Zc;