行为
支持 #464
打开多边形面积计算
状态:
新建
优先级:
普通
指派给:
-
目标版本:
-
开始日期:
2026-01-13
计划完成日期:
% 完成:
0%
预期时间:
#2:
描述
subroutine POLYICK ( NVERT, XYZ, AREA, IERRFL, TILT, I3DIM ) POLYI 77
c POLYI 78
c---- calculate area and coplanarity of a polygon. POLYI 79
c POLYI 80
c AREA .lt. 0 if Not CCW. ??? POLYI 81
c IERRFL =1 if NOT coplanar. POLYI 82
c =10+IERRFL if two vertices are identical. POLYI 83
c =100+IERRFL if 2-dim polygon is not CCW. POLYI 84
c TILT tilt in deg. POLYI 85
c I3DIM =1 if this is a 3d polygon. POLYI 86
c POLYI 87
real XYZ(3,NVERT) POLYI 88
c POLYI 89
IERRFL = 0 POLYI 90
c POLYI 91
c---- compute area , assume vertices are ordered ccw POLYI 92
c POLYI 93
XCOMP = 0.0 POLYI 94
YCOMP = 0.0 POLYI 95
ZCOMP = 0.0 POLYI 96
X0 = XYZ(1,NVERT) POLYI 97
Y0 = XYZ(2,NVERT) POLYI 98
Z0 = XYZ(3,NVERT) POLYI 99
c POLYI 100
do I = 1 , NVERT POLYI 101
X1 = XYZ(1,I) POLYI 102
Y1 = XYZ(2,I) POLYI 103
Z1 = XYZ(3,I) POLYI 104
XCOMP = XCOMP + Y0*Z1 - Y1*Z0 POLYI 105
YCOMP = YCOMP + Z0*X1 - Z1*X0 POLYI 106
ZCOMP = ZCOMP + X0*Y1 - X1*Y0 POLYI 107
X0 = X1 POLYI 108
Y0 = Y1 POLYI 109
Z0 = Z1 POLYI 110
enddo POLYI 111
AREA = 0.5 * sqrt( XCOMP*XCOMP + YCOMP*YCOMP + ZCOMP*ZCOMP ) POLYI 112
AREA = AREA/.092903 POLYI 113
if ( AREA .ne. 0.0 ) then POLYI 114
TILT = acos( ZCOMP/(2.*AREA) ) POLYI 115
PROJ = sqrt( XCOMP*XCOMP+YCOMP*YCOMP ) POLYI 116
POLYI 117
if ( PROJ .le. 0.0001*AREA ) then POLYI 118
AZIM = 0.0 POLYI 119
else POLYI 120
if ( XCOMP .lt. 0 ) then POLYI 121
if ( YCOMP .GE. 0 ) then POLYI 122
AZIM = 4.71238898 + asin(YCOMP/PROJ) POLYI 123
else POLYI 124
AZIM = 3.141592654 + asin(-XCOMP/PROJ) POLYI 125
endif POLYI 126
else POLYI 127
if( YCOMP .lt. 0 ) then POLYI 128
AZIM = 1.570796327 + asin(-YCOMP/PROJ) POLYI 129
else POLYI 130
AZIM = asin(XCOMP/PROJ) POLYI 131
endif POLYI 132
endif POLYI 133
endif POLYI 134
endif POLYI 135
c POLYI 136
c---- CK if all vertices are coplanar. POLYI 137
c POLYI 138
if ( NVERT .gt. 3 ) then POLYI 139
Z23 = XYZ(3,2) - XYZ(3,3) POLYI 140
Z13 = XYZ(3,1) - XYZ(3,3) POLYI 141
Z12 = XYZ(3,1) - XYZ(3,2) POLYI 142
A = XYZ(2,1) * Z23 - XYZ(2,2) * Z13 - XYZ(2,3) * Z12 POLYI 143
B = -( XYZ(1,1) * Z23 - XYZ(1,2) * Z13 - XYZ(1,3) * Z12 ) POLYI 144
C = XYZ(1,1) * ( XYZ(2,2) - XYZ(2,3) ) - POLYI 145
1 XYZ(1,2) * ( XYZ(2,1) - XYZ(2,3) ) + POLYI 146
2 XYZ(1,3) * ( XYZ(2,1) - XYZ(2,2) ) POLYI 147
D = -( XYZ(1,1) * ( XYZ(2,2)*XYZ(3,3) - XYZ(2,3)*XYZ(3,2) ) - POLYI 148
1 XYZ(2,1) * ( XYZ(1,2)*XYZ(3,3) - XYZ(1,3)*XYZ(3,2) ) + POLYI 149
2 XYZ(3,1) * ( XYZ(1,2)*XYZ(2,3) - XYZ(1,3)*XYZ(2,2) ) ) POLYI 150
do I = 4 , NVERT POLYI 151
if ( abs( A*XYZ(1,I) + B*XYZ(2,I) + C*XYZ(3,I) + D ) .gt. POLYI 152
1 0.001 ) IERRFL = 1 POLYI 153
enddo POLYI 154
endif POLYI 155
c POLYI 156
c---- CK if two vertices are identical. POLYI 157
c---- also set I3DIM POLYI 158
c POLYI 159
I3DIM = 0 POLYI 160
do I = 1 , NVERT POLYI 161
if ( XYZ(3,I) .ne. 0.0 ) I3DIM = 1 POLYI 162
if ( I .eq. NVERT ) goto 900 POLYI 163
if ( XYZ(1,I) .eq. XYZ(1,I+1) .and. POLYI 164
1 XYZ(2,I) .eq. XYZ(2,I+1) .and. POLYI 165
2 XYZ(3,I) .eq. XYZ(3,I+1) ) then POLYI 166
IERRFL = IERRFL + 10 POLYI 167
goto 900 POLYI 168
endif POLYI 169
enddo POLYI 170
900 continue POLYI 171
c POLYI 172
c---- for 2-dim polygons, CK if vertices are Not CCW. POLYI 173
c POLYI 174
if ( I3DIM .eq. 0 ) then POLYI 175
if ( XCOMP .lt. 0.0 .or. YCOMP .lt. 0.0 .or. POLYI 176
1 ZCOMP .lt. 0.0 ) IERRFL = IERRFL + 100 POLYI 177
endif POLYI 178
c POLYI 179
TILT = TILT * (180./3.141592654) POLYI 180
c POLYI 181
return POLYI 182
end POLYI 183
由 李立奎 更新于 6 天 之前
这是Fortran语言编写的 POLYICK 子程序，用于计算多边形面积并检查共面性。让我为您详细解释其工作原理：
功能概述¶
该子程序计算三维空间中多边形的：
- 面积（AREA）
- 倾斜角（TILT）- 多边形法向量与Z轴的夹角
- 共面性检查（IERRFL错误标志）
- 顶点顺序检查（是否逆时针CCW）
- 维度判断（I3DIM=1表示是3D多边形）
参数说明¶
subroutine POLYICK ( NVERT, XYZ, AREA, IERRFL, TILT, I3DIM )
-
NVERT: 顶点数量 -
XYZ: 顶点坐标数组，维度为XYZ(3,NVERT)，即 X=XYZ(1,I), Y=XYZ(2,I), Z=XYZ(3,I) -
AREA: 输出面积（单位：平方英尺，因为代码除以了0.092903转换因子） -
IERRFL: 错误标志位 -
TILT: 倾斜角（输出，单位：度） -
I3DIM: 输出标志，1=3D多边形，0=2D多边形
核心算法：面积计算（行93-113）¶
数学原理¶
使用多边形面积向量法（3D版格林公式）：
面积向量 = 0.5 × Σ(Vi × Vi+1) （Vi为顶点向量，最后一个顶点连接第一个）
代码实现¶
XCOMP = Σ(Y0×Z1 - Y1×Z0) ! 法向量的X分量
YCOMP = Σ(Z0×X1 - Z1×X0) ! 法向量的Y分量
ZCOMP = Σ(X0×Y1 - X1×Y0) ! 法向量的Z分量
AREA = 0.5 × √(XCOMP² + YCOMP² + ZCOMP²)
注意：(XCOMP, YCOMP, ZCOMP) 实际上是法向量的两倍（2 × 法向量），其模长等于 2 × AREA。
共面性检查（行139-155）¶
对于顶点数>3的多边形，通过平面方程验证：
-
计算平面系数（使用前3个点）：
Ax + By + Cz + D = 0代码中A,B,C,D的计算使用了行列式方法确定唯一平面。
-
验证所有点：检查第4个及之后的顶点是否满足
|Ax + By + Cz + D| ≤ 0.001- 若超出容差，设置
IERRFL = 1（不共面）
- 若超出容差，设置
错误标志位（IERRFL）解码¶
采用十进制位编码设计：
| 标志值 | 含义 |
|---|---|
1 |
多边形不共面 |
+10 |
存在两个相同顶点 |
+100 |
2D多边形顶点顺序不是逆时针（CCW） |
示例：
-
IERRFL = 101→ 不共面(1) + 非CCW(100) -
IERRFL = 11→ 不共面(1) + 重复顶点(10)
2D多边形逆时针检查（行175-178）¶
仅当 I3DIM=0（Z坐标全为0）时触发：
if (XCOMP < 0 或 YCOMP < 0 或 ZCOMP < 0) IERRFL += 100
原理：在2D情况下（XY平面），叉积的Z分量符号决定方向：
- ZCOMP > 0：逆时针（CCW）
- ZCOMP < 0：顺时针（CW）
倾斜角计算（行115, 180）¶
TILT = acos( ZCOMP / (2*AREA) ) ! 单位：弧度
TILT = TILT × (180/π) ! 转换为度
这是法向量与Z轴的夹角，即多边形的"倾斜"程度。
代码特点与注意事项¶
-
旧式Fortran风格：
- 使用
goto 900跳转 - 固定格式（每行末尾有行号如
POLYI 77） - 隐式类型声明（I-N开头的变量默认为整型）
- 使用
-
单位转换：
AREA = AREA/0.092903将面积从平方米转换为平方英尺（1 ft² = 0.092903 m²） -
重复顶点检查：仅检查相邻顶点是否完全相同
-
容差设置：
- 共面性检查：
0.001 - 投影判断：
0.0001 × AREA
- 共面性检查：
-
可能的bug：
TILT在未判断AREA ≠ 0的情况下计算，若面积为0可能导致除零错误（尽管前面有检查，但ACOS计算在ENDIF之后）。
流程总结¶
- 计算法向量和面积
- 计算倾斜角（弧度→度）
- 检查共面性（NVERT>3时）
- 检查重复顶点，设置I3DIM标志
- 2D情况下检查CCW顺序
- 返回综合错误标志
这是一个典型计算几何中的实用工具函数，常用于CAD、有限元前处理等工程软件。
行为