% Example program to calculate geoid heights using spherical harmonic
% coefficients from EGM96 to to degree and order nmax. It reads the coefficients
% from the EGM96 file.
%
% The geoid height values refer to the WGS-84 ellipsoid. In order to refer
% heights to the ellipsoid, it must first subtract the even-numbered zonal
% coefficients that go into the definition of the ellipsoid. See W. Torge
% (1980) Geodesy, section 5.2 for example, or page 20 of Lambeck (1988).
% For the WGS-84 ellipsoid, we remove the zonal coefficients J2 through J10
% based on the WGS-84 reference values.
%
% See http://cddisa.gsfc.nasa.gov/926/egm96/doc/S11.HTML for specific
% values of J2-J10. The values are just below equation 11-4.3. Note that there
% is an error in equation 11-4.3. The sqrt(2n+1) should be sqrt(4n+1).
%
% The equation for the geoid height anomaly is 
%
%    N = R * SUM(n=2,inf) [ SUM(m=0,n) [ Pnm(cos(theta)*
%                { Cnm*cos(m*phi) + Snm*(sin(m*phi) } ] ]
%
%  where the Cnm, Snm have had the zonal correction made, theta is
%  colatitude and phi is longitude.
%
% This program ignores the zero degree term, which amounts to a correction of
% -53 cm for EGM96 and WGS-84.
%
% This program makes two plots. The first is a rectangular plot of the
% potential (red=positive, blue=negative). The second plot shows
% the potential mapped onto a sphere.
%
%Author: Jeff Freymueller, Univ. of Alaska, Jan. 2004

% Set the maximum degree and order here:

nmax = 30;

% Read EGM96. We read the nominal values of GM, a (earth radius), and the
% spherical harmonic coefficients Cnm and Snm. Each of these are
% dimensioned Cnm(360,361) for EGM96, which is complete to degree and order
% 360.

[GM, a, Cnm, Snm] = readEGM96;

% Subtract the even zonals that are part of the ellipsoid definition for
% WGS-84

Cnm(2,1)  = Cnm(2,1)  + 0.108262982131*10^-2/sqrt(5);      % J2
Cnm(4,1)  = Cnm(4,1)  - 0.237091120053*10^-5/sqrt(9);      % J4
Cnm(6,1)  = Cnm(6,1)  + 0.608346498882*10^-8/sqrt(13);     % J6
Cnm(8,1)  = Cnm(8,1)  - 0.142681087920*10^-10/sqrt(17);    % J8
Cnm(10,1) = Cnm(10,1) + 0.121439275882*10^-13/sqrt(21);    % J10

% Define a regular grid of theta (colatitude) and phi (longitude)

theta = linspace(0, pi, 101);
phi = linspace(0, 2*pi, 201);

% 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. 

% Initialize vals with zeros. The summation of terms for geoid anomalies
% begins with n=2. The n=0 term (GM/r in the potential) is part of the
% definition of the ellipsoid, so it is subtracted implicitly along with
% the even zonals above.

N = zeros(length(theta), length(phi));

% Do the summation over degree and order
%    Note that if nmax<2, MATLAB will not execute this loop. Some languages
%    execute a loop once no matter what, so beware if this is translated.

for in = 2:nmax,

    Pnm = legendre(in, cos(theta), 'sch');
    
    for m = 0:in,

%       Renormalize the function in accord with Kaula standard used in geodesy
%       This is the difference between Kaula's fully normalized and the Schmidt
%       normalization

        if ( m == 0 )
            renorm = sqrt(2*(2*in+1));
        else
            renorm = sqrt(in+1/2);
        end

%       Add the contribution from this n, m
        
        N = N + renorm*a*Pnm(m+1,:)'*( Cnm(in,m+1)*cos(m*phi) ...
                                         + Snm(in,m+1)*sin(m*phi) );
    end
end

% Print out the range of geoid heights

disp('Minimum and maximum geoid height')
disp( [max(max(N)) min(min(N)) ])

% Now make figures

% first figure: a rectangular plot
%  Note the transformation of the axes.
%  x-axis: longitude given in degrees
%  y-axis: latitude in degrees
%          recall that latitude = 90 - colatitude

figure
pcolor(phi*180/pi, 90-theta*180/pi, N);
shading('flat');

% second figure: map the harmonic onto a sphere.
%   We need to reverse the order of the rows, so that the northern
%   hemisphere (theta = 0 to pi/2) ends up on the top of the sphere.

[nx ny] = size(N);
for i=1:nx,
    N2(i,:) = N(nx+1-i,:);
end

figure
sphere; h = findobj(gcf, 'Type', 'surface');
set(h, 'CData', N2, 'FaceColor', 'texturemap')
axis equal;
ylabel('long = 0')
