function [G,V,EL,EM,N,GM2AL,MTAP,IMTAP]=glmalpha(TH,L,sord,blox,upco,resc)
% [G,V,EL,EM,N,GM2AL,MTAP,IMTAP]=GLMALPHA(TH,L,sord,blox,upco,resc)
%
% Returns an (lm)X(alpha) matrix with spherical harmonic coefficients of
% the BANDLIMITED Slepian functions of the SINGLE or DOUBLE polar cap, or
% of a geographical region of interest. Only in this case are the
% eigenvalues automatically sorted. The matrix G is orthogonal, G'*G is
% the identity.
%
% INPUT:
%
% TH       Angular extent of the spherical cap, in degrees OR
%          'england', 'eurasia',  'namerica', 'australia', 'greenland'
%          'africa', 'samerica', 'amazon', 'orinoco'
% L        Bandwidth (maximum angular degree)
% sord     1 Single polar cap of diameter TH [default]
%          2 Double polar cap of diameter TH (each, and watch out!)
%          N Splining smoothness for geographical regions [default: 10]
%
% The following options are only for axisymmetric polar caps. 
%
% blox     0 Standard (lm) ordering, l=0:L, m=-l:l [default]
%          1 Block-diagonal ordering, m=[0 -1 1 -2 2 ... -L L], l=abs(m):L
% upco     +ve fraction of unit radius for upward continuation [default: 0]
%          -ve fraction of unit radius for downward continuation
% resc     0 Not rescaled [default]
%          1 Rescaled to have a unit integral over the unit sphere
%            (this is only relevant for the down/upward continued functions)
%
% OUTPUT:
%
% G        The unitary matrix of localization coefficients
% V        The eigenvalues in this ordering (not automatically sorted)
% EL       The vector with spherical harmonic degrees as first index of G
% EM       The vector with spherical harmonic orders as first index of G
% N        The Shannon number
% GM2AL    The sum over all orders of the squared coefficients, i.e. the
%          TOTAL power, NOT the power spectral density
% MTAP     The order of the eigentapers, if the region is axisymmetric
% IMTAP    The rank within that particular order of the eigentapers
%
% EXAMPLE:
%
% Note that using ADDMOUT you can get this back to block-diagonal form
% G=glmalpha; [a,b,c,bl]=addmout(18); imagesc(G(bl,:))
% difer(G'*G-eye(size(G)))
%
% SEE ALSO:
%
% ADDMOUT, ADDMON, KERNELC, LOCALIZATION, GALPHA, DLMLMP
%
% Last modified by fjsimons-at-alum.mit.edu, 04/13/2007

% Should be able to update this to retain the rank order per m as well as
% the global ordering

% Does this work for the whole-sphere? In that case, should really want G
% to be the identity - all though of course, anything works, too. You 
% don't get necessarily the spherical harmonics back... 

defval('TH',30)
defval('L',18)
defval('blox',0)
defval('upco',0)
defval('resc',0)

defval('mesg','GLMALPHA Check passed')
% Hold all messages
mesg=NaN;

% Just get the file names here
if upco==0 & resc==0
  if ~isstr(TH) % POLAR CAPS
    defval('sord',1) % SINGLE OR DOUBLE CAP
    fname=fullfile(getenv('IFILES'),'GLMALPHA',...
		   sprintf('glmalpha-%i-%i-%i-%i.mat',TH,L,sord,blox));
    % Initialize ordering matrices
    MTAP=repmat(0,1,(L+1)^2);
    IMTAP=repmat(0,1,(L+1)^2);

  else % GEOGRAPHICAL REGIONS
    defval('sord',10) % SPLINING SMOOTHNESS
    fname=fullfile(getenv('IFILES'),'GLMALPHA',...
		   sprintf('glmalpha-%s-%i.mat',TH,L));
    defval('MTAP',NaN) % If not, calculate order per taper
    defval('IMTAP',NaN) % And rank ordering within that taper
    defval('xver',0) % For excessive verification of the geographical case
  end
else
  fname='neveravailable';
  defval('xver',0) % For excessive verification of the upco'd case
end

if exist(fname,'file')==2 % & 1==3
  load(fname)
  % disp(sprintf('Loading %s',fname))
else
  % Initialize matrices
  G=repmat(0,(L+1)^2,(L+1)^2);
  V=repmat(0,1,(L+1)^2);
  
  % Find indices into G belonging to the orders
  [EM,EL,mz,blkm]=addmout(L);
  
  % Find increasing column index; that's how many belong to this order
  alpha=cumsum([1 L+1 gamini(L:-1:1,2)]);
    
  % For GEOGRAPHICAL REGIONS
  if isstr(TH)
    % Calculates the localization kernel for this domain
    [Klmlmp,XY]=kernelc(L,TH,sord);
    
    % Calculates the eigenfunctions/values for this localization problem
    [G,V]=eig(Klmlmp);
    [V,isrt]=sort(sum(real(V),1));
    V=fliplr(V);
    G=G(:,fliplr(isrt));
    
    [a,b,c,d,e,f,ems,els,R1,R2]=addmon(L);
    % This indexes the orders of G back as 0 -101 -2-1012 etc
    G=G(R1,:);
    % Check indexing
    difer(els(R1)-EL,[],[],mesg)
    difer(ems(R1)-EM,[],[],mesg)
    
    % Calculate Shannon number and compare this with the theory
    N=sum(V);

    % Is the Shannon number right?
    difer((L+1)^2*Klmlmp(1)-N,[],[],mesg)

    % Check if the expansion of a basis function is indeed either 1 or 0
    if xver==1
      disp('Excessive verification')
      % Is the area right?
      difer(Klmlmp(1)-spharea(TH),4,[],mesg)
      % This is a bit double up... but it's only for excessive verification
      [V1,C]=localization(L,TH,sord);
      difer(V-V1',[],[],mesg)
      for index=1:length(C)
	salpha=G'*C{index}(R2);
	% Only one of these functions should get "hit"
	difer(sum(abs(salpha)>1e-9)-1,[],[],mesg)
      end
    end
  else
    % For AXISYMMETRIC REGIONS
    if blox~=0 & blox~=1
      error('Specify valid block-sorting option ''blox''')
    end
    % For the SINGLE or DOUBLE POLAR CAPS
    for m=0:L
      % Same results for +/- m
      if sord==1
	[E,Vg,th,C,T,Vp]=grunbaum(TH,L,m,[]);
      elseif sord==2
	[E,Vg,th,C,T,Vp]=grunbaum2(TH,L,m,[]);
      else
	error('Specify single or double polar cap')
      end
      if upco~=0
	if upco>0
	  % The upward continuation matrix
	  A=diag((1+upco).^[-(m:L)-1]);
	elseif upco<0
	  % The downward continuation matrix
	  A=diag((1+abs(upco)).^[(m:L)+1]);
	end
	
	if xver==1
	  % This should give the same result, more or less, less accurate 
	  if sord==1
	    [a,Vs,c,d,Cs,e,f,g,h,j,D]=sdwcap(TH,L,m,0,-1);
	  else
	    [a,Vs,c,Cs,e,f,D]=sdwcap2(TH,L,m,0,-1);
	  end
	  % This should give the eigenvalues again, which we'd had from
          % orthocheck 
	  warning off
	  % Check difference integration eigenvalues and those from
          % kernel 
	  difer(Vp(:)-diag((C'*D*C)./(C'*C)),[],[],mesg)
	  % Check difference integration and diagonalization eigenvalues
	  difer(Vp(:)-Vs(:),[],[],mesg)
	  warning on
	  Vc=diag((C'*A*C*diag(Vp)*C'*A*C));
	  Vp0=Vp;
	end
	
	% Upward continuation from 1 to 1+a or from 1+a to 1:
	% New eigenfunctions, same name
	C=A*C;
	% Calculate new eigenvalues, same name
	[jk1,jk2,jk3,Vp,nofa]=orthocheck(C,[],TH/180*pi,m,sord,1);

	% Make sure these are sorted, since that's not automatically the case
	% [Vp,ind]=sort(Vp,'descend');
	% C=C(:,ind);
	% Current thinking is: do NOT resort, as you'll want to compare the
	% best at a=0 with whatever it becomes later!
	
	if xver==1
	  warning off
	  % Check difference integration eigenvalues and those from kernel
	  difer(Vp(:)-diag((C'*D*C)./(C'*C)),[],[],mesg)
	  warning on
	  % Check how many Vp>Vp0 
	  disp(sprintf('%i/%i eigenvalues greater',sum(Vp(:)>Vp0(:)), ...
		       length(Vp0)))
	  disp(sprintf('Shannon number greater: %i',sum(Vp)>sum(Vp0)))
	end
	if resc==1
	  % Rescale these functions to have an integral to unity over the
	  % sphere; note: this doesn't make the set orthonormal of course
	  C=C*diag(1./nofa);
	end
      end
      
      % Distribute this at the right point in the huge matrix
      if m>0
	G(EM==-m,alpha(2*m):alpha(2*m+1)-1)=C;
	V(alpha(2*m):alpha(2*m+1)-1)=Vp;
	MTAP(alpha(2*m):alpha(2*m+1)-1)=-m;
	% It's all neatly ordered here, downgoing within every order
	IMTAP(alpha(2*m):alpha(2*m+1)-1)=1:length(Vp);
      end
      % Duplicate for the other order in case the region is axisymmetric
      G(EM==m,alpha(2*m+1):alpha(2*m+2)-1)=C;
      V(alpha(2*m+1):alpha(2*m+2)-1)=Vp;
      MTAP(alpha(2*m+1):alpha(2*m+2)-1)=m;
      % It's all neatly ordered here, downgoing within every order
      IMTAP(alpha(2*m+1):alpha(2*m+2)-1)=1:length(Vp);;
    end
    % Calculate the Shannon number and compare it to the theory
    N=sum(V);
    difer(N-(L+1)^2*sord*(1-cos(TH/180*pi))/2,[],[],mesg);
    
    % Compute the sum over all orders of the squared coefficients
    % Thus works when they have not been blocksorted yet. 
    GM2AL=repmat(0,(L+1)^2,L+1);
    for l=0:L
      b=(l-1+1)^2+1;
      e=(l+1)^2;
      GM2AL(:,l+1)=sum(G(b:e,:).^2,1)';
    end
    % Make sure that the sum over all degrees is 1
    difer(sum(GM2AL,2)-1,[],[],mesg)

    % This is not blockdiagonal, unless you make it thus
    if blox==1
      G=G(blkm,:);
      EM=EM(blkm);
      EL=EL(blkm);
    end
  end
  % Save the results
  eval(sprintf(...
      'save %s G V EL EM N GM2AL MTAP IMTAP',fname))
end