subroutine aptil08 (levmax, ibeg, iend, ax, ay, az, bx, by, bz,
& cx, cy, cz, lev,
& nedges, ex, ey, ez, fx, fy, fz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTIL08
c
c call aptil08 (levmax, ibeg, iend, ax, ay, az, bx, by, bz,
c & cx, cy, cz, lev,
c & nedges, ex, ey, ez, fx, fy, fz, nerr)
c
c Version: aptil08 Updated 1993 May 3 18:20.
c aptil08 Originated 1990 December 17 16:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To tile a unit sphere with planar triangles, based on an
c initial regular octahedron, with 8 faces, each an equilateral
c triangle, with all vertices on the major axes, and all edges
c in one of the major planes through the origin.
c Each additional level has 4 times as many triangles, with all
c triangle vertices in the surface of the unit sphere.
c Each level is formed by projecting the midpoints of each
c triangle edge to the surface of the sphere, and joining the
c projected points, to form four new triangles.
c A table of unique edges is returned, for the final level.
c Flag nerr indicates any input error, if not zero.
c
c Input: levmax
c
c Output: ibeg, iend, ax, ay, az, bx, by, bz, cx, cy, cz, lev,
c nedges, ex, ey, ez, fx, fy, fz, nerr.
c
c Glossary:
c
c ax,ay,az Output The x, y, z coordinates of vertex "a" of a triangle.
c Size 8 * (4**levmax / 3) (integer arithmatic).
c
c bx,by,bz Output The x, y, z coordinates of vertex "b" of a triangle.
c Size 8 * (4**levmax / 3) (integer arithmatic).
c
c cx,cy,cz Output The x, y, z coordinates of vertex "c" of a triangle.
c Size 8 * (4**levmax / 3) (integer arithmatic).
c
c ex, fx Output The x coordinates of the beginning and end of edge n.
c Size 12 * 4**(levmax - 1).
c
c ey, fy Output The y coordinates of the beginning and end of edge n.
c Size 12 * 4**(levmax - 1).
c
c ez, fz Output The z coordinates of the beginning and end of edge n.
c Size 12 * 4**(levmax - 1).
c
c ibeg Output The index of the first triangle of the final level.
c ibeg = 1 + 8 * (4**(levmax - 1)/ 3) (integer
c arithmatic).
c
c iend Output The index of the last triangle of the final level.
c iend = 8 * (4**levmax / 3) (integer arithmatic).
c
c lev Output Level of tiling of triangle n, n = 1, iend.
c Size 8 * (4**levmax / 3) (integer arithmatic).
c
c levmax Input Maximum level of tiling. Number of triangles
c is 8 * 4**(levmax - 1).
c
c nedges Output The number of edges for the final tiling.
c nedges = 12 * 4**(levmax - 1).
c
c nerr Output Indicates an input error, if not 0.
c 1 if levmax is less than 1.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of vertex a of triangle.
dimension ax (8)
c---- Coordinate y of vertex a of triangle.
dimension ay (8)
c---- Coordinate z of vertex a of triangle.
dimension az (8)
c---- Coordinate x of vertex b of triangle.
dimension bx (8)
c---- Coordinate y of vertex b of triangle.
dimension by (8)
c---- Coordinate z of vertex b of triangle.
dimension bz (8)
c---- Coordinate x of vertex c of triangle.
dimension cx (8)
c---- Coordinate y of vertex c of triangle.
dimension cy (8)
c---- Coordinate z of vertex c of triangle.
dimension cz (8)
c---- Coordinate x of beginning of edge.
dimension ex (8)
c---- Coordinate x of end of edge.
dimension fx (8)
c---- Coordinate y of beginning of edge.
dimension ey (8)
c---- Coordinate y of end of edge.
dimension fy (8)
c---- Coordinate z of beginning of edge.
dimension ez (8)
c---- Coordinate z of end of edge.
dimension fz (8)
c---- Level of tiling of triangle.
dimension lev (1)
c.... Local variables.
c---- Index of triangle being quartered.
common /laptil08/ n
c---- Index of triangle being formed.
common /laptil08/ nn
c---- Last index of penultimate level.
common /laptil08/ nmax
c---- Radius of a vertex.
common /laptil08/ vrad
c---- Coordinate x of octahedron vertex.
dimension vx (6)
c---- Coordinate y of octahedron vertex.
dimension vy (6)
c---- Coordinate z of octahedron vertex.
dimension vz (6)
data (vx(n), n = 1,6) / 0.0, 1.0, 0.0, -1.0, 0.0, 0.0 /
data (vy(n), n = 1,6) / 0.0, 0.0, 1.0, 0.0, -1.0, 0.0 /
data (vz(n), n = 1,6) / 1.0, 0.0, 0.0, 0.0, 0.0, -1.0 /
cbugc***DEBUG begins.
cbug 8801 format (/ 'aptil08 covering a unit sphere with triangles.' /
cbug & ' levmax=',i3)
cbugc***DEBUG ends.
c.... Initialize.
nedges = 0
c.... Test for input errors.
nerr = 0
if (levmax .lt. 1) then
nerr = 1
go to 210
endif
c.... Generate the first 8 equilateral triangles (level 1).
c.... The eight faces have vertices 1-2-3, 1-3-4, 1-4-5, 1-5-2,
c.... 6-2-3, 6-3-4, 6-4-5 and 6-5-2.
do 110 n = 1, 8
lev(n) = 1
110 continue
ax(1) = vx(1)
ay(1) = vy(1)
az(1) = vz(1)
bx(1) = vx(2)
by(1) = vy(2)
bz(1) = vz(2)
cx(1) = vx(3)
cy(1) = vy(3)
cz(1) = vz(3)
nedges = nedges + 1
ex(nedges) = ax(1)
ey(nedges) = ay(1)
ez(nedges) = az(1)
fx(nedges) = bx(1)
fy(nedges) = by(1)
fz(nedges) = bz(1)
nedges = nedges + 1
ex(nedges) = bx(1)
ey(nedges) = by(1)
ez(nedges) = bz(1)
fx(nedges) = cx(1)
fy(nedges) = cy(1)
fz(nedges) = cz(1)
ax(2) = vx(1)
ay(2) = vy(1)
az(2) = vz(1)
bx(2) = vx(3)
by(2) = vy(3)
bz(2) = vz(3)
cx(2) = vx(4)
cy(2) = vy(4)
cz(2) = vz(4)
nedges = nedges + 1
ex(nedges) = ax(2)
ey(nedges) = ay(2)
ez(nedges) = az(2)
fx(nedges) = bx(2)
fy(nedges) = by(2)
fz(nedges) = bz(2)
nedges = nedges + 1
ex(nedges) = bx(2)
ey(nedges) = by(2)
ez(nedges) = bz(2)
fx(nedges) = cx(2)
fy(nedges) = cy(2)
fz(nedges) = cz(2)
ax(3) = vx(1)
ay(3) = vy(1)
az(3) = vz(1)
bx(3) = vx(4)
by(3) = vy(4)
bz(3) = vz(4)
cx(3) = vx(5)
cy(3) = vy(5)
cz(3) = vz(5)
nedges = nedges + 1
ex(nedges) = ax(3)
ey(nedges) = ay(3)
ez(nedges) = az(3)
fx(nedges) = bx(3)
fy(nedges) = by(3)
fz(nedges) = bz(3)
nedges = nedges + 1
ex(nedges) = bx(3)
ey(nedges) = by(3)
ez(nedges) = bz(3)
fx(nedges) = cx(3)
fy(nedges) = cy(3)
fz(nedges) = cz(3)
ax(4) = vx(1)
ay(4) = vy(1)
az(4) = vz(1)
bx(4) = vx(5)
by(4) = vy(5)
bz(4) = vz(5)
cx(4) = vx(2)
cy(4) = vy(2)
cz(4) = vz(2)
nedges = nedges + 1
ex(nedges) = ax(4)
ey(nedges) = ay(4)
ez(nedges) = az(4)
fx(nedges) = bx(4)
fy(nedges) = by(4)
fz(nedges) = bz(4)
nedges = nedges + 1
ex(nedges) = bx(4)
ey(nedges) = by(4)
ez(nedges) = bz(4)
fx(nedges) = cx(4)
fy(nedges) = cy(4)
fz(nedges) = cz(4)
ax(5) = vx(6)
ay(5) = vy(6)
az(5) = vz(6)
bx(5) = vx(2)
by(5) = vy(2)
bz(5) = vz(2)
cx(5) = vx(3)
cy(5) = vy(3)
cz(5) = vz(3)
nedges = nedges + 1
ex(nedges) = ax(5)
ey(nedges) = ay(5)
ez(nedges) = az(5)
fx(nedges) = bx(5)
fy(nedges) = by(5)
fz(nedges) = bz(5)
ax(6) = vx(6)
ay(6) = vy(6)
az(6) = vz(6)
bx(6) = vx(3)
by(6) = vy(3)
bz(6) = vz(3)
cx(6) = vx(4)
cy(6) = vy(4)
cz(6) = vz(4)
nedges = nedges + 1
ex(nedges) = ax(6)
ey(nedges) = ay(6)
ez(nedges) = az(6)
fx(nedges) = bx(6)
fy(nedges) = by(6)
fz(nedges) = bz(6)
ax(7) = vx(6)
ay(7) = vy(6)
az(7) = vz(6)
bx(7) = vx(4)
by(7) = vy(4)
bz(7) = vz(4)
cx(7) = vx(5)
cy(7) = vy(5)
cz(7) = vz(5)
nedges = nedges + 1
ex(nedges) = ax(7)
ey(nedges) = ay(7)
ez(nedges) = az(7)
fx(nedges) = bx(7)
fy(nedges) = by(7)
fz(nedges) = bz(7)
ax(8) = vx(6)
ay(8) = vy(6)
az(8) = vz(6)
bx(8) = vx(5)
by(8) = vy(5)
bz(8) = vz(5)
cx(8) = vx(2)
cy(8) = vy(2)
cz(8) = vz(2)
nedges = nedges + 1
ex(nedges) = ax(8)
ey(nedges) = ay(8)
ez(nedges) = az(8)
fx(nedges) = bx(8)
fy(nedges) = by(8)
fz(nedges) = bz(8)
c.... Generate triangles in all remaining levels.
c++++ Last index of penultimate level.
nmax = 8 * (4**(levmax - 1) / 3)
c---- First index of final level.
ibeg = nmax + 1
c---- Last index of final level.
iend = 8 * (4**levmax / 3)
if (nmax .eq. 0) go to 210
nn = 8
c---- Loop over pre-levmax triangles.
do 120 n = 1, nmax
c.... Move the midpoints of the triangle edges to the sphere surface,
c.... and join to make a new triangle.
nn = nn + 1
lev(nn) = lev(n) + 1
c.... Restart edge array for each new level of triangles.
if (lev(nn) .gt. lev(nn-1)) then
nedges = 0
endif
ax(nn) = 0.5 * (ax(n) + bx(n))
ay(nn) = 0.5 * (ay(n) + by(n))
az(nn) = 0.5 * (az(n) + bz(n))
vrad = sqrt (ax(nn)**2 + ay(nn)**2 + az(nn)**2)
ax(nn) = ax(nn) / vrad
ay(nn) = ay(nn) / vrad
az(nn) = az(nn) / vrad
bx(nn) = 0.5 * (bx(n) + cx(n))
by(nn) = 0.5 * (by(n) + cy(n))
bz(nn) = 0.5 * (bz(n) + cz(n))
vrad = sqrt (bx(nn)**2 + by(nn)**2 + bz(nn)**2)
bx(nn) = bx(nn) / vrad
by(nn) = by(nn) / vrad
bz(nn) = bz(nn) / vrad
cx(nn) = 0.5 * (ax(n) + cx(n))
cy(nn) = 0.5 * (ay(n) + cy(n))
cz(nn) = 0.5 * (az(n) + cz(n))
vrad = sqrt (cx(nn)**2 + cy(nn)**2 + cz(nn)**2)
cx(nn) = cx(nn) / vrad
cy(nn) = cy(nn) / vrad
cz(nn) = cz(nn) / vrad
nedges = nedges + 1
ex(nedges) = ax(nn)
ey(nedges) = ay(nn)
ez(nedges) = az(nn)
fx(nedges) = bx(nn)
fy(nedges) = by(nn)
fz(nedges) = bz(nn)
nedges = nedges + 1
ex(nedges) = bx(nn)
ey(nedges) = by(nn)
ez(nedges) = bz(nn)
fx(nedges) = cx(nn)
fy(nedges) = cy(nn)
fz(nedges) = cz(nn)
nedges = nedges + 1
ex(nedges) = cx(nn)
ey(nedges) = cy(nn)
ez(nedges) = cz(nn)
fx(nedges) = ax(nn)
fy(nedges) = ay(nn)
fz(nedges) = az(nn)
c.... Connect initial vertex a to the adjacent new vertices.
nn = nn + 1
lev(nn) = lev(n) + 1
ax(nn) = ax(n)
ay(nn) = ay(n)
az(nn) = az(n)
bx(nn) = ax(nn-1)
by(nn) = ay(nn-1)
bz(nn) = az(nn-1)
cx(nn) = cx(nn-1)
cy(nn) = cy(nn-1)
cz(nn) = cz(nn-1)
nedges = nedges + 1
ex(nedges) = ax(nn)
ey(nedges) = ay(nn)
ez(nedges) = az(nn)
fx(nedges) = bx(nn)
fy(nedges) = by(nn)
fz(nedges) = bz(nn)
c.... Connect initial vertex b to the adjacent new vertices.
nn = nn + 1
lev(nn) = lev(n) + 1
ax(nn) = bx(n)
ay(nn) = by(n)
az(nn) = bz(n)
bx(nn) = bx(nn-2)
by(nn) = by(nn-2)
bz(nn) = bz(nn-2)
cx(nn) = ax(nn-2)
cy(nn) = ay(nn-2)
cz(nn) = az(nn-2)
nedges = nedges + 1
ex(nedges) = ax(nn)
ey(nedges) = ay(nn)
ez(nedges) = az(nn)
fx(nedges) = bx(nn)
fy(nedges) = by(nn)
fz(nedges) = bz(nn)
c.... Connect initial vertex c to the adjacent new vertices.
nn = nn + 1
lev(nn) = lev(n) + 1
ax(nn) = cx(n)
ay(nn) = cy(n)
az(nn) = cz(n)
bx(nn) = cx(nn-3)
by(nn) = cy(nn-3)
bz(nn) = cz(nn-3)
cx(nn) = bx(nn-3)
cy(nn) = by(nn-3)
cz(nn) = bz(nn-3)
nedges = nedges + 1
ex(nedges) = ax(nn)
ey(nedges) = ay(nn)
ez(nedges) = az(nn)
fx(nedges) = bx(nn)
fy(nedges) = by(nn)
fz(nedges) = bz(nn)
c---- End of loop over pre-levmax triangles.
120 continue
210 continue
cbugc***DEBUG begins.
cbug 8802 format (/ 'aptil08 results: nerr =',i3,' iend=',i5,
cbug & ' nedges=',i5)
cbug write ( 3, 8802) nerr, iend, nedges
cbug if (nerr .ne. 0) return
cbug 8803 format (' n,l,a,b,c=',2i3,9f7.4)
cbug 8804 format (' n,e,f =', i3,6f7.4)
cbug write ( 3, 8803) (n, lev(n), ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(n),
cbug & n = 1, iend)
cbug write ( 3, 8804) (n, ex(n), ey(n), ez(n), fx(n), fy(n), fz(n),
cbug & n = 1, nedges)
cbugc***DEBUG ends.
return
c.... End of subroutine aptil08. (+1 line.)
end
UCRL-WEB-209832