Radius of inscribed sphere

To locate the edges of the Fourier box, insert the following code block immediately after your computation of G2 in ``setup.m'',

%# Locate edges (assume S's are even!) and determine max 'ok' G2
if any(rem(S,2)!=0)
  printf("Odd dimension in S, cannot continue...\n");
  return;
endif
eS=S/2+0.5;
edges=find(any(abs(M-ones(size(M,1),1)*eS')<1,2));

This block first verifies the assumption that the FFT box sizes are all even. After execution of this code block, edges contains a list of the indices of all rows of M with any element equal to S/2 or S/2+1, namely those with any component of $\vec G$ of largest possible absolute value.

Next, to find a list of the active rows of G (those corresponding to wave vectors within the Fourier sphere), first compute the minimum magnitude reciprocal lattice vector along the edges, and then find all vectors with magnitude less than one-half of that minimum magnitude (one-fourth of the square magnitude), by adding the code block

%# Compute active list and corresponding G2's
G2mx=min(G2(edges));
active=find(G2<G2mx/4); %# Sphere is 1/2 size (but looking at G^2!)
G2c=G2(active);

printf("Compression: %f (theoretical: %f)\n", ...
                     length(G2)/length(G2c), 1/(4*pi*(1/4)^3/3));

immediately after the block determining the edges. Note that this will produce in G2c a sequential list (without gaps!) of the values of $G^2$ associated with the rows listed in active.

Finally, be sure to make the values of G2c and active available to the various operators by getting up a global variable gbl_G2c in just the same way gbl_G and the other global variables are handled.