Home > voicebox > upolyhedron.m

upolyhedron

PURPOSE ^

UPOLYHEDRON calculate uniform polyhedron characteristics

SYNOPSIS ^

function [vlist,edgeq,flist,info]=upolyhedron(w,md)

DESCRIPTION ^

 UPOLYHEDRON calculate uniform polyhedron characteristics

 Inputs:  W    Specifies the desired polyhedron in one of three forms:
                 (a) name e.g. W='cube'; precede by 'dual' for dual or
                 'laevo' or 'dextro' [default] for chiral polyhedra
                 (b) index in the list given below e.g. W=6 is the cube; negative for dual
                     n.{1,2,3,4,5} gives n-sided prism, antiprism, grammic prism, grammic antiprism or grammic crossed antiprism
                 (c) Wythoff symbol (see below) e.g. W=[3 0 2 4] is the cube
                          |p q r, p|q r, p q|r or p q r| respectively with 0 for |
                         using -1 instead of 0 gives the dual of the polyhedron
                         using -2 (or -3 for dual) gives a reflected version of a snub polyhedron
          MD   specifies a mode string
                 'g' plot image
                 'w' plot wireframe
                 'f' plot coloured faces
                 'v' number vertices
                 't' plot vertex figure (a slice centred on a vertex)

                 ? homogeneous output coordinates
                 ? create net
                 ? segment faces to remove internal portions of the surfaces
                 ? size: [max] diameter=1, [longest] edge=1, include anisotropic normalization
                   to minimize variance of vertex radius or edge vector length
                 ? orientation: vertex at the top, largest stable base at bottom

 Outputs:
          VLIST(:,7)  gives the [x y z d n e t] for each vertex
                       x,y,z = position, d=distance from origin, n=valency, e=edge index, t=type (-ve for reflected)
          EDGEQ(:,9) has one row for each direction of each edge:
                         1  v1   first vertex (normally start)
                         2  v2   second vertex
                         3  f1   first face (normally on left)
                         4  f2   second face
                         5  ev1  next edge around vertex 1 (normally anticlockwise)
                         6  ef1  next edge around f1 (normally anticlockwise)
                         7  er   reverse edge
                         8  z    twisted edge: clockwise neighbours around v1 and v2 are on the same face
                         9  sf   swap face order: ???
                        10  sv   swap vertex order: v2 preceeds v1 around f1
          FLIST(:,7)  gives the [x y z d n e t] for each face
                       x,y,z = unit normal, d=distance from origin, n=valency, e=edge index, t=type (-ve for reflected)
          INFO        structure containing the following fields:
                         hemi      true if faces are hemispherical (i.e. pass through the origin)
                         onesided  true is one-sided (like a moebius strip)
                         snub      true if a snub polyhedron

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [vlist,edgeq,flist,info]=upolyhedron(w,md)
0002 % UPOLYHEDRON calculate uniform polyhedron characteristics
0003 %
0004 % Inputs:  W    Specifies the desired polyhedron in one of three forms:
0005 %                 (a) name e.g. W='cube'; precede by 'dual' for dual or
0006 %                 'laevo' or 'dextro' [default] for chiral polyhedra
0007 %                 (b) index in the list given below e.g. W=6 is the cube; negative for dual
0008 %                     n.{1,2,3,4,5} gives n-sided prism, antiprism, grammic prism, grammic antiprism or grammic crossed antiprism
0009 %                 (c) Wythoff symbol (see below) e.g. W=[3 0 2 4] is the cube
0010 %                          |p q r, p|q r, p q|r or p q r| respectively with 0 for |
0011 %                         using -1 instead of 0 gives the dual of the polyhedron
0012 %                         using -2 (or -3 for dual) gives a reflected version of a snub polyhedron
0013 %          MD   specifies a mode string
0014 %                 'g' plot image
0015 %                 'w' plot wireframe
0016 %                 'f' plot coloured faces
0017 %                 'v' number vertices
0018 %                 't' plot vertex figure (a slice centred on a vertex)
0019 %
0020 %                 ? homogeneous output coordinates
0021 %                 ? create net
0022 %                 ? segment faces to remove internal portions of the surfaces
0023 %                 ? size: [max] diameter=1, [longest] edge=1, include anisotropic normalization
0024 %                   to minimize variance of vertex radius or edge vector length
0025 %                 ? orientation: vertex at the top, largest stable base at bottom
0026 %
0027 % Outputs:
0028 %          VLIST(:,7)  gives the [x y z d n e t] for each vertex
0029 %                       x,y,z = position, d=distance from origin, n=valency, e=edge index, t=type (-ve for reflected)
0030 %          EDGEQ(:,9) has one row for each direction of each edge:
0031 %                         1  v1   first vertex (normally start)
0032 %                         2  v2   second vertex
0033 %                         3  f1   first face (normally on left)
0034 %                         4  f2   second face
0035 %                         5  ev1  next edge around vertex 1 (normally anticlockwise)
0036 %                         6  ef1  next edge around f1 (normally anticlockwise)
0037 %                         7  er   reverse edge
0038 %                         8  z    twisted edge: clockwise neighbours around v1 and v2 are on the same face
0039 %                         9  sf   swap face order: ???
0040 %                        10  sv   swap vertex order: v2 preceeds v1 around f1
0041 %          FLIST(:,7)  gives the [x y z d n e t] for each face
0042 %                       x,y,z = unit normal, d=distance from origin, n=valency, e=edge index, t=type (-ve for reflected)
0043 %          INFO        structure containing the following fields:
0044 %                         hemi      true if faces are hemispherical (i.e. pass through the origin)
0045 %                         onesided  true is one-sided (like a moebius strip)
0046 %                         snub      true if a snub polyhedron
0047 
0048 % This software is based closely on a Mathmatica program described in [1] which, in turn, was based on a
0049 % C program described in [2].
0050 %
0051 % Wythoff Symbol
0052 %    p,q,r define a spherical triangle whose angles at the corners are pi/p, pi/q and pi/r;
0053 %    this triangle tiles the sphere if repeatedly reflected in its sides. The
0054 %    polyhedron vertices are at the reflections of a seed vertex as follows
0055 %    where the Vertex configuration gives the polygon orders of the faces around each vertex:
0056 %       |p q r  : Vertex at a point such that when rotated around any of p,q,r by twice the angle at that
0057 %                 corner is displaced by the same distance for each corner. This is a snub polyhedron and
0058 %                 only even numbers of reflections are used to generate vertices. Configuration {3 p 3 q 3 r}
0059 %        p|q r  : Vertex at p. Configuration = {q r q r ... q r} with 2n terms where n is the numerator of p
0060 %        p q|r  : Vertex on pq and on the bisector of angle r. Configuration {p 2r q 2r}
0061 %        p q r| : Vertex at incentre: meeting point of all the angle bisectors. Configuration {2p 2q 2r}
0062 %        |3/2 5/3 3 5/2  This special case is the great dirhombicosidodecahedron. It is a bit weird because
0063 %                 many of the edges are shared by four faces instead of the usual two.
0064 %    If two of p,q,r = 2 then the third is arbitrary (prisms and antiprisms), otherwise only the numerators
0065 %    1:5 can occur, and 4 and 5 cannot occur together. If all are integers then the poyhedron is convex.
0066 %
0067 % References:
0068 %   [1] R. E. Maeder. Uniform polyhedra. The Mathematica Journal, 3 (4): 48–57, 1993.
0069 %   [2] Z. Har’El. Uniform solution for uniform polyhedra. Geometriae Dedicata, 47: 57–110, 1993.
0070 %   [3] H. S. M. Coxeter, M. S. Longuet-Higgins, and J. C. P. Miller. Uniform polyhedra.
0071 %       Philosophical Transactions of the Royal Society A, 246 (916): 401–450, May 1954.
0072 %   [4] P. W. Messer. Closed-form expressions for uniform polyhedra and their duals.
0073 %       Discrete and Computational Geometry, 27 (3): 353–375, Jan. 2002.
0074 
0075 %%%% BUGS and SUGGESTIONS %%%%%%
0076 % (1) we should ensure the "first" edges of the vertices and faces are consistent
0077 % (2) need to sort faces and vertices into a type order
0078 % (3) should ensure that for a non-chiral polyhedron, the vertex polarity alternates
0079 % (4) w=75 does not work
0080 % (5) dual of henispherical poyhedron
0081 % (6) flist not calculated correctly for duals
0082 % (7) add additional stuff into info.*
0083 % (8) vertex figures seem to include additional (duplicated) lines
0084 % (9) sort out when reflected vertices are really rotationally congruent
0085 % (10) could optionally colour by face type
0086 % (11) order alphabetically by noun
0087 % (12) allow abbreviated names + prisms with preceding decimal number
0088 % (13) calculate correct names
0089 % (14) include names for duals
0090 % (15) make "names" etc persistent
0091 
0092 % Example slide show:
0093 % for i=1:74, disp(num2str(i)); upolyhedron(i); pause(2); end
0094 % for i=5:10, disp(num2str(i)); for j=1:5, upolyhedron(i+j/10); pause(2); end, end
0095 
0096 %      Copyright (C) Mike Brookes 1997
0097 %      Version: $Id: upolyhedron.m 713 2011-10-16 14:45:43Z dmb $
0098 %
0099 %   VOICEBOX is a MATLAB toolbox for speech processing.
0100 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0101 %
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103 %   This program is free software; you can redistribute it and/or modify
0104 %   it under the terms of the GNU General Public License as published by
0105 %   the Free Software Foundation; either version 2 of the License, or
0106 %   (at your option) any later version.
0107 %
0108 %   This program is distributed in the hope that it will be useful,
0109 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0110 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0111 %   GNU General Public License for more details.
0112 %
0113 %   You can obtain a copy of the GNU General Public License from
0114 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0115 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0117 
0118 % [line nn]  referes to the kaleido Mathematica program on which this is based
0119 % see http://www.mathconsult.ch/showroom/unipoly/unipoly.html#Images
0120 
0121 % Variables used
0122 %   adjacent    adjacency matrix
0123 %   chi         characteristic
0124 %   cosa        cos of the angle subtended by a half edge
0125 %   d           density
0126 %   e           number of edges
0127 %   even        number of even faces to remove
0128 %   f           number of faces
0129 %   fi(n)       number of faces of type i
0130 %   g           order of group
0131 %   gamma(n)    included spherical triangle angle between centre of face,
0132 %               vertex and edge for face type i
0133 %   hemiQ       =1 for hemispherical faces (includes polyhedron centre)
0134 %   incidence   faces incident at vertex
0135 %   k           type of group (2=dihedral, 3=tetrahedral, 4=octahedral, 5=icosahedral)
0136 %   m           number of edges meeting at each vertex (valence)
0137 %   mi(n)       Number of faces of type i
0138 %   n           number of face types
0139 %   ni(n)       number of edges on i'th face type
0140 %   p           1st Wythoff number
0141 %   q           2nd Wythoff number
0142 %   r           3rd Wythoff number
0143 %   rot(m)      vertex configuration: anti-clockwise list of face types
0144 %               around a vertex
0145 %   snub(m)     =1 for snub triangles
0146 %   snubQ       =1 for snub polyhedron
0147 %   v           number of vertices
0148 %   vcoord      vertex coordinates
0149 %   wy(4)       Wythoff Symbol
0150 %   wyd(4)      Wythoff Symbol Denominators
0151 %   wyn(4)      Wythoff Symbol Numerators
0152 
0153 % names: [name abbreviation synonym dual dual-synonym]
0154 persistent names wys prefs
0155 if ~numel(names)
0156     prefs={'dual '; 'dextro '; 'laevo '};
0157     names={
0158         'Tetrahedron' 'Tet' '' '' '';                              %  1
0159         'truncated tetrahedron' 'Tut' '' '' '';                    %  2
0160         'octahemioctahedron' 'Oho' '' '' '';                       %  3
0161         'tetrahemihexahedron' 'Thah' '' '' '';                      %  4
0162         'Octahedron' 'Oct' '' 'Cube' '';                               %  5
0163         'cube' 'Cube' '' 'Octahedron' '';                                     %  6
0164         'cuboctahedron' 'Co' '' 'rhombic dodecahedron' '';                            %  7
0165         'truncated octahedron' 'Toe' '' 'tetrakishexahedron' '';                     %  8
0166         'truncated cube' 'Tic' '' 'triakisoctahedron' '';                           %  9
0167         'rhombicuboctahedron' 'Sirco' '' 'deltoidal icositetrahedron' '';                      % 10
0168         'truncated cuboctahedron' '' '' 'disdyakisdodecahedron' '';                  % 11
0169         'snub cube' 'Snic' '' 'pentagonal icositetrahedron' '';                                % 12
0170         'small cubicuboctahedron' 'Socco' '' 'small hexacronic icositetrahedron' '';                  % 13
0171         'great cubicuboctahedron' 'Gocco' '' 'great hexacronic icositetrahedron' '';                  % 14
0172         'cubohemioctahedron' 'Cho' '' 'hexahemioctacron' '';                       % 15
0173         'cubitruncated cuboctahedron' 'Cotco' '' 'tetradyakishexahedron' '';              % 16
0174         'great rhombicuboctahedron' 'Girco' '' 'great deltoidal icositetrahedron' '';                % 17 or Querco
0175         'small rhombihexahedron' 'Sroh' '' 'small rhombihexacron' '';                   % 18
0176         'stellated truncated hexahedron' 'Quith' '' 'great triakisoctahedron' '';              % 19
0177         'great truncated cuboctahedron' 'Quitco' '' 'great disdyakisdodecahedron' '';               % 20
0178         'great rhombihexahedron' 'Groh' '' 'great rhombihexacron' '';                   % 21
0179         'icosahedron' 'Ike' '' 'dodecahedron' '';                              % 22
0180         'dodecahedron' 'Doe' '' 'icosahedron' '';                             % 23
0181         'icosidodecahedron' 'Id' '' 'rhombic triacontahedron' '';                        % 24
0182         'truncated icosahedron' 'Ti' '' 'pentakisdodecahedron' '';                    % 25
0183         'truncated dodecahedron' 'Tid' '' 'triakisicosahedron' '';                   % 26
0184         'rhombicosidodecahedron' 'Srid' '' 'deltoidal hexecontahedron' '';                   % 27
0185         'truncated icosidodechedon' 'Grid' 'great rhombiicosidodecahedron' 'disdyakistriacontahedron' '';                % 28
0186         'snub dodecahedron' 'Snid' '' 'pentagonal hexecontahedron' '';                        % 29
0187         'small ditrigonal icosidodecahedron' 'Sidtid' '' 'small triambic icosahedron' '';         % 30
0188         'small icosicosidodecahedron' 'Siid' '' 'small icosacronic hexecontahedron' '';              % 31
0189         'small snub icosicosidodecahedron' 'Seside' '' 'small hexagonal hexecontahedron' '';            % 32
0190         'small dodecicosidodecahedron' 'Saddid' '' 'small dodecacronic hexecontahedron' '';                % 33
0191         'small stellated dodecahedron' 'Sissid' '' 'great dodecahedron' '';                % 34
0192         'great dodecahedron' 'Gad' '' 'small stellated dodecahedron' '';                       % 35
0193         'dodecadodecahedron' 'Did' '' 'medial rhombic triacontahedron' '';                       % 36
0194         'truncated great dodecahedron' 'Tigid' '' 'small stellapentakisdodecahedron' '';                % 37
0195         'rhombidodecadodecahedron' 'Raded' '' 'medial deltoidal hexecontahedron' '';                 % 38
0196         'small rhombidodecahedron' 'Sird' '' 'small rhombidodecacron' '';                 % 39
0197         'snub dodecadodecahedron' 'Siddid' '' 'medial pentagonal hexecontahedron' '';                  % 40
0198         'ditrigonal dodecadodecahedron' 'Ditdid' '' 'medial triambic icosahedron' '';              % 41
0199         'great ditrigonal dodecicosidodecahedron' 'Gidditdid' '' 'great ditrigonal dodecacronic hexecontahedron' '';    % 42
0200         'small ditrigonal dodecicosidodecahedron' 'Sidditdid' '' 'small ditrigonal dodecacronic hexecontahedron' '';    % 43
0201         'icosidodecadodecahedron' 'Ided' '' 'medial icosacronic hexecontahedron' '';                  % 44
0202         'icositruncated dodecadodecahedron' 'Idtid' '' 'tridyakisicosahedron' '';        % 45
0203         'snub icosidodecadodecahedron' 'Sided' '' 'medial hexagonal hexecontahedron' '';             % 46
0204         'great ditrigonal icosidodecahedron' 'Gidtid' '' 'great triambic icosahedron' '';       % 47
0205         'great icosicosidodecahedron' 'Giid' '' 'great icosacronic hexecontahedron' '';              % 48
0206         'small icosihemidodecahedron' 'Seihid' '' 'small icosihemidodecacron' '';              % 49
0207         'small dodecicosahedron' 'Siddy' '' 'small dodecicosacron' '';                   % 50
0208         'small dodecahemidodecahedron' 'Sidhid' '' 'small dodecahemidodecacron' '';             % 51
0209         'great stellated dodecahedron' 'Gissid' '' 'great icosahedron' '';               % 52
0210         'great icosahedron' 'Gike' '' 'great stellated dodecahedron' '';                        % 53
0211         'great icosidodecahedron' 'Gid' '' 'great rhombic triacontahedron' '';                  % 54
0212         'great truncated icosahedron' 'Tiggy' '' 'great stellapentakisdodecahedron' '';              % 55
0213         'rhombicosahedron' 'Ri' '' 'rhombicosacron' '';                         % 56
0214         'great snub icosidodecahedron' 'Gosid' '' 'great pentagonal hexecontahedron' '';             % 57
0215         'small stellated truncated dodecahedron' 'Quitsissid' '' 'great pentakisdodekahedron' '';     % 58
0216         'truncated dodecadodecahedron' 'Quitdid' '' 'medial disdyakistriacontahedron' '';               % 59
0217         'inverted snub dodecadodecahedron' 'Isdid' '' 'medial inverted pentagonal hexecontahedron' '';           % 60
0218         'great dodecicosidodecahedron' 'Gaddid' '' 'great dodecacronic hexecontahedron' '';               % 61
0219         'small dodecahemicosahedron' 'Sidhei' '' 'small dodecahemicosacron' '';               % 62
0220         'great dodecicosahedron' 'Giddy' '' 'great dodecicosacron' '';                   % 63
0221         'great snub dodecicosidodecahedron' 'Gisdid' '' 'great hexagonal hexecontahedron' '';          % 64
0222         'great dodecahemicosahedron' 'Gidhei' '' 'great dodecahemicosacron' '';               % 65
0223         'great stellated truncated dodecahedron' 'Quitgissid' '' 'great triakisicosahedron' '';      % 66
0224         'great rhombicosidodecahedron' 'Qrid' '' 'great deltoidal hexecontahedron' '';                % 67 or Qrid
0225         'great truncated icosidodecahedron' 'Gaquatid' '' 'great disdyakistriacontahedron' '';          % 68
0226         'great inverted snub icosidodecahedron' 'Gisid' '' 'great inverted pentagonal hexecontahedron' '';       % 69
0227         'great dodecahemidodecahedron' 'Gidhid' '' 'great dodecahemidodecacron' '';                % 70
0228         'great icosihemidodecahedron' 'Geihid' '' 'great icosihemidodecacron' '';              % 71
0229         'small retrosnub icosicosidodecahedron' 'Sirsid' '' 'small hexagrammic hexecontahedron' '';       % 72
0230         'great rhombidodecahedron' 'Gird' '' 'great rhombidodecacron' '';                 % 73
0231         'great retrosnub icosidodecahedron' 'Girsid' '' 'great pentagrammic hexecontahedron' '';        % 74
0232         'great dirhombicosidodecahedron' 'Gidrid' '' 'great dirhombicosidodecacron' '';             % 75
0233         'pentagonal prism' 'Pip' '' '' '';                         % 76
0234         'pentagonal antiprism' 'Pap' '' '' '';                     % 77
0235         'pentagrammic prism' 'Stip' '' '' '';                       % 78
0236         'pentagrammic antiprism' 'Stap' '' '' '';                   % 79
0237         'pentagrammic crossed antiprism'  'Starp'   '' '' '';          % 80
0238         'triangular prism' '' '' '' '';      % 81
0239         'triangular antiprism' '' '' '' '';     % 82
0240         'square prism' '' '' '' '';     % 83
0241         'square antiprism' '' '' '' ''};    % 84
0242 
0243     % could multiply by 12 to make integer
0244     wys=[
0245         2, 3, 2, 3;        %  1   'tetrahedron' 'Tet';
0246         3, 2, 3, 3;        %  2   'truncated tetrahedron' 'Tut';
0247         3, 3/2, 3, 3;      %  3   'octahemioctahedron' 'Oho';
0248         3, 3/2, 3, 2;      %  4   'tetrahemihexahedron' 'Thah';
0249         2, 4, 2, 3;        %  5   'octahedron' 'Oct';
0250         2, 3, 2, 4;        %  6   'cube' 'Cube';
0251         2, 2, 3, 4;        %  7   'cuboctahedron' 'Co';
0252         3, 2, 4, 3;        %  8   'truncated octahedron' 'Toe';
0253         3, 2, 3, 4;        %  9   'truncated cube' 'Tic';
0254         3, 3, 4, 2;        % 10   'rhombicuboctahedron' 'Sirco';
0255         4, 2, 3, 4;        % 11   'truncated cuboctahedron' '?';
0256         1, 2, 3, 4;        % 12   'snub cube' 'Snic';
0257         3, 3/2, 4, 4;      % 13   'small cubicuboctahedron' 'Socco';
0258         3, 3, 4, 4/3;      % 14   'great cubicuboctahedron' 'Gocco';
0259         3, 4/3, 4, 3;      % 15   'cubohemioctahedron' 'Cho';
0260         4, 4/3, 3, 4;      % 16   'cubitruncated cuboctahedron' 'Cotco';
0261         3, 3/2, 4, 2;      % 17   'great rhombicuboctahedron' 'Girco?';
0262         4, 3/2, 2, 4;      % 18   'small rhombihexahedron' 'Sroh';
0263         3, 2, 3, 4/3;      % 19   'stellated truncated hexahedron' 'Quith';
0264         4, 4/3, 2, 3;      % 20   'great truncated cuboctahedron' 'Quitco';
0265         4, 4/3, 3/2, 2;    % 21   'great rhombihexahedron' 'Groh';
0266         2, 5, 2, 3;        % 22   'icosahedron' 'Ike';
0267         2, 3, 2, 5;        % 23   'dodecahedron' 'Doe';
0268         2, 2, 3, 5;        % 24   'icosidodecahedron' 'Id';
0269         3, 2, 5, 3;        % 25   'truncated icosahedron' 'Ti';
0270         3, 2, 3, 5;        % 26   'truncated dodecahedron' 'Tid';
0271         3, 3, 5, 2;        % 27   'rhombicosidodecahedron' 'Srid';
0272         4, 2, 3, 5;        % 28   'truncated icosidodecahedron' 'Grid'; or 'great rhombiicosidodecahedron' Grid
0273         1, 2, 3, 5;        % 29   'snub dodecahedron' 'Snid';
0274         2, 3, 5/2, 3;      % 30   'small ditrigonal icosidodecahedron' 'Sidtid';
0275         3, 5/2, 3, 3;      % 31   'small icosicosidodecahedron' 'Siid';
0276         1, 5/2, 3, 3;      % 32   'small snub icosicosidodecahedron' 'Seside';
0277         3, 3/2, 5, 5;      % 33   'small dodecicosidodecahedron' 'Saddid';
0278         2, 5, 2, 5/2;      % 34   'small stellated dodecahedron' 'Sissid';
0279         2, 5/2, 2, 5;      % 35   'great dodecahedron' 'Gad';
0280         2, 2, 5/2, 5;      % 36   'dodecadodecahedron' 'Did';
0281         3, 2, 5/2, 5;      % 37   'truncated great dodecahedron' 'Tigid';
0282         3, 5/2, 5, 2;      % 38   'rhombidodecadodecahedron' 'Raded';
0283         4, 2, 5/2, 5;      % 39   'small rhombidodecahedron' 'Sird';
0284         1, 2, 5/2, 5;      % 40   'snub dodecadodecahedron' 'Siddid';
0285         2, 3, 5/3, 5;      % 41   'ditrigonal dodecadodecahedron' 'Ditdid';
0286         3, 3, 5, 5/3;      % 42   'great ditrigonal dodecicosidodecahedron' 'Gidditdid'
0287         3, 5/3, 3, 5;      % 43   'small ditrigonal dodecicosidodecahedron' 'Sidditdid'
0288         3, 5/3, 5, 3;      % 44   'icosidodecadodecahedron' 'Ided';
0289         4, 5/3, 3, 5;      % 45   'icositruncated dodecadodecahedron' 'Idtid';
0290         1, 5/3, 3, 5;      % 46   'snub icosidodecadodecahedron' 'Sided';
0291         2, 3/2, 3, 5;      % 47   'great ditrigonal icosidodecahedron' 'Gidtid';
0292         3, 3/2, 5, 3;      % 48   'great icosicosidodecahedron' 'Giid';
0293         3, 3/2, 3, 5;      % 49   'small icosihemidodecahedron' 'Seihid';
0294         4, 3/2, 3, 5;      % 50   'small dodecicosahedron' 'Siddy';
0295         3, 5/4, 5, 5;      % 51   'small dodecahemidodecahedron' 'Sidhid';
0296         2, 3, 2, 5/2;      % 52   'great stellated dodecahedron' 'Gissid';
0297         2, 5/2, 2, 3;      % 53   'great icosahedron' 'Gike';
0298         2, 2, 5/2, 3;      % 54   'great icosidodecahedron' 'Gid';
0299         3, 2, 5/2, 3;      % 55   'great truncated icosahedron' 'Tiggy';
0300         4, 2, 5/2, 3;      % 56   'rhombicosahedron' 'Ri';
0301         1, 2, 5/2, 3;      % 57   'great snub icosidodecahedron' 'Gosid';
0302         3, 2, 5, 5/3;      % 58   'small stellated truncated dodecahedron' 'Quitsissid'
0303         4, 5/3, 2, 5;      % 59   'truncated dodecadodecahedron' 'Quitdid';
0304         1, 5/3, 2, 5;      % 60   'inverted snub dodecadodecahedron' 'Isdid';
0305         3, 5/2, 3, 5/3;    % 61   'great dodecicosidodecahedron' 'Gaddid';
0306         3, 5/3, 5/2, 3;    % 62   'small dodecahemicosahedron' 'Sidhei';
0307         4, 5/3, 5/2, 3;    % 63   'great dodecicosahedron' 'Giddy';
0308         1, 5/3, 5/2, 3;    % 64   'great snub dodecicosidodecahedron' 'Gisdid';
0309         3, 5/4, 5, 3;      % 65   'great dodecahemicosahedron' 'Gidhei';
0310         3, 2, 3, 5/3;      % 66   'great stellated truncated dodecahedron' 'Quitgissid'
0311         3, 5/3, 3, 2;      % 67   'quasirhombicosidodecahedron' 'Qrid'; 'great rhombicosidodecahedron' 'Nonconvex great rhombicosidodecahedron
0312         4, 5/3, 2, 3;      % 68   'great truncated icosidodecahedron' 'Gaquatid';
0313         1, 5/3, 2, 3;      % 69   'great inverted snub icosidodecahedron' 'Gisid';
0314         3, 5/3, 5/2, 5/3;  % 70   'great dodecahemidodecahedron' 'Gidhid';
0315         3, 3/2, 3, 5/3;    % 71   'great icosihemidodecahedron' 'Geihid';
0316         1, 3/2, 3/2, 5/2;  % 72   'small retrosnub icosicosidodecahedron' 'Sirsid';
0317         4, 3/2, 5/3, 2;    % 73   'great rhombidodecahedron' 'Gird';
0318         1, 3/2, 5/3, 2;    % 74   'great retrosnub icosidodecahedron' 'Girsid';
0319         5, 3/2, 5/3, 3;    % 75   'great dirhombicosidodecahedron' 'Gidrid';
0320         3, 2, 5, 2;        % 76   'pentagonal prism' 'Pip';
0321         1, 2, 2, 5;        % 77   'pentagonal antiprism' 'Pap';
0322         3, 2, 5/2, 2;      % 78   'pentagrammic prism' 'Stip';
0323         1, 2, 2, 5/2;      % 79   'pentagrammic antiprism' 'Stap';
0324         1, 2, 2, 5/3;     % 80   'pentagrammic crossed antiprism'  'Starp'
0325         3, 2, 3, 2;        % 81   'triangular prism' ;
0326         1, 2, 2, 3;        % 82   'triangular antiprism' ;
0327         3, 2, 4, 2;        % 83   'square prism' ;
0328         1, 2, 2, 4];        % 84   'square antiprism' ;
0329 end
0330 
0331 if nargin<2
0332     md='';
0333 end
0334 dual=0;  % dual polyhedron
0335 dextro=0; % reflected (for snub polyhedra only)
0336 if ischar(w)
0337     i=0;
0338     while i<length(prefs)
0339         i=i+1;
0340         if length(w)>length(prefs{i}) && strcmpi(w(1:length(prefs{i})),prefs{i})
0341             switch i
0342                 case 1
0343                     dual=1;
0344                 case 2
0345                     dextro=1;
0346             end
0347             w=w(length(prefs{i})+1:end);
0348             i=0;
0349         end
0350     end
0351 
0352     wyidx=find(strcmpi(w,names(:)),1);
0353     if isempty(wyidx)
0354         % check for prism names here
0355         error('Cannot find %s',w);
0356     end
0357     if wyidx>3*size(names,1)
0358         dual=1-dual;
0359     end
0360     wyidx=1+rem(wyidx-1,size(names,1));
0361     wy=wys(wyidx,:);
0362 else
0363     if length(w)==1
0364         dual=w<0;
0365         wyidx=abs(w);
0366         wyjdx=floor(wyidx);
0367         wykdx=round(10*(wyidx-wyjdx));
0368         switch wykdx
0369             case 1  % prism
0370                 wy=[3 2 wyjdx 2];
0371             case 2  % antiprism
0372                 wy=[1 2 2 wyjdx];
0373             case 3  % grammic prism
0374                 wy=[3 2 wyjdx/2 2];
0375             case 4  % grammic antiprism
0376                 wy=[1 2 2 wyjdx/2];
0377             case 5  % grammic crossed antiprism
0378                 wy=[1 2 2 wyjdx/3];
0379             otherwise
0380                 if wyjdx>=1 && wyjdx<=size(wys,1)
0381                     wy=wys(wyjdx,:);
0382                 else
0383                     error('Polyhedron number out of range');
0384                 end
0385         end
0386     elseif length(w)==4
0387         vbar=find(w<=0,1);
0388         if ~numel(vbar)
0389             error('Invalid Wythoff symbol: %g %g %g %g %g',w);
0390         end
0391         dual=mod(w(vbar),2);    % least significant bit indicates dual
0392         dextro=mod(1+w(vbar),4)>=2; % second bit indicates reflected version
0393         wy=[vbar w(1:vbar-1) w(vbar+1:4)];
0394         [xx,wyidx]=min(sum((wys-repmat(wy,size(wys,1),1)).^2,2));
0395         if abs(xx)>1e-3
0396             if sum(wy==2)<2
0397                 error('Invalid Wythoff symbol: %g %g %g %g %g',w);
0398             else
0399                 % need to sort out prism names here
0400             end
0401         end
0402     elseif length(w)==5
0403         if w(1)<=0 && all(round(12*w(2:end))==[18 40 36 30])   % [3/2 5/3 3 5/2]
0404             dual=w(1)<0;
0405             wyidx=75;
0406             wy=wys(wyidx,:);
0407         else
0408             error('Invalid Wythoff symbol: %g %g %g %g %g',w);
0409         end
0410     else
0411         error('Invalid polyhedron specification');
0412     end
0413 end
0414 
0415 
0416 [wyn,wyd]=rat(wy(2:4));    % convert to rational numbers
0417 wy(2:4)=wyn./wyd;           % force exact rational values
0418 p=wy(2);
0419 q=wy(3);
0420 r=wy(4);
0421 hemiQ = 0;          % includes hemispherical faces that go through polyhedron centre [ p q | r] with pq=p+q
0422 onesidedQ = 0;      % one-sided polyhedron (aka: non-orientable) [ p q r |] with exactly one of p,q,r having an even denominator
0423 evenQ = 0;           % identifies which of p, q, r has an even denominator
0424 even=1;
0425 snubQ = 0;          % snub polyhedron: [ | p q r ]
0426 allrot = 1;         % all vertices a congruent with rotations (no reflections needed)
0427 
0428 % call to AnalyseWythoff[line 105]
0429 
0430 k= max(wyn);
0431 if sum(wy(2:4)==2)>=2      % check if it is a prism
0432     % this is a prism - treat specially [note mathematica line 86 has redundant chack for 2 >5]
0433     g=4*k;      % order of group
0434     k=2;
0435 else                       % not a prism: only numerators 1,2,3,4,5 are allowed
0436     if (any(wy(2:4)<=1) || (k>5) || any(wyn==4) && any(wyn==5))
0437         error('Invalid Wythoff numbers [%d/%d %d/%d %d/%d]',[wyn;wyd]);
0438     end
0439     g=24*k/(6-k);       % order of group
0440 end
0441 if abs(wy(1))==5 % special case
0442     chi = 0;
0443     d = 0;
0444 else
0445     chi = (sum(wyn.^(-1))-1)*g/2;
0446     d = (sum(wy(2:4).^(-1))-1)*g/4;
0447     if d<0
0448         error('density < 0');
0449     end
0450 end
0451 if abs(wy(1))==5
0452     error('great dirhombicosidodecahedron not implemented');
0453 end
0454 switch abs(wy(1))
0455     case {1,5}    % [ | p q r ] snub polyhedron
0456         n=4;            % number of face types
0457         m=6;            % valence of vertices
0458         v=g/2;          % number of vertices
0459 %         ni=[3 wy(2:4)]; % type of each face
0460 %         mi=[3 1 1 1];   % number of faces of each type
0461         nimi=[3 wy(2:4); 3 1 1 1]';
0462 %         rot = [1 2+dextro 1 3-dextro 1 4];
0463         snubQ=1;
0464 %         snub=[1 0 1 0 1 0];
0465         rotsnub=[1 2+dextro 1 3-dextro 1 4; repmat([1 0],1,3)]';
0466     case 2  % [ p | q r ]
0467         n=2;
0468         m=2*wyn(1);
0469         v=g/m;
0470 %         ni=wy(3:4);
0471 %         mi=[p p];         % these face counts may be fractional: m counts their numerators
0472 % we could ill in the numerator denominator columns, nimi(:,3:4) here instead of later
0473         nimi=[wy(3:4); p p]';
0474 %         rot = repmat(1:2,1,wyn(1));
0475         rotsnub = repmat([1 0; 2 0],wyn(1),1);
0476     case 3  % [ p q | r ]
0477         n=3;
0478         m=4;
0479         v=g/2;
0480 %         ni=[2*r wy(2:3)];
0481 %         mi = [2 1 1];
0482         nimi=[2*r wy(2:3); 2 1 1]';
0483 %         rot = [1 2 1 3];
0484         rotsnub = [1 2 1 3; zeros(1,4)]';
0485         if abs(p-q/(q-1))<1e-6
0486             hemiQ=1;
0487             d=0;
0488             if (p~=2 && ~(wy(4)==3 && any(wy(2:3)==3)))
0489                 onesidedQ=1;
0490                 v=v/2;
0491                 chi=chi/2;
0492             end
0493         end
0494     case 4  % [ p q r | ]
0495         n=3;
0496         m=3;
0497         v=g;
0498 %         ni=2*wy(2:4);
0499 %         mi=ones(1,3);
0500         nimi=[2*wy(2:4); ones(1,3)]';
0501 %         rot=1:3;
0502         rotsnub=[1:3; zeros(1,3)]';
0503         allrot=(wy(2)==wy(3)) || (wy(2)==wy(4)) || (wy(3)==wy(4)); % isosceles so all
0504         evenden=find(~rem(wyd,2)); % check for even denominators
0505         if (length(evenden)>1)
0506             error('Multiple even denominators are not allowed');
0507         end
0508         if ~isempty(evenden)
0509             even = evenden;
0510             evenQ=1;
0511             onesidedQ = 1;
0512             d=0;
0513             v=v/2;
0514             chi = chi - g/wyn(even)/2;  % check this *****
0515         end
0516     otherwise
0517         error('Invalid polyhedron type %d',wy(1));
0518 end
0519 % [line 155] call sortAndMerge [sortAndMerge: line 220]
0520 % save stuff to prevent overwriting
0521 n_1=n;
0522 m_1=m;
0523 % now do sort merge based on rotsnub and nimi
0524 % rotsnub(:,2) = [rot snub]
0525 % nimi(:,4) = [ni mi mi-num mi-den]
0526 
0527 nimi(nimi(:,1)==2,1)=-2;            % make equal to -2 to force to bottom of list
0528 [nimi,inim]=sortrows(nimi,-1);
0529 jnim=zeros(n,1);
0530 jnim(inim)=1:n;
0531 msk=[~0; nimi(2:end,1)<nimi(1:end-1,1)]; % identify distinct values
0532 cmsk=cumsum(msk); % destination of each value
0533 nimi(msk,2)=sparse(1,cmsk,nimi(:,2));  % add up repreated counts
0534 nimi=nimi(msk,:);       % select just the distinct values
0535 % cmsk(jnim) maps original to new positions
0536 [nimi(:,3),nimi(:,4)]=rat(nimi(:,2));
0537 rotsnub(:,1)=cmsk(jnim(rotsnub(:,1)));
0538 even=cmsk(jnim(even));
0539 if nimi(end,1)<0        % remove digons
0540     if size(nimi,1)==1
0541         error('Degenerate polyhedron (digons only)');
0542     end
0543     m=m-nimi(end,3);  % note that m=sum(nimi(:,3)) not sum(nimi(:,2)) as you might expect
0544     nimi(end,:)=[];
0545     n=size(nimi,1);
0546     rotsnub(rotsnub(:,1)>n,:)=[];         % abolish references to digons
0547 end
0548 
0549 % original sort and merge
0550 
0551 % [ss,tr]=sort(-ni);          % sort into descending order
0552 % itr=zeros(1,n);
0553 % itr(tr)=1:n;
0554 % ni=ni(tr);
0555 % mi=mi(tr);
0556 % rot = itr(rot);
0557 % if evenQ
0558 %     even = itr(even);
0559 % end
0560 % % now merge equal faces
0561 % i=1;
0562 % while i<n
0563 %     if ni(i)==ni(i+1)
0564 %         mi(i)=mi(i)+mi(i+1);
0565 %         mi(i+1)=[];
0566 %         ni(i+1)=[];
0567 %         n=n-1;
0568 %         even=even-(even>i);
0569 %         rot = rot - (rot>i);
0570 %     else
0571 %         i=i+1;
0572 %     end
0573 % end
0574 % [mii,mij]=rat(mi);    % find numerator
0575 % i=find(ni==2);  % find digons
0576 % if ~isempty(i)
0577 %     if n==1
0578 %         error('Degenerate (digons only)');
0579 %     end
0580 %     i=i(1);
0581 %     m=m-mii(i);      % reduce total valance
0582 %     ni(i)=[];
0583 %     mi(i)=[];
0584 %     mii(i)=[];
0585 %     even=even-(even>i);
0586 %     if snubQ
0587 %         snub(rot==i)=[];
0588 %     end
0589 %     rot(rot==i)=[];
0590 %     rot = rot - (rot>i);
0591 %     n=n-1;
0592 % end
0593 %
0594 % % check that new version is correct
0595 %
0596 %
0597 % if any(nimi(:,1:3)~=[ni' mi' mii'])
0598 %     error('nimi has errors');
0599 % end
0600 % if snubQ
0601 %     if any(rotsnub~=[rot' snub'])
0602 %         error('rotsnub has errors');
0603 %     end
0604 % else
0605 %     if any(rotsnub(:,1)~=rot')
0606 %         error('rotsnub has errors');
0607 %     end
0608 % end
0609 % if evenQ && even~=even_1
0610 %     error('even has errors');
0611 % end
0612 
0613 % [line 159] Solve fundamental equations [line 266]
0614 % doesn't converge well for a dodecahedron
0615 % we could avoid the iteration when n=1
0616 %
0617 % We want to solve the following set of equations for cosa and gammai:
0618 %  (a) sum(mi.*gammai)=pi
0619 %  (b) cos(alphai) = cosa*sin(gammai)
0620 %      where ni, mi, alphai and gammai are vectors of length m, the number of edges at a vertex
0621 %      ni gives the number of sides in each face type in decreasing order (possibly fractional)
0622 %      mi is the repretition count of each face type (possibly fractional)
0623 %      alphai=pi * ni.^(-1) is the angle subtended by a half edge at the centre of a face
0624 %      gammai (> pi/2 - alphai) is half the spherical angle subtended by a face at the vertex
0625 %      cosa (< 1) is the distance of the edge centre from the origin if the radius of the vertex is unity
0626 % we solve this iteratively using gammai(1) as the independent variable
0627 
0628 
0629 alphai = pi*nimi(:,1).^(-1); % half of external angle of polygon
0630 calphai = cos(alphai);
0631 ca1 = calphai(1);
0632 calphai(1)=[];
0633 % initial values
0634 gammai = pi/2 - alphai; % half of internal angle of polygon
0635 delta = pi - nimi(:,2)'*gammai;
0636 % we could initialize it to be exact for a single face type
0637 % g1n=gammai(1)*pi/(mi*gammai')
0638 % cosa = ca1/sin(g1n);
0639 % gammai = [g1n asin(calphai/cosa)];
0640 % delta = pi - mi*gammai';
0641 % should check here to see if delta is better than before and g1n is valid
0642 iter=0;
0643 while (abs(delta) > 1e-10)
0644     g1n = gammai(1) + delta*tan(gammai(1))/(nimi(:,2)'*tan(gammai));
0645     if (g1n<0 || g1n>pi)
0646         error ('Gamma out of range');
0647     end
0648     cosa = ca1/sin(g1n);
0649     gammai = [g1n; asin(calphai/cosa)];
0650     delta = pi - nimi(:,2)'*gammai;
0651     iter=iter+1;
0652 end
0653 gamma = gammai;
0654 cosa = ca1/sin(gammai(1));
0655 
0656 % [line 165] postprocess special cases
0657 
0658 if evenQ
0659     nimi(even,:)=[];
0660     gamma(even)=[];
0661     nh=n-1;
0662     n=2*nh;
0663     m=4;
0664 %     ni=[ni 1+fliplr(ni-1).^(-1)];
0665 nimi=repmat(nimi,2,1);
0666 nimi(nh+1:end,1)=1+(nimi(nh:-1:1,1)-1).^(-1);
0667     gamma=[gamma; -flipud(gamma)];
0668 %     mi=repmat(3-nh,1,n);
0669 %     mii=mi;
0670 nimi(:,2:3)=repmat(3-nh,n,2);
0671     rotsnub=[1 nh n 1+nh; zeros(1,4)]';
0672 end
0673 if wy(1)==5
0674     error('not yet implemented');
0675     % needs to use rotsnub and nimi
0676     n=5;
0677     m=8;
0678     hemiQ = 1;
0679     d=0;
0680     mi=[4 1 1 1 1];
0681     mii=mi;
0682     ni=[4 ni(1:3) ni(1)/(ni(1)-1)];
0683     gamma=[pi/2 gamma(1:3) -gamma(1)];
0684     rot=[rot+[0 1 0 1 0 1] 1 5]; % [1 4 1 3 1 2 1 5]
0685     snub=[snub 1 0];
0686 end
0687 
0688 
0689 
0690 % [line 197] count vertices and faces [count: line 286]
0691 
0692 e = m*v/2;          % number of edges = valency * vertices /2
0693 [nii,nij]=rat(nimi(:,1));  % find numerators
0694 fi = v*nimi(:,3)./nii;
0695 f = sum(fi);
0696 if (d>0) && (gamma(1)> pi/2) % [might be better to use cycles rather than radians for gamma]
0697     d = fi(1)-d;
0698 end
0699 if wy(1)==5
0700     chi = v - e + f; % Euler characteristic
0701 end
0702 
0703 % [line 199] generate vertex coordinates and adjacency matrix [vertices: line 297]
0704 % we currently scale it so the edge mid-points have unit radius
0705 
0706 adj=zeros(v,m);         % j'th vertex adjacent to vertex i anti-clockwise;
0707 frot = zeros(1,v);
0708 rev = zeros(v,1);
0709 vlist = zeros(v+1,7); % vertex information: [x y z d n e t] for each vertex
0710 %                       x,y,z = position, d=distance from origin, n=valency, e=edge index, t=type (-ve for reflected)
0711 vlist(:,5)=m;           % all vertices have the same valency
0712 v1 = [0 0 1]/cosa;        % first vertex
0713 frot(1) = 1;
0714 adj(1,1) = 2;
0715 v2 = [2*cosa*sqrt(1-cosa^2), 0, 2*cosa^2-1]/cosa;   % second vertex
0716 if snubQ
0717     frot(2) = m*(1-rotsnub(m,2))+rotsnub(m,2); % inefficient
0718     adj(2,1)=1;
0719 else                    % reflexible
0720     frot(2) = 1;
0721     rev(2) = 1;
0722     adj(2,m)=1;
0723 end
0724 vlist(1,1:3)=v1;
0725 vlist(2,1:3)=v2;
0726 nv = 2;
0727 i = 1;
0728 skew=zeros(3,3);
0729 skewp=[6 7 2];  % index of positive entries
0730 skewn=[8 3 4];  % index of negative entries
0731 cosg=cos(2*gamma);
0732 sing=sin(2*gamma);
0733 eye3=eye(3);
0734 % veqth=cos(acos(cosa)/(d+1));   % threshold for vertex equality test = cosa * |vlist(1,1:3)|^2
0735 veqth=0.999999/(cosa^2);  % threshold for vertex equality test = cosa * |vlist(1,1:3)|^2
0736 while i<=nv  % loop for each vertex
0737     if rev(i)
0738         one = -1;
0739         start = m-1;
0740         limit = 1;
0741     else
0742         one = 1;
0743         start = 2;
0744         limit = m;
0745     end
0746     k = frot(i);    % rotation to use first
0747     v1 = vlist(i,1:3);    % the centre of rotation
0748     v1 = v1/sqrt(v1*v1');   % normalize to unit length [not clear why we don't make unit length in the first place]
0749     v2 = vlist(adj(i,start-one),1:3);
0750     for j=start:one:limit
0751         %             R=wwT+cos(x)(I-wwT)+sin(x)SKEW(w)
0752         %             SKEW(a) =  [0 -a3 a2; a3 0 -a1; -a2 a1 0]
0753         skew(skewp)=v1;
0754         skew(skewn)=-v1;
0755         rsym=v1'*v1;
0756         rotk=rotsnub(k,1);
0757         rotmat=rsym+cosg(rotk)*(eye3-rsym)-one*sing(rotk)*skew;
0758         v2=v2*rotmat;  % rotate v2 poition around v1
0759         vlist(nv+1,1:3)=v2;   % add into list in case it is good
0760         nvi=find(vlist(:,1:3)*v2'>veqth,1); % take the first matching vertex
0761         %         cmatch=1-vlist(nvi,:)*v2'*cosa^2      % diagnostic printout for vertex match test (0 for perfect match)
0762         adj(i,j)=nvi; % save as next vertex
0763         lastk = k;
0764         k=1+mod(k,m); % increment k circularly in the range 1:m
0765         if nvi>nv       % we have a new vertex
0766             if snubQ % if a snub polyhedron
0767                 % Mathematica: frot[[nvi]] = If[ !sn[[lastk]], lastk, If[ !sn[[k]], next[k, m], k ] ];
0768                 % A snub triangle plays the role of the next available snub triangle for the next vertex
0769                 frot(nvi)=lastk+rotsnub(lastk,2)*(rotsnub(k,2)*k-lastk+(1-rotsnub(k,2))*(1+mod(k,m)));
0770                 rev(nvi)=0;   % snub polyhedra always have rev=0
0771                 adj(nvi,1)=i;
0772             else
0773                 frot(nvi)=k;
0774                 rev(nvi)=(1+one)/2; % = 1-rev(i)
0775                 adj(nvi,1+(m-1)*rev(nvi))=i;
0776             end
0777             nv=nvi;
0778             if nv>v
0779                 error('Too many vertices found');
0780             end
0781         end
0782     end
0783     i=i+1;
0784 end
0785 if (nv~=v)
0786     error('Not all vertices found');
0787 end
0788 vlist(v+1,:)=[];    % remove the dummy extra vertex
0789 vlist(:,7)=1-2*rev;     % vertex types all all +1 or -1
0790 
0791 % construct edge map
0792 % edgeq: 1=v1 2=v2 3=f1 4=f2 5=ev1 6=ef1 7=er 8=z 9=sf 10=sv]
0793 %    1  v1   first vertex (normally start)
0794 %    2  v2   second vertex
0795 %    3  f1   first face (normally on left)
0796 %    4  f2   second face
0797 %    5  ev1  next edge around vertex 1 (normally anticlockwise)
0798 %    6  ef1  next edge around f1 (normally anticlockwise)
0799 %    7  er   reverse edge
0800 %    8  z    twisted edge: clockwise neighbours around v1 and v2 are on the same face
0801 %    9  sf   swap face order: ???
0802 %   10  sv   swap vertex order: v2 preceeds v1 around f1
0803 
0804 edgeq=zeros(v*m,10); % OLD *** ===[reverse_edge start_vertex left_face next_edge_at_vertex next_edge_on_face]
0805 edgeq(:,1:2)=[repmat((1:v)',m,1) adj(:)];  % 2=v_start 3=v_end
0806 edgeq(:,3:4)=edgeq(:,1:2)*[v 1; 1 v];      % encode start and end as a single integer
0807 [xx,ia]=sort(edgeq(:,3));
0808 [xx,ib]=sort(edgeq(:,4));
0809 edgeq(:,3:4)=0;
0810 ib(ib)=1:v*m;               % make reverse index
0811 edgeq(:,7)=ia(ib);          % index to reverse direction edge: 1=e_reverse
0812 edgeq(:,5)=[v+1:v*m 1:v];   % adjacent edge anti-clockwise around vertex: 6=e_next_at_start
0813 ia(edgeq(:,5))=1:v*m;       % adjacent edge clockwise around vertex; alternatively ia = mod((1:v*m)-v,v*m)'
0814 edgeq(:,6)=ia(edgeq(:,7));  % adjacent edge anti-clockwise around face [but not necessarily correct
0815 fpt=zeros(f,2);             % edge count and pointer to an edge on a face
0816 iff=0;                      % face index
0817 if onesidedQ
0818     edgeq(:,8)=reshape(mod(repmat(frot',1,m)+repmat((1:m),v,1),2),v*m,1); %face type modulo 2
0819     edgeq(:,8)=abs(edgeq(:,8)-edgeq(edgeq(:,6),8)); % 1 = twisted edge
0820     while iff<f
0821         iee=find(edgeq(:,3)==0,1); % find an edge without a left face
0822         if ~numel(iee)
0823             error('Not enough faces found');
0824         end
0825         iff=iff+1;
0826         fpt(iff,2)=iee;
0827         nee=0;
0828         flip=abs((edgeq(iee,6)==edgeq(edgeq(iee,7),5))-edgeq(iee,8)); % check if this initial edge is flipped on output
0829         while edgeq(iee,3)~=iff
0830             if edgeq(iee,3)  % use the reverse edge if this one is already taken
0831                 iee=edgeq(iee,7);
0832                 edgeq(pee,6)=iee;           % correct the previous exit edge taken
0833                 edgeq(iee,10)=1;     % use it in the reverse direction
0834                 edgeq(iee,6)=ia(iee); % normal exit for this direction
0835                 jee=edgeq(iee,5);  % alternate exit edge if flipping
0836             else
0837                 if flip && ~edgeq(iee,4) % we are entering on the usual reverse exit
0838                     edgeq(edgeq(iee,7),6)=edgeq(iee,5); % So give the reverse exit a valid path
0839                 end
0840                 jee=edgeq(edgeq(iee,7),5); % alternate exit edge if flipping
0841             end
0842             if edgeq(iee,3)
0843                 error('Edge already in use');
0844             end
0845             edgeq(iee,3)=iff;
0846             edgeq(edgeq(iee,7),4)=iff;      % mark right face in reverse edge
0847             nee=nee+1;
0848             flip=abs(flip-edgeq(iee,8));
0849             if flip   % use alternate edge
0850                 edgeq(iee,6)=jee;           % mark actual edge used
0851             end
0852             pee=iee;            % save previous iee value
0853             iee=edgeq(iee,6);   % step around the face
0854         end
0855         fpt(iff,1)=nee;
0856     end
0857 else  % a normal two-sided polygon
0858     while iff<f
0859         iee=find(edgeq(:,3)==0,1); % find an edge without a left face
0860         iff=iff+1;
0861         fpt(iff,2)=iee;
0862         nee=0;
0863         while edgeq(iee,3)~=iff
0864             if edgeq(iee,3)
0865                 error('Edge already in use');
0866             end
0867             edgeq(iee,3)=iff;
0868             edgeq(edgeq(iee,7),4)=iff;      % mark right face in reverse edge
0869             nee=nee+1;
0870             iee=edgeq(iee,6);   % step around the face
0871         end
0872         fpt(iff,1)=nee;
0873     end
0874 end
0875 % [(1:v*m)' edgeq]
0876 
0877 vlist(:,4)=sqrt(sum(vlist(:,1:3).^2,2));        % fill in the radii
0878 vlist(:,6)=(1:v)';
0879 % now fill in the flist array
0880 % [x y z d n e t] for each face
0881 flist=zeros(f,7);
0882 [ia,ib]=sort(edgeq(:,3));           % finst all the edges belonging to each face
0883 ix=[1; 1+find(ia(2:end)>ia(1:end-1)); 2*e+1];
0884 flist(:,6)=ib(ix(1:f));         % pointer to first edge for each face
0885 flist(:,5)=ix(2:end)-ix(1:end-1);       % size of each face
0886 % edgeq: 1=v1 2=v2 3=f1 4=f2 5=ev1 6=ef1 7=er 8=z 9=sf 10=sv]
0887 tedge=flist(:,6);           % this edge
0888 nedge=edgeq(tedge,6);       % next edge
0889 vmid=vlist(edgeq(tedge+2*e*(1-edgeq(tedge,10))),1:3);
0890 flist(:,1:3)=cross(vlist(edgeq(nedge+2*e*(1-edgeq(nedge,10))),1:3)-vmid,vlist(edgeq(tedge+2*e*edgeq(tedge,10)),1:3)-vmid,2);
0891 flist(:,1:3)=flist(:,1:3)./repmat(sqrt(sum(flist(:,1:3).^2,2)),1,3);
0892 flist(:,4)=sum(flist(:,1:3).*vmid,2);
0893 
0894 % now check for dual
0895 
0896 if dual
0897     vlisto=vlist;        % save vlist
0898     vlist=flist;         % old faces are new vertices
0899     flist=vlisto;
0900     % edgeq: 1=v1 2=v2 3=f1 4=f2 5=ev1 6=ef1 7=er 8=z 9=sf 10=sv]
0901     edgeq=edgeq(:,[3 4 2 1 6 5 7 8 10 9]);  % swap vertices and faces in the edge list
0902     erev=edgeq(:,7);        % reverse edge
0903     edgeq(:,6)=erev(edgeq(erev,6));
0904     flist(:,6)=erev(flist(:,6));
0905     if hemiQ             % cannot invert through centre
0906         error('Cannot take dual');
0907     else
0908         vlist(:,4)=vlist(:,4).^(-1);
0909     end
0910     vlist(:,1:3)=vlist(:,1:3).*repmat(vlist(:,4),1,3);
0911     v=size(vlist,1);
0912     vlist(:,1:3)=vlist(:,1:3)-repmat(mean(vlist(:,1:3),1),v,1);
0913     f=size(flist,1);            % recalculate face positions from scratch
0914     tedge=flist(:,6);           % this edge
0915     nedge=edgeq(tedge,6);       % next edge
0916     vmid=vlist(edgeq(tedge+2*e*(1-edgeq(tedge,10))),1:3);
0917     flist(:,1:3)=cross(vlist(edgeq(nedge+2*e*(1-edgeq(nedge,10))),1:3)-vmid,vlist(edgeq(tedge+2*e*edgeq(tedge,10)),1:3)-vmid,2);
0918     flist(:,1:3)=flist(:,1:3)./repmat(sqrt(sum(flist(:,1:3).^2,2)),1,3);
0919     flist(:,4)=sum(flist(:,1:3).*vmid,2);
0920 end
0921 
0922 info.snub=snubQ>0;
0923 info.onesided=onesidedQ>0;
0924 info.hemi=hemiQ>0;
0925 
0926 % create Wythoff string
0927 wystr='';
0928 j=0;
0929 for i=1:4
0930     if i==abs(wy(1))
0931         wystr=[wystr '| '];
0932         j=1;
0933     else
0934         if (wyd(i-j)==1)
0935             wystr=[wystr num2str(wyn(i-j)) ' '];
0936         else
0937             wystr=[wystr num2str(wyn(i-j)) '/' num2str(wyd(i-j)) ' '];
0938         end
0939     end
0940 end
0941 info.wythoff=wystr(1:end-1);
0942 info.vef=[v e f];
0943 info.chi=chi;   % Euler Characteristic
0944 
0945 % now draw an image if requested
0946 
0947 if ~nargout || any(md=='g')
0948     clf;
0949     if any(md=='t')
0950         [nii,nij]=rat(nimi(:,1));      % face sizes as rational numbers
0951         veca=0:1;
0952         plot(0,0,'ok');
0953         hold on
0954         for j=1:m
0955             i=rotsnub(j);
0956             veca=veca(2)*[1 exp(2i*gamma(i))];
0957             th=(0:nii(i))*2*pi/nimi(i,1);
0958             sc=[sin(th); 1-cos(th)];
0959             scv=sc(:,[2 end-1]);  % first and last vertices
0960             veci=veca/scv;
0961             plot(real(veci*sc), imag(veci*sc),'-k',real(veci*scv), imag(veci*scv),':b');
0962         end
0963         hold off
0964     else
0965         if (all(wyd==1) && ~any(md=='w')) || any(md=='f')
0966 
0967             % now render the polyhedron
0968 
0969             for i=1:f
0970                 %     patch('Faces',fa{i},'Vertices',vlist,'FaceVertexCData',rgb(i,:),'FaceColor','flat')
0971                 fa=zeros(flist(i,5),1);
0972                 ix=flist(i,6);
0973                 for j=1:numel(fa)
0974                     fa(j)=edgeq(ix,1+edgeq(ix,10));   % first vertex of this edge
0975                     ix=edgeq(ix,6); % next edge around the face
0976                 end
0977                 patch(vlist(fa,1),vlist(fa,2),vlist(fa,3),1-((flist(i,1:3)+1)/2).^2,'facealpha',0.95);
0978             end
0979         else
0980             edges=edgeq(edgeq(:,1)<edgeq(:,2),1:2);
0981             plot3(reshape(vlist(edges,1),size(edges))',reshape(vlist(edges,2),size(edges))',reshape(vlist(edges,3),size(edges))','k');
0982         end
0983         if any(md=='v')
0984             for i=1:v
0985                 text(vlist(i,1),vlist(i,2),vlist(i,3),sprintf('%d',i));
0986             end
0987         end
0988     end
0989     axis equal
0990     if dual
0991         tdual='Dual ';
0992     else
0993         tdual='';
0994     end
0995     if wyidx==round(wyidx) && wyidx>=1 && wyidx<=size(wys,1)
0996         title(sprintf('%s%s: %s (%s)',tdual,info.wythoff,names{wyidx,1:2}));
0997     else
0998         title(sprintf('%s%s',tdual,info.wythoff));
0999     end
1000 end

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003