[Fortran] 纯文本查看 复制代码
Real Function planeeq(xyz) !通过三点来计算平面的方程，方程为aa*x+bb*y+cc*z+dd=0
Implicit None
Real,Parameter :: e = 0.0001 !设置误差范围
Real :: x1,y1,z1,x2,y2,z2,x3,y3,z3
Real, Dimension(3,3) :: xyz
Real, Dimension(4) :: planeeq
Real :: xnorm,det3
xyz(1,1)=x1
xyz(1,2)=y1
xyz(1,3)=z1
xyz(2,1)=x2
xyz(2,2)=y2
xyz(2,3)=z2
xyz(3,1)=x3
xyz(3,2)=y3
xyz(3,3)=z3
If (Abs(z1-z2)<= e .AND. Abs(z2-z3)<= e) then
aa = 0.0
bb = 0.0
cc = -1.0
dd = z1
Else if (Abs(x1-x2)<= e .AND. Abs(x2-x3)<= e) then
aa = -1.0
bb = 0.0
cc = 0.0
dd = x1
Else if (Abs(y1-y2)<= e .AND. Abs(y2-y3)<= e) then
aa = 0.0
bb = -1.0
cc = 0.0
dd = y1
Else
cc = 1.0
det3 = (x2-x1)*(y3-y2)-(x3-x2)*(y2-y1)
If (det3 <= e) then
cc = 0.0
aa = 1.0
If (y1 /= y2) then
bb = -1.0*(x2-x1)/(y2-y1)
dd = -1.0*x1-bb*y1
Else if (y2 /= y3) then
bb = -1.0*(x3-x2)/(y3-y2)
dd = -1.0*x2-bb*y2
Else if (y3 /= y1) then
bb = -1.0*(x1-x3)/(y1-y3)
dd = -1.0*x3-bb*y3
End if
Else
bb = -1.0*((x2-x1)*(z3-z2)-(x3-x2)*(z2-z1))/((x2-x1)*(y3-y2)-(x3-x2)*(y2-y1))
End if
If (x2 /= x1) then
aa = -1.0*(bb*(y2-y1)+cc*(z2-z1))/(x2-x1)
Else if (x3 /= x2) then
aa = -1.0*(bb*(y3-y2)+cc*(z3-z2))/(x3-x2)
End if
End if
xnorm = Sqrt(aa*aa+bb*bb+cc*cc)
aa = aa/xnorm
bb = bb/xnorm
cc = cc/xnorm
dd = dd/xnorm
End Function planeeq
下面我又定义了一个函数，需要调用函数planeeq计算的aa,bb,cc,dd，需要说明的是，这是在一个module中，aa,bb,cc,dd已经定义成了“全局变量”
Real Function calh(xyz_tetra) !来计算四面体的高
Implicit None
Real :: x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4
Real, Dimension(4,3) :: xyz_tetra
xyz_tetra(1,1)=x1
xyz_tetra(1,2)=y1
xyz_tetra(1,3)=z1
xyz_tetra(2,1)=x2
xyz_tetra(2,2)=y2
xyz_tetra(2,3)=z2
xyz_tetra(3,1)=x3
xyz_tetra(3,2)=y3
xyz_tetra(3,3)=z3
xyz_tetra(4,1)=x4
xyz_tetra(4,2)=y4
xyz_tetra(4,3)=z4
calh = Abs(aa*x4+bb*y4+dd*z4+dd)
End Function calh