# v_upolyhedron

## PURPOSE

V_UPOLYHEDRON calculate uniform polyhedron characteristics

## SYNOPSIS

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

## DESCRIPTION

    V_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]=v_upolyhedron(w,md)
    0002 %V_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)); v_upolyhedron(i); pause(2); end
    0094 % for i=5:10, disp(num2str(i)); for j=1:5, v_upolyhedron(i+j/10); pause(2); end, end
    0095
    0096 %      Copyright (C) Mike Brookes 1997
    0097 %      Version: $Id: v_upolyhedron.m 10865 2018-09-21 17:22:45Z 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 Lesser General Public License as published by
    0105 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0112 %
    0113 %   You can obtain a copy of the GNU Lesser General Public License from
    0114 %   https://www.gnu.org/licenses/ .
    0115 %
    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);    % v_rotation to use first
    0747     v1 = vlist(i,1:3);    % the centre of v_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