function [K2,K1,K]=klmlmp2rot(Klmlmp,lonc,latc)
% [K2,K1,K]=KLMLMP2ROT(Klmlmp,lonc,latc)
%
% 
% INPUT:
%
% Klmlmp       A localization kernel coming out of, e.g. KERNELC
% lonc         A "longitudinal" rotation parameter 
% latc         A "latitudinal" rotation parameter
%
% OUTPUT:
%
% K2           The same kernel rotated such that its eigenfunctions are
%              rotated versions of the eigenfunctions of the unrotated
%              kernel
% K1           The same kernel in an intermediate step that should reveal
%              the geometry of the concentration region when plotted
% K            This should be the same as the original kernel down to eps
%
% EXAMPLE:
%
% 
%
% Last modified by fjsimons-at-alum.mit.edu, 09/28/2009

% Here, we will rotate the kernel rather than rotate the solutions in
% LOCALIZATION, after the fact, which should be perfectly equivalent. The
% advantage here is that we can combine this with other localization
% kernels that do not require such a special counterrotation procedure.

% Get the maximum degree for this square, symmetric kernel
Lmax=sqrt(size(Klmlmp,1))-1;
% Get the right indices ready
[dems,dels,mz,lmc,mzin,mzo]=addmon(Lmax);

% Apply the rotation on the ROWS
alp=lonc;
bta=-latc;
gam=0;
% Do this once
disp('First rotation')
K1=rotateonce(Klmlmp,alp,bta,gam,dems,dels,mzin,mzo);

% For a near-zonal function the structure should now be apparent

% Check the norm, it should not have changed
difer(norm3d(K1)-norm3d(Klmlmp),[],[],NaN)
% Nor should the norm per column have changed
difer(sum(K1.^2,1)-sum(Klmlmp.^2,1),[],[],NaN)

% Apply the rotation on the COLUMNS thus the transpose
disp('Second rotation')
K2=rotateonce(K1',alp,bta,gam,dems,dels,mzin,mzo);

% Check the norm, it should not have changed
difer(norm3d(K2)-norm3d(K1),[],[],NaN)
% Nor should the norm per column have changed
difer(sum(K2.^2,1)-sum(K1'.^2,1),[],[],NaN)

% Check the symmetry of the output
difer(K2-K2',[],[],NaN)

% Only if you want to test
if nargout>2
  % Apply a DIFFERENT rotation on the ROWS of K1
  alp=0;
  bta=latc;
  gam=-lonc;
  disp('Third rotation (verification round)')
  K=rotateonce(K1,alp,bta,gam,dems,dels,mzin,mzo);
  % And check that you get the input back
  difer(K-Klmlmp,[],[],NaN)
  % but it ain't going to be identical, for sure
  % minmax(K-Klmlmp)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function K=rotateonce(K,alp,bta,gam,dems,dels,mzin,mzo)
% Initialize temporary structure
Ktmp=cellnan(size(K,1),length(dels),2);

tic
% First, rearrange the ROWS into LMCOSI format for each COLUMN
for index=1:size(K,1)
  Ktmp{index}=reshape(insert(K(:,index),0,mzin),2,length(dems))';
end
% Now perform the counterrotation as per alp, bta and gam
% Should really write something that avoids having to loop here
h=waitbar(0,sprintf('Performing spherical harmonic rotations'));

for index=1:size(K,1)
  % It is this PLM2ROT call that we need for the functions also, see
  % e.g. LOCALIZATION, where it is applied to the eigenfunctions instead 
  %  Ktmp{index}=kindeks(plm2rot([dels dems Ktmp{index}],...
  %			      alp,bta,gam),3:4);
  % And then stick these guys back into the same format as they came in
  % See PLOTSLEP where I've done this before
  %  K(:,index)=Ktmp{index}(mzo);
  % Or all at once, really, which is slightly faster (but not much)
  K(:,index)=indeks(kindeks(plm2rot([dels dems Ktmp{index}],...
			     alp,bta,gam),3:4),mzo);
  waitbar(index/size(K,1),h)
end
close(h)
toc