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