% projective prolet

disp('Takes about 54 minutes to do NW=12 (24 calcs) on a Decstation 5000/200')
disp('could speed up by not interpolating for the zeros')

begtime = fix(clock)

NW = 12;

M = 100;    % size of final output image
N = 100;

if ~exist('V')
  [V d] = slepian(200,NW);
end%if

SS_all = zeros(M,N);   % size of image to be 100 by 100
t = linspace(-1,1,600).';
tproj = project(t+1.34,pi/4,1,2); % 1.34 chosen to keep within [-1,1] and get 
                                  % projlet close to center

for nw     = 1:2*NW   % only bother with the big fish
%for nw     = 24:24 % diagnostic quick test
  s_part   = V(:,nw); % slepian is one column of V
  s        = [zeros(s_part);s_part;zeros(s_part)];
  c        = chirp(-.25,.25,length(s)); % shear in freq by .3
  s        = s(:).*c(:);
chirp(-.25,.25,length(s)); % shear in freq by .3
  disp(sprintf('%g  <---------------------------- PROJECTING SLEPIAN %g',nw,nw))
  s_projt = interp2(t(:),s(:),tproj(:));  % interpolate new projected pts
  S_projt = fftshift(fft(s_projt(:)));
  C = chirp(.25,-.25,length(S_projt));
  S_projt = S_projt(:).*C(:);% chirp_in_freq=shear_in_time
  S_projtf = interp2(t(:),S_projt(:),tproj(:));%  interpolate but this time in F
  s_projtf = ifft(fftshift(S_projtf));
  SS_temp  = tf(s_projtf,M,N);
  SS_all = SS_all + SS_temp.*conj(SS_temp);  % must use SQR magnitude
end%for

endtime = fix(clock)
disp('took this much time to generate the one projected prolet:')
endtime-begtime

disp('now about to display the projected prolet')
dispmatb(SS_all)

