% Example program to call spharm.m
%
% This program makes two plots. The first is a rectangular plot of the
% selected harmonic (red=positive, blue=negative). The second plot shows
% the function mapped onto a sphere.
%
% Note that in this example, Cnm and Snm are now matrices. The (n,m)
% harmonic coeffficients are stored in Cnm(n,m+1) and Snm(n,m+1) because
% MATLAB can't handle an index of zero. We specify C0 (for P00) as a
% separate variable.
%
%Author: Jeff Freymueller, Univ. of Alaska, Jan. 2004

% Use the n=3, m=2 harmonic
%  C3,2 = 1 ; all Snm = 0

C0 = 0;
Cnm = zeros(3,4);
Snm = zeros(3,4);

n = 3;
m = 2;
Cnm(n,m+1) = 1;

% Uncomment this to add another harmonic
%Snm(2,2) = 1;

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

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

% Calculate the harmonic function a grid of points defined by the theta and
% phi values

vals = spharm_to_degree(n, C0, Cnm, Snm, theta, phi);

% 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, vals);
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(vals);
for i=1:nx,
    vals2(i,:) = vals(nx+1-i,:);
end

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