function [V,C,dels,dems,XY]=localization(L,dom,N,J,rotb,anti)
% [V,C,dels,dems,XY]=LOCALIZATION(L,dom,N,J,rotb,anti)
%
% Returns bandlimited spectral eigenfunctions and their associated
% eigenvalues localized to a closed domain on the unit sphere.
%
% INPUT:
%
% L        Bandwidth, maximum angular spherical harmonic degree
% dom      'patch', 'sqpatch', 'africa',  'eurasia',  'namerica', 
%          'australia', 'greenland', 'samerica', 'amazon', 'orinoco',
%          'gpsnamerica', 'antarctica'
%          OR: [lon lat] an ordered list defining a closed curve [degrees]
% N        The splining smoothness [default: 10] or
%          The 3 or 4 parameters for 'patch' and 'sqpatch' passed to KERNELC,
%          go look there for the definitions
% J        Loads only so many eigenfuctions as required [default: all]
% rotb     1 rotate back as required, e.g. for Antarctica [default: 0]
%          but see now also KERNELC and KLMLMP2ROT
% anti     Get us the "anti-"region instead of the requested region
%
% OUTPUT:
%
% V        A list of all the eigenvalues 
% C        A cell array with cosine/sine coefficients eigenfunctions
% dels     Spherical harmonic degrees 
% dems     Spherical harmonic orders 
% XY       Coordinates of the outlines of the region
%
% EXAMPLE:
%
% [V,C,dels,dems]=localization(18,'amazon');
% plotplm([dels dems C{1}])
%
% th0=114; ph0=134; thR=40;
% [V,C,dels,dems]=localization(6,'patch',[th0 ph0 thR]*pi/180);
% plotplm([dels dems C{1}],[],[],4,2); hold on; plot(ph0,90-th0,'w')
%
% [V,C,dels,dems]=localization(60,'antarctica',[],12,1);
% plotplm([dels dems C{12}],[],[],2,0.5)
%
% SEE ALSO:
%
% KERNELC, SDWCAP, SDWCAP2, CDLKERNEL, SDWCAPT, SDWCAPT2, LOCALIZATION2D,
% PLOTPLM, PLM2XYZ, PLOTSLEP, KLMLMP2ROT
% 
% Last modified by fjsimons-at-alum.mit.edu, 11/22/2009

% And study covariance?

% Default inputs
defval('L',18)
defval('dom','australia')
defval('N',10)
defval('rotb',0)
defval('anti',0)

% Calculates the localization kernel for this domain
[Klmlmp,XY]=kernelc(L,dom,N);

defval('J',length(Klmlmp))

% I should build in the fact that if XY is just coordinates we should
% come up with a name to store the eigenfunctions.
if ~isstr(dom)
  doms=input...
       ('We still need a string for the filename:\n','s');
else
  doms=dom;
end

% See if the diagonalization has been done before
switch dom  
 case 'sqpatch'
  fnpl=sprintf('%s/%s-%i-%i-%i-%i-%i.mat',...
	       fullfile(getenv('IFILES'),'LOCALIZE'),doms,L,...
	       round(N(1)*180/pi),round(N(2)*180/pi),...
	       round(N(3)*180/pi),round(N(4)*180/pi),J);
 case 'patch'
  fnpl=sprintf('%s/%s-%i-%i-%i-%i-%i.mat',...
	       fullfile(getenv('IFILES'),'LOCALIZE'),doms,L,...
	       round(N(1)*180/pi),round(N(2)*180/pi),...
	       round(N(3)*180/pi),J);
 otherwise
  fnpl=sprintf('%s/%s-%i-%i-%i.mat',...
	       fullfile(getenv('IFILES'),'LOCALIZE'),doms,L,N,J);
end

if anti==1
  Klmlmp=eye(size(Klmlmp))-Klmlmp;
  fnpl=sprintf('%s/%s-%i-%i-%i-%i.mat',...
	       fullfile(getenv('IFILES'),'LOCALIZE'),doms,L,N,J,anti);
end

if exist(fnpl,'file')==2
  disp(sprintf('Loading %s',fnpl))
  load(fnpl)
else
  % Calculates the eigenfunctions/values for this localization problem
  if J<length(Klmlmp)-1 & 1==3
    % Note that this often does not work very well at all, so don't try
    OPTS.disp=0;
    [C,V]=eigs(Klmlmp,J,'LA',OPTS);
    [V,isrt]=sort(sum(real(V),1),'descend');
    C=C(:,isrt(1:J));
  else
    % This is slow but better
    [C,V]=eig(Klmlmp);
    [V,isrt]=sort(sum(real(V),1),'descend');
    C=C(:,isrt(1:J));
    save(fnpl,'C','V','dom','L','N','J')
  end
end
  
% Sticks the cosine/sine coefficients back
% into the right place of LMCOSI
% Preallocate/initialize cell array
[dems,dels,mz,lmc,mzin]=addmon(L);
%CC=cell(1,J);
%[CC{:}]=deal(nan(length(dels),2));
CC=cellnan(J,length(dels),2);
for index=1:J
  CC{index}=reshape(insert(C(:,index),0,mzin),2,length(dems))';
end

% You can plot this here, if you want, by doing, e.g.
% plotplm([dels dems CC{10}],[],[],4,1); hold on
% [X,Y,Z]=sph2cart(XY(:,1)*pi/180,XY(:,2)*pi/180,1);
% plot3(X,Y,Z)
% But for Antarctica, this will show up in the "wrong" spot

% It it's Antarctica, need to rotate it back in shape
if strcmp(dom,'antarctica') && rotb==1
  [XY,lonc,latc]=antarctica;
  for index=1:J
    CC{index}=kindeks(plm2rot([dels dems CC{index}],...
			      lonc,-latc,0),3:4);
  end
end

% You can plot this here, if you want, by doing, e.g.
% plotplm([dels dems CC{10}],[],[],2,0.5)
% This now does show up in the right spot

% Make standard output
V=V(:);
C=CC;