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