function [E,V,N,th,C,ngl1,ngl2,unc,com,sdl,K]=sdwcap(TH,L,m,nth,vcut,grd)
% [E,V,N,th,C,ngl1,ngl2,unc,com,sdl,K]=SDWCAP(TH,L,m,nth,vcut,grd)
%
% Spherical harmonic localization to a SINGLE spherical polar cap: 
% bandlimited and optimally spatially concentrated solutions.
%
% INPUT:
%
% TH          Angular extent of the spherical cap, in degrees
% L           Bandwidth (maximum angular degree)
% m           Angular order of the required data window, -l<m<l
% nth         Number of colatitudes to evaluate (not if 0) 
% vcut        Cut-off eigenvalue [default= eps*10]
%             If vcut==-1 automatically recalculates everything
% grd         1 Colatitudes only; returns matrix E
%             2 Colatitude/Longitude; returns cell E
%
% OUTPUT:
%
% E           Optimally concentrated tapers, expanded to space domain 
% V           Eigenvalues, sorted
% N           Sum of all the eigenvalues (Shannon number)
% th          Colatitudes at which the functions are evaluated, in degrees
% C           Optimally concentrated tapers, spherical harmonics
%             coefficients, normalized to unity on the unit sphere
% ngl1, ngl2  Number of GL points on unit sphere and domain, respectively
% unc         Uncertainty product
% com         Center of mass
% sdl         Spherical harmonics standard deviation
% K           The matrix that is diagonalized to produce the eigenfunctions
%
% See also DOUBLECAP, GRUNBAUM
%
% Last modified by fjsimons-at-alum.mit.edu, 04.08.2005

defval('TH',40)
defval('L',18)
defval('m',0)
defval('nth',720)
defval('vcut',eps*10)
defval('grd',1)

% Work with the absolute value of m
mor=m;
m=abs(m);

if(m>L)
  error('Order cannot exceed degree')
end

% Filename of saved data
fnpl=sprintf('%s/SDW-%i-%i-%i-%i.mat',...
	     fullfile(getenv('IFILES'),'SDWCAP'),TH,L,nth,m);

% In 7.0.0.19901 (R14), in very rare cases (*-180-0-0) this
% save is buggy, and cannot be subsequently loaded
% Used matzerofix to remediate this
if exist(fnpl,'file')==2 & ~(L==180 & m==0) & (vcut~=-1)
  % eval(sprintf('load %s',fnpl))
  load(fnpl)
  disp(sprintf('%s loaded by SDWCAP',fnpl))
else
  N=(L+1)^2*(1-cos(TH/180*pi))/2;
  
  % Approximate Nyquist degree is $\pi/d\theta$
  if nth < (L+1) & nth~=0
    error('Sample finer to avoid aliasing')
  end
  
  % Convert to radians
  TH=TH*pi/180;

  % Initialize kernel
  K=repmat(NaN,L+1-m,L+1-m);
  % Construct kernel: integral of Ylm over the cap
  for lr=m:L
    for lc=lr:L
      % Orthonormalization of Ylm is to unity over unit sphere
      % When TH=180, K should be the identity matrix
      % The pi*(1+(m==0)) is the longitudinal integral
      K(lr+1-m,lc+1-m)=legendreprodint(lr,m,lc,m,cos(TH))...
	  *sqrt(2*lr+1)*sqrt(2*lc+1)/(4*pi)*pi*(1+(m==0));
      % Symmetrize
      K(lc+1-m,lr+1-m)=K(lr+1-m,lc+1-m);
    end
  end

  % Calculate eigenvalues and eigenvectors; C'*C=I
  [C,V]=eig(K);
  
  % Check normalization and get number of GL points
  [ngl1,ngl2,com]=orthocheck(C,V,TH,m);

  if length(C)>1 & m==0
    % Calculate uncertainty relation
    sdl=sqrt(sum(repmat([0:L]'.*[1:L+1]',1,size(C,2)).*C.^2,1));
    N=1;
    sdv=sqrt(N*(1+(1-2/N)*com.^2));
    unc=sdv./com.*sdl;
  else
    com=0;
    unc=0;
    sdl=0;
    sdv=0;
  end

  % Order eigenvalues and eigenfunctions downgoing
  [V,isrt]=sort(sum(V,1));
  V=fliplr(V);
  C=C(:,fliplr(isrt));
  unc=fliplr(unc); sdl=fliplr(sdl); com=fliplr(com); sdv=fliplr(sdv); 

  % Only return nonzero "useful" eigenvalues
  C=C(:,V>vcut); 
  V=V(V>vcut);

  % Compute spatial functions, colatitudinal part only
  if nth~=0
    % Zonal functions only 
    if m==0
      % Make spatial functions
      % This is SDW (2005) equation (5.10) combined with the sqrt(2-dom) of
      % (5.12) already included!
      [E,th]=pl2th(C,nth,1);
      th=th*180/pi;
      nlon=2*nth-1;
    else
      % This is SDW (2005) equation (5.10) combined with the sqrt(2-dom) of
      % (5.12) already included!
      [E,nlon,lat]=plm2th(C,nth,m,1);
      th=linspace(0,180,size(E,1));
    end    
    % Make E start with a positive lobe and ajust C too
    % Don't take first sample as for m~=0 it is 0
    for index=1:size(E,2)
      C(:,index)=C(:,index)*sign(E(2,index));
      E(:,index)=E(:,index)*sign(E(2,index));
    end
  else
    E=0;
    th=0;
    nlon=0;
  end
  % In 7.0.0.19901 (R14), in very rare cases (3-180-0-0) this
  % save is buggy, and cannot be subsequently loaded
  % Output in degrees
  eval(sprintf('save %s E V L N TH C th nlon ngl1 ngl2 unc com sdl K',fnpl))  
end

if nth~=0 & grd==2
  % Output on full grid; watch the sign of m
  % Note that sqrt(2-dom) is already part of the harmonic
  if mor<=0
    EE=E; clear E
    for index=1:size(EE,2)
      E{index}=EE(:,index)*cos(m*linspace(0,2*pi,nlon));
    end
  end
  if mor>0
    EE=E; clear E
    for index=1:size(EE,2)
      E{index}=EE(:,index)*sin(m*linspace(0,2*pi,nlon));
    end
  end
end