function [Harm,GradHarm] = HarmsAtPt(MaxDeg,Colat,Long)
%HarmsAtPt Values  of Harmonic Functions and their Gradients.
%   Y = HarmsAtPt(MaxDeg,Colat,Long) returns the column of values at (Colat,Long)
%   for lexicographically ordered Harmonic Functions 
%   with subscript < MaxDeg
%   
%   [Y,V] = HarmsAtPt(...) also returns the 3D gradient vector
%   Arguments are in degree measure

%   author: Andrey Khokhlov. 
%   $Revision: 0.25 $  $Date: 2011/02/05 08:15:33 $

if Colat <= 0.0 || Colat >= 180.0 
    error('MATLAB:HarmsAtPt:InvalidInput', 'Input argument Colat must be between 0 and 180.')
end

NSz=(MaxDeg+1)^2-1;
Harm=zeros(NSz,1);
if nargout == 2
  GradHarm=ones(NSz,3);
end
NormCff=1.;
TOGRAD=180/pi;
TORAD=1./TOGRAD;
   Sin_Theta=sin(Colat*TORAD);
   Cos_Theta=cos(Colat*TORAD); 
   [LL,DL]=LegDerLeg(MaxDeg,Cos_Theta);
   ilex=1;
   for subsc=1: MaxDeg
       subsc_shift=(subsc*(subsc+1))/2;
       for suprs=0:subsc
           PPind = suprs+subsc_shift+1;
           Sin_mm_Phi=sin(suprs*(Long*TORAD));
           Cos_mm_Phi=cos(suprs*(Long*TORAD)); 
           if nargout == 2
               GradHarm(ilex,1)= -Cos_mm_Phi*Sin_Theta*DL(PPind );
               GradHarm(ilex,2)=  (suprs*LL(PPind )*Sin_mm_Phi)/Sin_Theta;
               GradHarm(ilex,3)= -(1 + subsc)*Cos_mm_Phi*LL(PPind );
           end
           Harm(ilex)=Cos_mm_Phi*LL(PPind );
           ilex=ilex+1;
           if suprs == 0  
               continue;
           end
           if nargout == 2
               GradHarm(ilex,1)= -Sin_mm_Phi*Sin_Theta*DL(PPind );
               GradHarm(ilex,2)= -(suprs*LL(PPind )*Cos_mm_Phi)/Sin_Theta;
               GradHarm(ilex,3)= -(1 + subsc)*Sin_mm_Phi*LL(PPind );
           end
           Harm(ilex)=Sin_mm_Phi*LL(PPind );
           ilex=ilex+1;           
       end
   end
   
