inoutgradvecsdwcap.m

function [V,C]=inoutgradvecsdwcap(TH,Lin,Lout,m)
% [V,C]=inoutgradvecsdwcap(TH,Lin,Lout,m)
%
% INPUT:
%
% TH         Angular extent of the spherical cap, in degrees
%            OR: Angles of two spherical caps and we want the ring between 
%            them [TH1 TH2] 
% Lin        Maximum spherical-harmonic degree for inner sources
%            or passband [Lintmin Lintmax]
% Lout       Maximum spherical-harmonic degree for outer sources
%
% OUTPUT:
%
% C     Matrix containing vector Slepian function coefficients for the Elm
%       and Flm vector spherical harmonics IN ADDMOUT FORMAT
% V     Eigenvalues (conditioning values)
%
% See also gradvecsdwcap, inoutgradvecglmalpha
%
% Last modified by plattner-at-alumni.ethz.ch, 5/5/2017

if length(Lout)>1
    error('No Bandpass for external field')
end


maxLin=max(Lin);
lp=length(Lin)==1;
bp=length(Lin)==2;

if bp
    Lin=[min(Lin) max(Lin)];
end

defval('nth',0)%,720)
defval('vcut',-1)%,eps*10)
defval('grd',1)
defval('method','gl')

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

if(m>max(maxLin,Lout))
  error('Order cannot exceed degree')
end


% Filename of saved data
dirname=fullfile(getenv('IFILES'),'INOUTGRADVECSDWCAP');
if length(TH)==1
    if lp    
        fnpl=fullfile(dirname,sprintf(...
          'INOUTGRADVECSDWCAP-%g-%i-%i-%i.mat',...
          TH,Lin,Lout,m));
    elseif bp
      fnpl=fullfile(dirname,sprintf(...
          'INOUTGRADVECSDWCAP-%g-%i_%i-%i-%i.mat',...
          TH,Lin(1),Lin(2),Lout,m));
    else
      error('The degree range should be either one or two numbers')
    end 
elseif length(TH)==2
    if lp
        fnpl=fullfile(dirname,sprintf(...
          'INOUTGRADVECSDW-%g-%g-%i-%i-%i.mat',...
          max(TH),min(TH),Lin,Lout,m));
    elseif bp
        fnpl=fullfile(dirname,sprintf(...
          'INOUTGRADVECSDW-%g-%g-%i_%i-%i-%i.mat',...
          max(TH),min(TH),Lin(1),Lin(2),Lout,m));
    else
      error('The degree range should be either one or two numbers')
    end 
else
    error('Bad choice for TH')
end    
    

if exist(fnpl,'file')==2   %&& (vcut>0) & 1==3
  load(fnpl)
  disp(sprintf('%s loaded by inoutgradvecsdwcap.m',fnpl))
else
    % First load Pm and Bm for L=max(Lin,Lout)
    L=max(maxLin,Lout);    
    if length(TH)==2
        %Pm=kernelpm(max(TH),L,m)-kernelpm(min(TH),L,m);
        [~,~,~,~,~,~,~,~,~,~,Pm1]=sdwcap(max(TH),L,m,0,-1);
        [~,~,~,~,~,~,~,~,~,~,Pm2]=sdwcap(min(TH),L,m,0,-1);        
        Pm=Pm1-Pm2;
        clear Pm1
        clear Pm2
        Bm=kernelbm(max(TH),L,m)-kernelbm(min(TH),L,m);                
    else       
        %Pm=kernelpm(TH,L,m);
        [~,~,~,~,~,~,~,~,~,~,Pm]=sdwcap(max(TH),L,m,0,-1);
        Bm=kernelbm(TH,L,m);
    end    
             
    % Now take different matrix subsets for EE, FF, and EF  
    % If m>Lin or m>Lout, then some of these matrices will simply be empty. 
    
    % For EE
    PmEE=Pm(1:maxLin-m+1,1:maxLin-m+1);
    if m==0
        % In this case we need to add the m=0 zero rows and columns
        BmEE=zeros(size(PmEE));
        BmEE(2:end,2:end)=Bm(1:maxLin,1:maxLin);         
    else
        BmEE=Bm(1:maxLin-m+1,1:maxLin-m+1);
    end

    % For FF
    if m==0
        PmFF=Pm(2:Lout+1,2:Lout+1);
        BmFF=Bm(1:Lout,1:Lout);
    else
        PmFF=Pm(1:Lout-m+1,1:Lout-m+1);
        BmFF=Bm(1:Lout-m+1,1:Lout-m+1);
    end
    
    
    % For EF
    if m==0
        PmEF=Pm(1:maxLin+1,2:Lout+1);
        BmEF=zeros(maxLin+1,Lout);
        BmEF(2:end,:)=Bm(1:maxLin,1:Lout);
    else
        PmEF=Pm(1:maxLin-m+1,1:Lout-m+1);
        BmEF=Bm(1:maxLin-m+1,1:Lout-m+1);
    end        
    

    
    % Now calculate factor matrices for P and B to calculate E and F
    theLE=(m:maxLin)';
    theLF=(max(m,1):Lout)';
    
    facPE= sqrt( (theLE+1)./(2*theLE+1) );
    facPF= sqrt( (theLF  )./(2*theLF+1) );
    facBE=-sqrt( (theLE  )./(2*theLE+1) );
    facBF= sqrt( (theLF+1)./(2*theLF+1) );
    
    facPEmat=spdiags(facPE,0,length(facPE),length(facPE));
    facPFmat=spdiags(facPF,0,length(facPF),length(facPF));
    facBEmat=spdiags(facBE,0,length(facBE),length(facBE));
    facBFmat=spdiags(facBF,0,length(facBF),length(facBF));
    
    KEE=  facPEmat*PmEE*facPEmat + facBEmat*BmEE*facBEmat;
    KFF=  facPFmat*PmFF*facPFmat + facBFmat*BmFF*facBFmat;
    KEF=  facPEmat*PmEF*facPFmat + facBEmat*BmEF*facBFmat;
    
    
    % Plattner 5/11/2016: Now delete the low degrees
    if length(Lin)==2 & m<min(Lin);
        KEE=KEE(min(Lin)-m+1:end,min(Lin)-m+1:end);
        KEF=KEF(min(Lin)-m+1:end,:);
    end    
    
    % Now assemble the entire matrix
    K = [KEE KEF;KEF' KFF];           
    
    % And make sure it is symmetric
    try
        fprintf('Numerical asymmetry %g\n',norm(K-K'));
    catch
        fprintf('Numerical asymmetry %g\n',norm(full(K-K')));
    end
    K=(K+K')/2;
    
    % Calculate eigenvalues and eigenvectors; C'*C=I
    [C,V]=eig(full(K));
    
    [V,isrt]=sort(sum(V,1),'descend');
    
    C=C(:,isrt);
    
    % Plattner 5/11/2016: Now add zeros for the low degrees
    if length(Lin)==2 & m<min(Lin);
        C=[zeros(min(Lin)-m,size(C,2));C];
    end
    
    save(fnpl,'C','V');
    
end