subroutine aptpoly (ax, ay, az, bx, by, bz, cx, cy, cz, np, tol,
& dx, dy, dz, sl, ri, rc, ap, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPOLY
c
c call aptpoly (ax, ay, az, bx, by, bz, cx, cy, cz, np, tol,
c & dx, dy, dz, sl, ri, rc, ap, nerr)
c
c Version: aptpoly Updated 1990 September 19 17:40.
c aptpoly Originated 1989 September 19 17:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the np vertices d = (dx, dy, dz) of the regular polygon
c of np sides, with the center at point a = (ax, ay, az), the
c first vertex at point b = (bx, by, bz), and in the plane
c passing through the points "a", "b" and c = (cx, cy, cz).
c Also, to find the edge length sl, the radius of the inscribed
c circle ri, the radius of the circumscribed circle rc, and the
c area of the polygon ap.
c Flag nerr indicates any input error.
c
c Note: Other subroutines may be used to translate (apttran),
c rotate (aptrota, aptrotp, aptrots, aptrott, aptrotv),
c or scale (aptsclu) the regular polygon generated by aptpoly.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, np, tol.
c
c Output: dx, dy, dz, sl, ri, rc, ap, nerr.
c
c Calls: aptvpln, aptrota, aptmopv
c
c
c Glossary:
c
c ap Output The area of the polygon.
c ap = 0.5 * np * ri * sl.
c ap = 0.5 * np * rc**2 * sin (2.0 * pi / np).
c
c ax,ay,az Input The x, y and z coordinates of the center of the
c regular polygon with np sides.
c
c bx,by,bz Input The x, y and z coordinates of the first vertex of the
c regular polygon with np sides. Point "b" must be
c distinct from point "a", based on tol.
c
c cx,cy,cz Input The x, y and z coordinates of any point distinct from,
c and not colinear with, points "a" and "b", based on
c tol, and in the plane of the regular polygon.
c If a vector "v" normal to the plane of the polygon
c is known (the dot product of "v" with vector "ab"
c must be zero), than one such vector "c" is
c vector "a" plus the cross product ab x v.
c
c dx,dy,dz Output The x, y and z coordinates of the vertices of the
c regular polygon centered at point "a", with the first
c vertex at point "b", and in the plane containing
c points "a", "b" and "c". Size np.
c The vertices will be ordered around the center in the
c direction as points "a", "b" and "c".
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not three or more.
c 2 if points "a", "b" and "c" are not distinct,
c or are colinear.
c
c np Input The number of sides or vertices on the regular polygon.
c Must be at least 3.
c
c rc Output The radius of the circumscribed circle, or distance
c from point "a" to any vertex "d", including "b".
c rc = sqrt ((bx-ax)**2 + (by-ay)**2 + (bz-az)**2).
c rc = ri * sec (pi / np).
c
c ri Output The radius of the inscribed circle, or perpendicular
c distance from point "a" to any edge of the polygon.
c ri = rc * cos (pi / np).
c
c sl Output The length of an edge of the regular polygon.
c sl = 2.0 * rc * sin (pi / np).
c sl = 2.0 * ri * tan (pi / np).
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of a polygon vertex.
dimension dx (1)
c---- Coordinate y of a polygon vertex.
dimension dy (1)
c---- Coordinate z of a polygon vertex.
dimension dz (1)
c.... Local variables.
c---- Cosine of angle pi / np.
common /laptpoly/ costh
c---- Component x of normal to polygon.
common /laptpoly/ ex
c---- Component y of normal to polygon.
common /laptpoly/ ey
c---- Component z of normal to polygon.
common /laptpoly/ ez
c---- Index in local array.
common /laptpoly/ n
c---- Mathematical constant, pi.
common /laptpoly/ pi
c---- Coordinate x of a vertex.
common /laptpoly/ px
c---- Coordinate y of a vertex.
common /laptpoly/ py
c---- Coordinate z of a vertex.
common /laptpoly/ pz
c---- Rotation matrix, theta around "e".
common /laptpoly/ rotm (3,3)
c---- Sine of angle pi / np.
common /laptpoly/ sinth
c---- Angle 2.0 * pi / np.
common /laptpoly/ theta
c---- Length of normal vector "e".
common /laptpoly/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpoly. Find the vertices of a regular polygon.',
cbug & ' np=',i5,' tol=',1pe13.5 /
cbug & ' Center at a, first vertex at b, plane thru abc:' /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14)
cbug write ( 3, 9901) np, tol, ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... Initialize.
c---- Mathematical constant, pi.
pi = 3.14159265358979323
nerr = 0
c.... Test for input errors.
c---- No polygon.
if (np .lt. 3) then
nerr = 1
go to 210
endif
c.... Find the symmetry axis of the polygon.
call aptvpln (ax, ay, az, bx, by, bz, cx, cy, cz, 1, tol,
& ex, ey, ez, vlen, nerr)
c---- Plane undefined.
if (vlen .le. 0.0) then
nerr = 2
go to 210
endif
c.... Find the circumscribed and inscribed circle radii, the edge length,
c.... and the polygon area.
theta = 2.0 * pi / np
costh = cos (0.5 * theta)
sinth = sin (0.5 * theta)
rc = sqrt ((bx - ax)**2 + (by - ay)**2 + (bz - az)**2)
ri = rc * costh
sl = 2.0 * rc * sinth
ap = 0.5 * np * ri * sl
c.... Generate the vertices by rotation around axis "e" by angle theta.
call aptrota (theta, 1, ex, ey, ez, tol, rotm, nerr)
px = bx
py = by
pz = bz
c---- Loop over polygon vertices.
do 110 n = 1, np
dx(n) = px
dy(n) = py
dz(n) = pz
call aptmopv (rotm, 0, ax, ay, az, px, py, pz, 1, tol, nerr)
c---- End of loop over polygon vertices.
110 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptpoly results:' /
cbug & ' ri,rc= ',1p2e22.14 /
cbug & ' sl,ap= ',1p2e22.14 /
cbug & (i3,' dx,dy,dz=',1p3e22.14))
cbug 9903 format (/ ' closure:' /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' px,py,pz=',1p3e22.14)
cbug write ( 3, 9902) ri, rc, sl, ap,
cbug & (n, dx(n), dy(n), dz(n), n = 1, np)
cbug write ( 3, 9903) bx, by, bz, px, py, pz
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptpoly. (+1 line.)
end
UCRL-WEB-209832