subroutine aptil12 (levmax, ibeg, iend, ax, ay, az, bx, by, bz,
& cx, cy, cz, lev,
& nedges, ex1, ey1, ez1, ex2, ey2, ez2, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTIL12
c
c call aptil12 (levmax, ibeg, iend, ax, ay, az, bx, by, bz,
c & cx, cy, cz, lev,
c & nedges, ex1, ey1, ez1, ex2, ey2, ez2, nerr)
c
c Version: aptil12 Updated 1989 December 27 15:50.
c aptil12 Originated 1989 November 2 14:10.
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,
c based on an initial regular dodecahedron, with 12 pentagons.
c Each pentagon is divided into 5 triangles.
c Each additional level has 4 times as many triangles.
c All triangle vertices are in the surface of the unit sphere.
c For level 2 tiling and higher, a table of unique edges is
c returned, for plotting. For level 1, plot the first two
c edges of each triangle, plus the third edge of triangle 20.
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, ex1, ey1, ez1, ex2, ey2, ez2, nerr.
c
c Calls: aptpent
c
c Glossary:
c
c ax,ay,az Output The x, y, z coordinates of vertex "a" of a triangle.
c Size 60 * (4**levmax / 3) (integer arithmatic).
c
c bx,by,bz Output The x, y, z coordinates of vertex "b" of a triangle.
c Size 60 * (4**levmax / 3) (integer arithmatic).
c
c cx,cy,cz Output The x, y, z coordinates of vertex "c" of a triangle.
c Size 60 * (4**levmax / 3) (integer arithmatic).
c
c ex1, ex2 Output The x coordinates of the beginning and end of edge n.
c Size 90 * 4**(levmax - 1).
c
c ey1, ey2 Output The y coordinates of the beginning and end of edge n.
c Size 90 * 4**(levmax - 1).
c
c ez1, ez2 Output The z coordinates of the beginning and end of edge n.
c Size 90 * 4**(levmax - 1).
c
c ibeg Output The first index of final level.
c
c iend Output Last index of final level.
c
c lev Output Level of tiling of triangle n.
c Size 60 * (4**levmax / 3) (integer arithmatic).
c
c levmax Input Max level of tiling. Number of triangles
c is 60 * 4**(levmax - 1).
c
c nedges Output Number of edges = 90 * 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 (1)
c---- Coordinate y of vertex a of triangle.
dimension ay (1)
c---- Coordinate z of vertex a of triangle.
dimension az (1)
c---- Coordinate x of vertex b of triangle.
dimension bx (1)
c---- Coordinate y of vertex b of triangle.
dimension by (1)
c---- Coordinate z of vertex b of triangle.
dimension bz (1)
c---- Coordinate x of vertex c of triangle.
dimension cx (1)
c---- Coordinate y of vertex c of triangle.
dimension cy (1)
c---- Coordinate z of vertex c of triangle.
dimension cz (1)
c---- Coordinate x of beginning of edge.
dimension ex1 (1)
c---- Coordinate x of end of edge.
dimension ex2 (1)
c---- Coordinate y of beginning of edge.
dimension ey1 (1)
c---- Coordinate y of end of edge.
dimension ey2 (1)
c---- Coordinate z of beginning of edge.
dimension ez1 (1)
c---- Coordinate z of end of edge.
dimension ez2 (1)
c---- Level of tiling of triangle.
dimension lev (1)
c.... Local variables.
c---- Cosine of angle from z axis.
common /laptil12/ cosph1
c---- Cosine of angle from z axis.
common /laptil12/ cosph2
c---- Index of triangle being quartered.
common /laptil12/ n
c---- Index of triangle being formed.
common /laptil12/ nn
c---- Last index of penultimate level.
common /laptil12/ nmax
c---- Mathematical constant, pi.
common /laptil12/ pi
c---- Sine of angle from z axis.
common /laptil12/ sinph1
c---- Sine of angle from z axis.
common /laptil12/ sinph2
c---- Radius of a vertex.
common /laptil12/ vrad
c---- Coordinate x of dodecahedron vertex.
common /laptil12/ vx (20)
c---- Coordinate y of dodecahedron vertex.
common /laptil12/ vy (20)
c---- Coordinate z of dodecahedron vertex.
common /laptil12/ vz (20)
c---- Mathematical constant, pi.
pi = 3.14159265358979323
c.... Test for input errors.
nerr = 0
if (levmax .lt. 1) then
nerr = 1
go to 210
endif
c.... Generate the 20 vertices of a regular dodecahedron, on a unit sphere.
c.... Orientation is symmetric with regular icosahedron.
cosph1 = 0.794654472292
sinph1 = 0.607061998206
vx( 1) = cos (0.2 * pi) * sinph1
vy( 1) = sin (0.2 * pi) * sinph1
vz( 1) = cosph1
vx( 2) = cos (0.6 * pi) * sinph1
vy( 2) = sin (0.6 * pi) * sinph1
vz( 2) = vz(1)
vx( 3) = -sinph1
vy( 3) = 0.0
vz( 3) = vz(1)
vx( 4) = vx(2)
vy( 4) = -vy(2)
vz( 4) = vz(1)
vx( 5) = vx(1)
vy( 5) = -vy(1)
vz( 5) = vz(1)
cosph2 = 0.187592474085
sinph2 = 0.982246946377
vx( 6) = cosph1
vy( 6) = vy(2)
vz( 6) = cosph2
vx( 7) = -0.5 * vx(3)
vy( 7) = sin (0.4 * pi) * sinph2
vz( 7) = -cosph2
vx( 8) = 0.5 * vx(3)
vy( 8) = vy(7)
vz( 8) = cosph2
vx( 9) = -vx(6)
vy( 9) = vy(2)
vz( 9) = -cosph2
vx(10) = -sinph2
vy(10) = vy(3)
vz(10) = cosph2
vx(11) = -vx(6)
vy(11) = -vy(2)
vz(11) = -cosph2
vx(12) = 0.5 * vx(3)
vy(12) = -vy(7)
vz(12) = cosph2
vx(13) = -0.5 * vx(3)
vy(13) = -vy(7)
vz(13) = -cosph2
vx(14) = vx(6)
vy(14) = -vy(2)
vz(14) = cosph2
vx(15) = -vx(10)
vy(15) = vy(3)
vz(15) = -cosph2
vx(16) = -vx(2)
vy(16) = vy(2)
vz(16) = -vz(1)
vx(17) = -vx(1)
vy(17) = vy(1)
vz(17) = -vz(1)
vx(18) = -vx(1)
vy(18) = -vy(1)
vz(18) = -vz(1)
vx(19) = -vx(2)
vy(19) = -vy(2)
vz(19) = -vz(1)
vx(20) = -vx(3)
vy(20) = vy(3)
vz(20) = -vz(1)
cbugc***DEBUG begins.
cbug 9901 format (// 'aptil12 finding dodecahedron vertices.' //
cbug & ' n x',15x,'y',15x,'z')
cbug 9902 format (i3,3f16.12)
cbug write ( 3, 9901)
cbug write ( 3, 9902) (n, vx(n), vy(n), vz(n), n = 1, 20)
cbugc***DEBUG ends.
c.... Subdivide each dedecahedron face (a pentagon) into 5 triangles.
do 110 n = 1, 60
lev(n) = 1
110 continue
n = 0
call aptpent (n, 1, 2, 3, 4, 5, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 1, 5, 14, 15, 6, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 1, 6, 7, 8, 2, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 2, 8, 9, 10, 3, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 3, 10, 11, 12, 4, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 4, 12, 13, 14, 5, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 6, 15, 20, 16, 7, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 7, 16, 17, 9, 8, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 9, 17, 18, 11, 10, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 11, 18, 19, 13, 12, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 13, 19, 20, 15, 14, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
call aptpent (n, 16, 20, 19, 18, 17, vx, vy, vz,
& ax, ay, az, bx, by, bz, cx, cy, cz)
c.... Generate triangles in all remaining levels.
c++++ Last index of penultimate level.
nmax = 60 * (4**(levmax - 1) / 3)
c---- First index of final level.
ibeg = nmax + 1
c---- Last index of final level.
iend = 60 * (4**levmax / 3)
if (nmax .eq. 0) go to 210
nn = 60
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
cbugc***DEBUG begins.
cbug if (lev(nn) .gt. lev(nn-1)) then
cbug nedges = 0
cbug endif
cbugc***DEBUG ends.
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
ex1(nedges) = ax(nn)
ey1(nedges) = ay(nn)
ez1(nedges) = az(nn)
ex2(nedges) = bx(nn)
ey2(nedges) = by(nn)
ez2(nedges) = bz(nn)
nedges = nedges + 1
ex1(nedges) = bx(nn)
ey1(nedges) = by(nn)
ez1(nedges) = bz(nn)
ex2(nedges) = cx(nn)
ey2(nedges) = cy(nn)
ez2(nedges) = cz(nn)
nedges = nedges + 1
ex1(nedges) = cx(nn)
ey1(nedges) = cy(nn)
ez1(nedges) = cz(nn)
ex2(nedges) = ax(nn)
ey2(nedges) = ay(nn)
ez2(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
ex1(nedges) = ax(nn)
ey1(nedges) = ay(nn)
ez1(nedges) = az(nn)
ex2(nedges) = bx(nn)
ey2(nedges) = by(nn)
ez2(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
ex1(nedges) = ax(nn)
ey1(nedges) = ay(nn)
ez1(nedges) = az(nn)
ex2(nedges) = bx(nn)
ey2(nedges) = by(nn)
ez2(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
ex1(nedges) = ax(nn)
ey1(nedges) = ay(nn)
ez1(nedges) = az(nn)
ex2(nedges) = bx(nn)
ey2(nedges) = by(nn)
ez2(nedges) = bz(nn)
120 continue
cbugc***DEBUG begins.
cbug 9903 format (// 'n=',i5,' nn=',i6,' ibeg=',i6,' iend=',i6)
cbug write ( 3, 9903) n, nn, ibeg, iend
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptil12. (+1 line.)
end
UCRL-WEB-209832