sphgrid.m: Wrong number of triangles

[florianquintes]

When calculating the Delaunay triangles, an error occurs in the case of point groups with openphi==True. The number of triangles depends on the number of octants and nodes of the grid, whereby (nKnots-1)**2 triangles are obtained for one octant.

As can be seen for D2h, this is true for openphi==False. For the example C4h, however, too many triangles are obtained. Several errors lead to this:

1. for triangles near the pole, due to the symmetry conditions, two points of the triangle have the identical angles, so that a line is obtained instead of a triangle (e.g. triangle [1 2 2]). These triangles are then sorted out on the basis of their area (line 190-192), which means that these angles do not (fully) contribute to the spectrum.
2. since too many triangles are obtained despite the sorting out of "triangles" with an area of 0, there must still be triangles that are obtained more than once. I have not yet investigated this in more detail, for this you would have to look at the indices of all triangles obtained.

This bug should only have a minor effect, as the majority of triangles are correct and the proportion decreases as the grid size increases.

D2h:

Sys = struct('S',1/2,'g',[2 2 2]);
Sys.lw = 0.5;
Exp = struct('CenterSweep',[325 80],'mwFreq',9.5,'Harmonic',1);
Opt = struct('Verbosity',2);
Opt.GridSize = [10 0];
Opt.GridSymmetry = 'D2h'

[x,y1] = pepper(Sys,Exp,Opt);

Output:

-absorption spectrum construction---------------------- interpolation off triangle/segment projection spectrum array size: 1x1024 total 81 triangles (55 orientations), 1 transitions

D2h_matlab.pdf

C4h

Sys = struct('S',1/2,'g',[2 2 2]);
Sys.lw = 0.5;
Exp = struct('CenterSweep',[325 80],'mwFreq',9.5,'Harmonic',1);
Opt = struct('Verbosity',2);
Opt.GridSize = [10 0];
Opt.GridSymmetry = 'C4h'

[x,y1] = pepper(Sys,Exp,Opt);

Output:

-absorption spectrum construction---------------------- interpolation off triangle/segment projection spectrum array size: 1x1024 total 85 triangles (46 orientations), 1 transitions

C4h.pdf

Both plots were obtained by adding the following code in sphgrid.m (line 190-194):

  rmv = Areas==0;
  Tri(rmv,:) = [];
  phx = [phi maxPhi*ones(1,GridSize-1)];           % added lines
  thx = [theta pi/2*(1:GridSize-1)/(GridSize-1)];  % copied from l. 139-142
  x = thx.*cos(phx);                               %
  y = thx.*sin(phx);                               %
  triplot(Tri, x, y)                               % Plot function
  Areas(rmv) = [];
  Tri = sort(Tri,2);
  Tri = uint32(Tri);