function vals = spharm(n, m, Cnm, Snm, theta, phi);
%
%vals = function spharm(n, m, Cnm, Snm, theta, phi);
%
%Calculate spherical harmonic function of degree n, order m, with
%coefficients Cnm, Snm, at value(s) (theta, phi).
%    theta = colatitude, 0 at north pole, radians
%    phi   = longitude, radians
%
%This function uses the normalization that is standard in geodesy,
%following Kaula (1966). See also Lambeck (1988) Geophysical Geodesy.
%
%To calculate the value at a specific point, pass scalar values of theta
%and phi. In this case, val returns a scalar value.
%
%To calculate the value at a grid of points, pass vectors of theta and phi.
%In this case, val returns a matrix of values at the grid defined by the
%vectors theta and phi. The indices for val(i,j) are:
%      i = index on values of theta
%      j = index on values of phi
%
%Author: Jeff Freymueller, Univ. of Alaska, Jan. 2004

% Calculate the Pnm. Note MATLAB's normalization is for the associated
% Legendre functions only: (Pnm, Pnm) = 1. When multiplied by the
% cos(m*phi) and sin(m*phi) terms, MATLAB's version does not result in a
% fully normalized spherical harmonic. Thus we add a renormalization to
% give the standard geodetic usage. It seems easiest to calculate the
% Schmidt normalization and renormalize that. 
%
% Renormalize the function in accord with Kaula standard used in geodesy
% This is the difference between Kaula's fully normalized and the Schmidt
% normalization.
%
% See for example, Hofmann-Wellenhof and Moritz, Physical Geodesy, 2nd ed.,
%  page 23 (section 1.10)

Pnm = legendre(n, cos(theta), 'sch');
if ( m == 0 )
    renorm = sqrt(2*(2*n+1));
else
    renorm = sqrt(n+1/2);
end

% Compute the values of the harmonic at all "gridpoints" as defined above.
%   Note that this expression is conveniently vectorized by MATLAB.

vals = renorm*Pnm(m+1,:)'*(Cnm*cos(m*phi) + Snm*sin(m*phi));
