``cI.m''

In lecture, we found that the I() operator becomes ${\bf I}={\bf {\bf F}_{(3)}} {\bf B}$, where ``${\bf F}_{(3)}$'' is the three-dimensional Fourier transform kernel and ``${\bf B}$'' maps data from sphere of active basis functions into the Fourier box. To implement this, include ``global gbl_active'' at the top of cI(). Note that, regardless of the number of rows of the input data, the output always has the same size of data, the full FFT box. Thus, you should change your initial allocation of the output to

out=zeros(prod(gbl_S),size(in,2));

Finally, to perform the mapping ${\bf B}$ when needed, you should use a code block along the lines of

if size(in,1)==prod(gbl_S)
  out(:,col)= ... %# previous code
else
  full=zeros(prod(gbl_S),1); full(gbl_active)=in(:,col);
  out(:,col)= ... %# previous code with "full" replacing "in(:,col)"
end

Note that, here and below, ``along the lines of'' indicates that the names of some variables (especially internal working variables of a function which you designed) in the example may not correspond exactly to those in your routine.