三维空间随机旋转

echo "1
mol
H 0 0 5" | awk '
BEGIN{ Pi=4*atan2(1,1); TwoPi=2*Pi
	method=3; Nmax=1000
}

{	Natm=$1; getline
	for(i=1; i<=Natm; i++) { getline; Satm[i]=$1; Xatm[i]=$2; Yatm[i]=$3; Zatm[i]=$4 }

	print Nmax*Natm"\nRot. Method: "method
	for(n=1; n<=Nmax; n++) {
		isang=getAxi(method, axi, n, Nmax)

		if(isang) {
			tht=Pi
			x=.5*axi[1]
			y=.5*axi[2]
			z=.5*(axi[3]+1)
		} else {
			tht=axi[0]
			x=axi[1]
			y=axi[2]
			z=axi[3]
		}

		AxiRotMat(isang, tht, x, y, z)
		for(i=1; i<=Natm; i++) {
			V[1]=Xatm[i]; V[2]=Yatm[i]; V[3]=Zatm[i]
			RotVect(RotMat, V)
			printf "%s %9.6f %9.6f %9.6f\n", Satm[i], V[1], V[2], V[3]
		}
	}
	exit
}

function acos(x) { return atan2(sqrt(1-x*x),x) }

function getAxi(method, axi, Ipnt, Npnt,  i,r1,r2,tht1,tht2) {
	isang=1; ;axi[0]=0
	if(method==1) { # Marsaglia Method
		while(1) {
			r1=2*rand()-1
			r2=2*rand()-1
			if(r1*r1+r2*r2<1) break
		}
		tht1=r1*r1+r2*r2
		axi[1]=2*r1*sqrt(1-tht1)
		axi[2]=2*r2*sqrt(1-tht1)
		axi[3]=1-2*tht1
	} else if(method==2) {    # theta distributuion Method
		tht1=acos(2*rand()-1) # tht1=rand()*Pi # can NOT use
		tht2=rand()*TwoPi
		axi[1]=sin(tht1)*cos(tht2)
		axi[2]=sin(tht1)*sin(tht2)
		axi[3]=cos(tht1)
	} else if(method==3) { isang=0 # 四元数方法
		tht1=rand()*TwoPi
		tht2=rand()*TwoPi
		r1=sqrt(1-rand())
		r2=sqrt(1-r1*r1)
		axi[0]=sin(tht1)*r1
		axi[1]=cos(tht1)*r1
		axi[2]=sin(tht2)*r2
		axi[3]=cos(tht2)*r2
	} else if(method==4) { # Sprial Method, for regular equidistributed
		Ipnt--
		r1   = -1+2*(Ipnt+0.5)/Npnt # +0.5 表示在厚度的中点处取点
		tht2 = Pi*(sqrt(5)-1)*Ipnt
		r2   = sqrt(1-r1*r1)
		axi[1] = r2*cos(tht2)
		axi[2] = r2*sin(tht2)
		axi[3] = r1
	}
	return isang
}

function RotVect(RotMat, V) { # 按照矩阵RotMat旋转矢量V
	Vx = V[1]; Vy = V[2]; Vz = V[3]
	V[1]= RotMat[1,1]*Vx + RotMat[1,2]*Vy + RotMat[1,3]*Vz
	V[2]= RotMat[2,1]*Vx + RotMat[2,2]*Vy + RotMat[2,3]*Vz
	V[3]= RotMat[3,1]*Vx + RotMat[3,2]*Vy + RotMat[3,3]*Vz
}

function AxiRotMat(isang, tht, x, y, z,   w,rsq) { # 绕矢量(x,y,z)旋转tht弧度，四元数方法
	if(isang) { # w 为角度, 使用 轴-角 计算旋转矩阵
		w=cos(tht/2)
		rsq=sqrt(x*x+y*y+z*z); if(rsq>0) rsq = sin(tht/2)/rsq
		x *= rsq; y *= rsq; z *= rsq
	} else {    # 否则直接使用传入的四元数构建
		w=tht
		rsq=sqrt(x*x+y*y+z*z+w*w); if(rsq>0) rsq = 1/rsq
		x *= rsq; y *= rsq; z *= rsq; w *= rsq
	}
	RotMat[1,1]= 1-2*(y*y+z*z); RotMat[1,2]=   2*(x*y-w*z); RotMat[1,3]=   2*(x*z+w*y);
	RotMat[2,1]=   2*(x*y+w*z); RotMat[2,2]= 1-2*(x*x+z*z); RotMat[2,3]=   2*(y*z-w*x);
	RotMat[3,1]=   2*(x*z-w*y); RotMat[3,2]=   2*(y*z+w*x); RotMat[3,3]= 1-2*(x*x+y*y);
	if(x*x+y*y+z*z+w*w==0) { RotMat[1,1]=-1; RotMat[2,2]=-1; RotMat[3,3]=-1 }
}
'