function [m,rimnodes] = newcirclemesh( sz, circum, nrings, centre, inner, dealign, sector, coneangle )
%[m,rimnodes] = newcirclemesh( sz, circum, nrings, centre, inner, dealign, sector, coneangle )

    if nargin < 8
        if sz(3) ~= 0
            coneangle = pi/2;
        else
            coneangle = 0;
        end
    end
    if nargin < 7
        sector = 1;
    end
    if nargin < 6
        dealign = false;
    end
    if nargin < 5
        inner = 0;
    end
    if nargin < 4
        centre = [0,0,0];
    end
    if (nargin < 2) || (circum < 4)
        circum = 4;
    end
    if (nargin < 3) || (nrings <= 0)
        nrings = max(1,floor(circum/6));
    end
    if nrings==1
        circums = circum;
    else
        if inner==0
            inner = floor(circum/nrings);
        end
        inner = max( [ 3, ceil(4*sector), inner ] );
        circums = arithprog( inner, circum, nrings );
    end
    if nargin < 1
        sz = [2 2 1];
    end
    r = sqrt(sz(1)*sz(2))/2;
    if length(sz)==2 || (sz(3)==0)
        if coneangle==0
            sz(3) = 0;
        else
            if abs(coneangle) < 0.001
                h = coneangle*r/2;
            else
                h = r*(1-cos(coneangle))/sin(coneangle);
                if coneangle < 0
                    h = -h;
                end
            end
            sz(3) = h;
            coneangle = abs(coneangle);
        end
    elseif (sz(3) ~= 0) && (coneangle == 0)
        coneangle = pi/2;
    end
    
    pts = [0.0 0.0 sz(3)];
    np = 1;
    if dealign
        offsets = zeros( 1, nrings );
        for i=(nrings-1):-1:1
            m = circums(i);
            n = circums(i+1);
            g = gcd(m,n);
            desynch = randint(1,n)*m/n;
            desynch = desynch - floor(desynch);
            offsets(i) = ((-1)^i)*g/(2*n) + desynch;
        end
      % offsets
    else
        offsets = zeros( 1, nrings );
    end
    angles = cell(1,nrings);
    sectorAngle = sector*pi*2;
    if coneangle ~= 0
        csc_ca = 1/sin(coneangle);
        cos_ca = cos(coneangle);
    end
    for i=1:nrings
        cn = circums(i);
        if sector==1
            angleStep = sectorAngle/cn;
        else
            angleStep = sectorAngle/(cn-1);
        end
        angles{i} = (((0:(cn-1)) + offsets(i))*angleStep)';
        if sz(3) ~= 0
            phi = coneangle*(i/nrings);
            h0 = csc_ca*(1-cos_ca);
            cphi = cos(phi);
            sphi = sin(phi);
            hphi = csc_ca*(cphi-cos_ca);
            rphi = csc_ca*sphi;
            pts((np+1):(np+cn),[1 2]) = ...
                [ cos(angles{i}), sin(angles{i}) ] * rphi;
            pts((np+1):(np+cn),3) = sz(3) * hphi / h0;
        else
            radius = single(i)/nrings;
            pts((np+1):(np+cn),[1 2]) = ...
                [ cos(angles{i}), sin(angles{i}) ]*radius;
        end
        np = np+cn;
    end
    tri = zeros(0,3);
    if sector==1
        tri(1:circums(1),1:3) = [ ones(circums(1),1), ...
                                  (2:(circums(1)+1))', ...
                                  [(3:(circums(1)+1))'; 2] ];
    else
        tri(1:(circums(1)-1),1:3) = [ ones(circums(1)-1,1), ...
                                  (2:(circums(1)))', ...
                                  (3:(circums(1)+1))' ];
    end
    pis_inner = 2:(circums(1)+1);
    for i = 2:nrings
        outerStart = pis_inner(length(pis_inner))+1;
        pis_outer = outerStart : (outerStart - 1 + circums(i));
        iistart = 1;
        ii = iistart;
        imax = length(pis_inner);
        oistart = 1;
        oi = oistart;
        omax = length(pis_outer);
        detectend = false;
        while true
            if (sector ~= 1) && ((ii==imax) || (oi==omax))
                break;
            end
            ii1 = ii+1; if ii1 > imax, ii1 = 1; end
            oi1 = oi+1; if oi1 > omax, oi1 = 1; end
            outerPt = pts( pis_outer(oi), : );
            innerPt = pts( pis_inner(ii), : );
            outerPt1 = pts( pis_outer(oi1), : );
            innerPt1 = pts( pis_inner(ii1), : );
            outerDiag = norm( outerPt1 - innerPt );
            innerDiag = norm( outerPt - innerPt1 );
            nextOuter = outerDiag <= innerDiag;
            if nextOuter
                nextIndex = pis_outer(oi1);
            else
                nextIndex = pis_inner(ii1);
            end
            tri(size(tri,1)+1,:) = [ pis_inner(ii), pis_outer(oi), nextIndex ];
            if nextOuter
                oi = oi1;
            else
                ii = ii1;
            end
          % fprintf( 1, 'ii %d oi %d\n', ii, oi );
            if detectend && ((ii==iistart) || (oi==oistart))
                break;
            end
            if (ii ~= iistart) && (oi ~= oistart)
                detectend = true;
            end
        end
        if sector ~= 1
            if ii==imax
                while oi < omax
                    oi1 = oi+1;
                    tri(size(tri,1)+1,:) = [ pis_inner(ii), pis_outer(oi), pis_outer(oi1) ];
                    oi = oi1;
                end
            else
                % oi==omax
                ii1 = ii+1;
                tri(size(tri,1)+1,:) = [ pis_inner(ii1), pis_inner(ii), pis_outer(oi) ];
            end
        else
            if ii==iistart
                % oi ~= oistart
              % fprintf( 1, 'ii == iistart == %d, oi %d oistart %d\n', ...
              %     ii, oi, oistart );
                while oi ~= oistart
                    oi1 = oi+1;  if oi1 > omax, oi1 = 1; end
                    tri(size(tri,1)+1,:) = [ pis_inner(ii), pis_outer(oi), pis_outer(oi1) ];
                    oi = oi1;
                end
            else
                % oi == oistart
              % fprintf( 1, 'oi == oistart == %d, ii %d iistart %d\n', ...
              %     oi, ii, iistart );
                while ii ~= iistart
                    ii1 = ii+1;  if ii1 > imax, ii1 = 1; end
                    tri(size(tri,1)+1,:) = [ pis_inner(ii1), pis_inner(ii), pis_outer(oi) ];
                    ii = ii1;
                end
            end
        end
        
        pis_inner = pis_outer;
    end
    pts(:,1) = pts(:,1) * (sz(1)/2);
    pts(:,2) = pts(:,2) * (sz(2)/2);
    pts = pts + repmat( centre, size(pts,1), 1 );
    m = struct( 'nodes', pts, 'tricellvxs', tri );
    m.globalProps.trinodesvalid = true;
    m.globalProps.prismnodesvalid = false;
    
    numnodes = size(m.nodes,1);
    rimnodes = [ (numnodes-circums(nrings)+1); ...
                 (numnodes:-1:(numnodes-circums(nrings)+2))' ];
end
