function node_xyz = sphere_gridpoints_icos2 ( factor, node_num ) %*****************************************************************************80 % %% SPHERE_GRIDPOINTS_ICOS2 returns icosahedral grid points on a sphere. % % Discussion: % % With FACTOR = 1, the grid has 20 triangular faces and 12 nodes. % % With FACTOR = 2, each triangle of the icosahedron is subdivided into % 2x2 subtriangles, resulting in 80 faces and % 42 = 12 + 20 * 3 * (1)/2 + 20 * 0 ) nodes. % % With FACTOR = 3, each triangle of the icosahedron is subdivided into % 3x3 subtriangles, resulting in 180 faces and % 92 ( = 12 + 20 * 3 * (2)/2 + 20 * 1 ) nodes. % % In general, each triangle is subdivided into FACTOR*FACTOR subtriangles, % resulting in 20 * FACTOR * FACTOR faces and % 12 % + 20 * 3 * (FACTOR-1) / 2 % + 20 * (FACTOR-2) * (FACTOR-1) / 2 nodes. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 28 August 2010 % % Author: % % John Burkardt % % Parameters: % % Input, integer FACTOR, the subdivision factor, which must % be at least 1. % % Input, integer NODE_NUM, the number of nodes, as reported % by SPHERE_IMP_GRID_ICOS_NUM. % % Output, real NODE_XYZ(3,NODE_NUM), the node coordinates. % % Local Parameters: % % POINT_NUM, EDGE_NUM, FACE_NUM and FACE_ORDER_MAX are counters % associated with the icosahedron, and POINT_COORD, EDGE_POINT, % FACE_ORDER and FACE_POINT are data associated with the icosahedron. % We need to refer to this data to generate the grid. % % NODE counts the number of nodes we have generated so far. At the % end of the routine, it should be equal to NODE_NUM. % % % Size the icosahedron. % [ point_num, edge_num, face_num, face_order_max ] = icos_num ( ); % % Set the icosahedron. % [ point_coord, edge_point, face_order, face_point ] = icos_shape ( ... point_num, edge_num, face_num, face_order_max ); % % Generate the point coordinates. % % A. Points that are the icosahedral vertices. % node = 0; node_xyz(1:3,1:point_num) = point_coord(1:3,1:point_num); % % B. Points in the icosahedral edges, at % 1/FACTOR, 2/FACTOR, ..., (FACTOR-1)/FACTOR. % node = 12; for edge = 1 : edge_num a = edge_point(1,edge); b = edge_point(2,edge); % % Determine the "distance" = angle between points A and B. % theta = sphere_distance_xyz ( point_coord(1:3,a), point_coord(1:3,b) ); % % Polarize B into BP + BN and normalize BN. % [ bn, bp ] = r8vec_polarize ( 3, point_coord(1:3,b), point_coord(1:3,a) ); bn(1:3,1) = bn(1:3,1) / norm ( bn ); % % March from A to B, by taking equally spaced angles from 0 to THETA. % F = 0 => ANGLE = 0 => A % F = FACTOR => ANGLE = THETA => B % for f = 1 : factor - 1 node = node + 1; angle = ( f * theta ) / factor; node_xyz(1:3,node) = cos ( angle ) * point_coord(1:3,a) ... + sin ( angle ) * bn(1:3,1); end end % % C. Points in the icosahedral faces. % for face = 1 : face_num a = face_point(1,face); b = face_point(2,face); c = face_point(3,face); % % Determine the "distance" = angle between points A and B, A and C. % theta_ab = sphere_distance_xyz ( point_coord(1:3,a), point_coord(1:3,b) ); theta_ac = sphere_distance_xyz ( point_coord(1:3,a), point_coord(1:3,c) ); % % Polarize B = BP + BN and normalize BN, C = CP + CN, and normalize CN. % [ bn, bp ] = r8vec_polarize ( 3, point_coord(1:3,b), point_coord(1:3,a) ); bn(1:3,1) = bn(1:3,1) / norm ( bn ); [ cn, cp ] = r8vec_polarize ( 3, point_coord(1:3,c), point_coord(1:3,a) ); cn(1:3,1) = cn(1:3,1) / norm ( cn ); % % March AB from A to B: % FA = 0 => ANGLE = 0 => AB = A % FA = FACTOR => ANGLE = THETA_AB => AB = B % % March AC from A to C: % FA = 0 => ANGLE = 0 => AC = A % FA = FACTOR => ANGLE = THETA_AC => AC = C % for fa = 2 : factor - 1 % % Determine points AB and AC that use cos ( FA / FACTOR ) of A % and cos ( ( FACTOR - FA ) / FACTOR ) of B or C. % angle = ( fa * theta_ab ) / factor; ab(1:3,1) = cos ( angle ) * point_coord(1:3,a) + sin ( angle ) * bn(1:3,1); angle = ( fa * theta_ac ) / factor; ac(1:3,1) = cos ( angle ) * point_coord(1:3,a) + sin ( angle ) * cn(1:3,1); % % Determine the "distance" = angle between points AB and AC. % theta_bc = sphere_distance_xyz ( ab, ac ); % % Polarize AC into ACP + ACN and normalize ACN. % [ acn, acp ] = r8vec_polarize ( 3, ac, ab ); acn(1:3,1) = acn(1:3,1) / norm ( acn ); % % The interval between AB and AC is broken into FA intervals. % Go from 1 to FA - 1. % for fbc = 1 : fa - 1 node = node + 1; angle = fbc * theta_bc / fa; node_xyz(1:3,node) = cos ( angle ) * ab(1:3,1) ... + sin ( angle ) * acn(1:3,1); end end end return end