function [phi1a,phi1b,phi2,phi3]=eigbas(m,n,x);
% finds vector eigenfunctions of K on sphere for n,m as function of x=cos(theta).
% phi1a multiplies (a/r)^n.
% phi1b multiplies (a/r)^(n+2)
% phi2 multiplies (a/r)^(n+2)
% phi3 multiplies (a/r)^(n+1)

ind=1:length(x);
ind2=ind(abs(1-x.^2)>1e-6); %handle x=1 and x=-1 separately
indp=ind((abs(1-x.^2)<=1e-6)&(x>0));
indm=ind((abs(1-x.^2)<=1e-6)&(x<0));


    P=legendre(n,x);
    eta=sqrt(4*pi*gamma(n+m+1)/gamma(n-m+1)/(2*n+1));
    Ynm=(-1)^m*P(m+1,:);
    dYnm=zeros(size(Ynm));
    pYnm=zeros(size(Ynm));
    if(m==0)
        if(n==0) ; 
        else
            dYnm=P(2,:);
        end
    else
        dYnm(ind2)=-(-1)^m*(-m*x(ind2).*P(m+1,ind2)./sqrt(1-x(ind2).^2));
        if(m<n)
            dYnm(ind2)=dYnm(ind2)+(-1)^m*P(m+2,ind2);
        end
        pYnm(ind2)=Ynm(ind2)./sqrt(1-x(ind2).^2);
        if(m==1)
            dYnm(indp)=dYnm(indp)+n*(n+1)/2*ones(size(indp));
            dYnm(indm)=dYnm(indm)-(-1)^(n+1)*n*(n+1)/2*ones(size(indm));
            pYnm(indp)=dYnm(indp);
            pYnm(indm)=-dYnm(indm);
        end
        
    end
    Anm=[n*Ynm; dYnm; sqrt(-1)*m*pYnm];
    Bnm=[-(n+1)*Ynm; dYnm; sqrt(-1)*m*pYnm];
    Cnm=[zeros(size(Ynm)); sqrt(-1)*m*pYnm; -dYnm];
    
phi1a=sqrt(n*(2*n+1))/eta*(1/n/(2*n+1)*Anm-(2*n-1)/(2*n+1)/2/(n+1)*Bnm);
phi1b=sqrt(n*(2*n+1))/eta*(2*n-1)/(2*n+1)/2/(n+1)*Bnm;
phi2=1/eta/sqrt((n+1)*(2*n+1))*Bnm;
phi3=sqrt(-1)/eta/sqrt(n*(n+1))*Cnm;