subroutine aptil20 (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 APTIL20
c
c call aptil20 (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: aptil20 Updated 1989 December 27 15:50.
c aptil20 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 icosahedron, with 20 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 Glossary:
c
c ax,ay,az Output The x, y, z coordinates of vertex "a" of a triangle.
c Size 20 * (4**levmax / 3) (integer arithmatic).
c
c bx,by,bz Output The x, y, z coordinates of vertex "b" of a triangle.
c Size 20 * (4**levmax / 3) (integer arithmatic).
c
c cx,cy,cz Output The x, y, z coordinates of vertex "c" of a triangle.
c Size 20 * (4**levmax / 3) (integer arithmatic).
c
c ex1, ex2 Output The x coordinates of the beginning and end of edge n.
c Size 30 * 4**(levmax - 1).
c
c ey1, ey2 Output The y coordinates of the beginning and end of edge n.
c Size 30 * 4**(levmax - 1).
c
c ez1, ez2 Output The z coordinates of the beginning and end of edge n.
c Size 30 * 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 20 * (4**levmax / 3) (integer arithmatic).
c
c levmax Input Max level of tiling. Number of triangles
c is 20 * 4**(levmax - 1).
c
c nedges Output Number of edges = 30 * 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 (20)
c---- Coordinate y of vertex a of triangle.
dimension ay (20)
c---- Coordinate z of vertex a of triangle.
dimension az (20)
c---- Coordinate x of vertex b of triangle.
dimension bx (20)
c---- Coordinate y of vertex b of triangle.
dimension by (20)
c---- Coordinate z of vertex b of triangle.
dimension bz (20)
c---- Coordinate x of vertex c of triangle.
dimension cx (20)
c---- Coordinate y of vertex c of triangle.
dimension cy (20)
c---- Coordinate z of vertex c of triangle.
dimension cz (20)
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 central angle of icosa edge.
common /laptil20/ cosph
c---- Index of triangle being quartered.
common /laptil20/ n
c---- Index of triangle being formed.
common /laptil20/ nn
c---- Last index of penultimate level.
common /laptil20/ nmax
c---- Mathematical constant, pi.
common /laptil20/ pi
c---- Sine of central angle of edge.
common /laptil20/ sinph
c---- Radius of a vertex.
common /laptil20/ vrad
c---- Coordinate x of icosahedron vertex.
common /laptil20/ vx (12)
c---- Coordinate y of icosahedron vertex.
common /laptil20/ vy (12)
c---- Coordinate z of icosahedron vertex.
common /laptil20/ vz (12)
c.... Initialize.
c---- Mathematical constant, pi.
pi = 3.14159265358979323
nedges = 0
c.... Test for input errors.
nerr = 0
if (levmax .lt. 1) then
nerr = 1
go to 210
endif
c.... Generate the 12 vertices of the regular icosahedron, on a unit sphere.
cosph = 1.0 / sqrt (5.0)
sinph = 2.0 / sqrt (5.0)
vx( 1) = 0.0
vy( 1) = 0.0
vz( 1) = 1.0
vx( 2) = sinph
vy( 2) = vy(1)
vz( 2) = cosph
vx( 3) = sinph * cos (0.4 * pi)
vy( 3) = sinph * sin (0.4 * pi)
vz( 3) = vz(2)
vx( 4) = sinph * cos (0.8 * pi)
vy( 4) = sinph * sin (0.8 * pi)
vz( 4) = vz(2)
vx( 5) = vx(4)
vy( 5) = -vy(4)
vz( 5) = vz(2)
vx( 6) = vx(3)
vy( 6) = -vy(3)
vz( 6) = vz(2)
vx( 7) = -vx(4)
vy( 7) = vy(4)
vz( 7) = -vz(2)
vx( 8) = -vx(3)
vy( 8) = vy(3)
vz( 8) = -vz(2)
vx( 9) = -vx(2)
vy( 9) = vy(2)
vz( 9) = -vz(2)
vx(10) = -vx(3)
vy(10) = -vy(3)
vz(10) = -vz(2)
vx(11) = -vx(4)
vy(11) = -vy(4)
vz(11) = -vz(2)
vx(12) = vx(1)
vy(12) = vy(1)
vz(12) = -vz(1)
c.... Generate the first 20 equilateral triangles (level 1).
do 110 n = 1, 20
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)
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)
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)
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(6)
cy(4) = vy(6)
cz(4) = vz(6)
ax(5) = vx(1)
ay(5) = vy(1)
az(5) = vz(1)
bx(5) = vx(6)
by(5) = vy(6)
bz(5) = vz(6)
cx(5) = vx(2)
cy(5) = vy(2)
cz(5) = vz(2)
ax(6) = vx(2)
ay(6) = vy(2)
az(6) = vz(2)
bx(6) = vx(7)
by(6) = vy(7)
bz(6) = vz(7)
cx(6) = vx(3)
cy(6) = vy(3)
cz(6) = vz(3)
ax(7) = vx(3)
ay(7) = vy(3)
az(7) = vz(3)
bx(7) = vx(7)
by(7) = vy(7)
bz(7) = vz(7)
cx(7) = vx(8)
cy(7) = vy(8)
cz(7) = vz(8)
ax(8) = vx(3)
ay(8) = vy(3)
az(8) = vz(3)
bx(8) = vx(8)
by(8) = vy(8)
bz(8) = vz(8)
cx(8) = vx(4)
cy(8) = vy(4)
cz(8) = vz(4)
ax(9) = vx(4)
ay(9) = vy(4)
az(9) = vz(4)
bx(9) = vx(8)
by(9) = vy(8)
bz(9) = vz(8)
cx(9) = vx(9)
cy(9) = vy(9)
cz(9) = vz(9)
ax(10) = vx(4)
ay(10) = vy(4)
az(10) = vz(4)
bx(10) = vx(9)
by(10) = vy(9)
bz(10) = vz(9)
cx(10) = vx(5)
cy(10) = vy(5)
cz(10) = vz(5)
ax(11) = vx(5)
ay(11) = vy(5)
az(11) = vz(5)
bx(11) = vx(9)
by(11) = vy(9)
bz(11) = vz(9)
cx(11) = vx(10)
cy(11) = vy(10)
cz(11) = vz(10)
ax(12) = vx(5)
ay(12) = vy(5)
az(12) = vz(5)
bx(12) = vx(10)
by(12) = vy(10)
bz(12) = vz(10)
cx(12) = vx(6)
cy(12) = vy(6)
cz(12) = vz(6)
ax(13) = vx(6)
ay(13) = vy(6)
az(13) = vz(6)
bx(13) = vx(10)
by(13) = vy(10)
bz(13) = vz(10)
cx(13) = vx(11)
cy(13) = vy(11)
cz(13) = vz(11)
ax(14) = vx(6)
ay(14) = vy(6)
az(14) = vz(6)
bx(14) = vx(11)
by(14) = vy(11)
bz(14) = vz(11)
cx(14) = vx(2)
cy(14) = vy(2)
cz(14) = vz(2)
ax(15) = vx(2)
ay(15) = vy(2)
az(15) = vz(2)
bx(15) = vx(11)
by(15) = vy(11)
bz(15) = vz(11)
cx(15) = vx(7)
cy(15) = vy(7)
cz(15) = vz(7)
ax(16) = vx(7)
ay(16) = vy(7)
az(16) = vz(7)
bx(16) = vx(12)
by(16) = vy(12)
bz(16) = vz(12)
cx(16) = vx(8)
cy(16) = vy(8)
cz(16) = vz(8)
ax(17) = vx(8)
ay(17) = vy(8)
az(17) = vz(8)
bx(17) = vx(12)
by(17) = vy(12)
bz(17) = vz(12)
cx(17) = vx(9)
cy(17) = vy(9)
cz(17) = vz(9)
ax(18) = vx(9)
ay(18) = vy(9)
az(18) = vz(9)
bx(18) = vx(12)
by(18) = vy(12)
bz(18) = vz(12)
cx(18) = vx(10)
cy(18) = vy(10)
cz(18) = vz(10)
ax(19) = vx(10)
ay(19) = vy(10)
az(19) = vz(10)
bx(19) = vx(12)
by(19) = vy(12)
bz(19) = vz(12)
cx(19) = vx(11)
cy(19) = vy(11)
cz(19) = vz(11)
ax(20) = vx(11)
ay(20) = vy(11)
az(20) = vz(11)
bx(20) = vx(12)
by(20) = vy(12)
bz(20) = vz(12)
cx(20) = vx(7)
cy(20) = vy(7)
cz(20) = vz(7)
c.... Generate triangles in all remaining levels.
c++++ Last index of penultimate level.
nmax = 20 * (4**(levmax - 1) / 3)
c---- First index of final level.
ibeg = nmax + 1
c---- Last index of final level.
iend = 20 * (4**levmax / 3)
if (nmax .eq. 0) go to 210
nn = 20
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
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)
c---- End of loop over pre-levmax triangles.
120 continue
210 return
c.... End of subroutine aptil20. (+1 line.)
end
UCRL-WEB-209832