行为
支持 #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
行为