Issue with C2h symmetry for sphgrid and workaround

[Murad]

Hi,

I've looked through the forum and couldn't see if this was previously reported.

When running sphgrid using C2h symmetry the values for phi are [0,pi); whereas when running sphtri the values of phi have been assumed to be [0,pi]. This means when trying to run the code from http://easyspin.org/documentation/sphtri.html

Code: Select all

x = sphgrid('D2h',nKnots);
g2 = [1.9*x(1,:); 2*x(2,:); 2.1*x(3,:)].^2;
g = sqrt(sum(g2));

tri = sphtri('D2h',nKnots);
trisurf(tri,x(1,:),x(2,:),x(3,:),g);
view(135,39); colorbar

generates a warning and does not plot. This is because the largest index generated by sphtri is (nKnots-1) larger than the number of points generated by sphgrid. So, I've added in a few lines to run after sphgrid to add in phi=pi at each knot (except the first one where theta=0). Hope this helps:

Code: Select all

% want to add in phi=pi for each knot: error is easyspin
nKnots=22;

[phi,theta]=sphgrid('C2h',nKnots);
[th_unique]=unique(theta);
th_unique=th_unique(2:end); % ignore first point because theta=0, so phi doesn't matter
insertion=[th_unique;pi*ones(1,numel(th_unique))]'; % points to be inserted
    
x_polar=[[theta',phi'];insertion];
[~,ind]=sort(x_polar(:,1));
x_polar=x_polar(ind,:); %x_polar gives the correct theta and phi for sphtri

for n= 1:numel(x_polar(:,1))
   x(:,n)=[sin(x_polar(n,1))*cos(x_polar(n,2));sin(x_polar(n,1))*sin(x_polar(n,2));cos(x_polar(n,1))];
end
% x gives the correct Cartesian values for sphtri

% now copying the code from http://easyspin.org/documentation/sphtri.html
g2 = [1.9*x(1,:); 2*x(2,:); 2.1*x(3,:)].^2;
g = sqrt(sum(g2));

tri = sphtri('C2h',nknots);
trisurf(tri,x(1,:),x(2,:),x(3,:),g);
view(135,39); colorbar