function geoworldmap (projection, fig, color)
% GEOWORLDMAP   Plot a global map using specified projection
%
% geoworldmap (projection, fig, color)
% Input:    projection  - Name of a supported projection
%                         Hammer, Sinusoidal, or Wagner
%           fig         - Number of figure window to use
%           color       - Color to use for the landmasses
%
% Example: geoworldmap ('Wagner', 3, 'r')

% Load in the landmasses we wish to plot
load Eurasiafrica.txt
load Americas.txt
load Australia.txt
load Greenland.txt
load Antarctica.txt
% Set the requested projection
geoinit (projection);
% Initialize the plot figure
figure (fig)
clf
% Paint the landmasses, one polygon at the time
geofill (Eurasiafrica(:,1), Eurasiafrica(:,2), color)
hold on
geofill (Americas(:,1), Americas(:,2), color)
geofill (Australia(:,1), Australia(:,2), color)
geofill (Greenland(:,1), Greenland(:,2), color)
geofill (Antarctica(:,1), Antarctica(:,2), color)
title (['The world according to the ' projection ' projection'])
grid on
axis equal

%%%%%%%%%%%% SUB FUNCTIONS FOR GEOWORLDMAP %%%%%%%%%%%%%%%%%%%

function geoinit (projection)
% GEOINIT   Sets the current map projection
%
% geoinit (projection)
% Input:    projection is the unique name of one of the supported
%           map projections such as 'Wagner', 'Sinusoidal', or 'Hammer'

global geoprojection georadius
georadius = 6371.008;   % Radius of spherical Earth in km

switch projection   % Assign the handle to the appropriate function
    case 'Sinusoidal'
        geoprojection = @geosinusoidal;
    case 'Wagner'
        geoprojection = @geowagner;
    case 'Hammer'
        geoprojection = @geohammer;
    otherwise
        disp (['Projection ' projection ' not implemented.  We default to sinusoidal'])
        geoprojection = @geosinusoidal;
end

function geofill (lon, lat, color)
% GEOPLOT   Fill a colored polygon on a map using current projection
% geofill (lon, lat, color)
% Input: lon and lat in degrees for the polygon
%        color to use

global geoprojection georadius

[x, y] = feval (geoprojection, lon, lat);
fill (x, y, color)

% HERE ARE THE ACTUAL MAP PROJECTION FUNCTIONS
% Maps are centered on Greenwich and all longitudes are assumed to be in
% the [-180,+180] range.

function [x, y] = geowagner (lon, lat)
% GEOWAGNER  Convert geographic positions to Wagner map coordinates
% [x, y] = geowagner (lon, lat)
% Input: lon, lat are geographic positions in degrees
% Output: x, y are the map positions using an Wagner projection

global georadius

S = 0.90631 .* sind (lat);
C0 = sqrt (1.0 - S.^2);
C1 = sqrt (2./ (1 + C0 .* cosd (lon ./ 3)));
x = 2.66723 .* georadius .* C0 .* C1 .* sind (lon ./ 3);
y = 1.24104 .* georadius .* S .* C1;


function [x, y] = geohammer (lon, lat)
% GEOHAMMER  geographic positions to Hammer map coordinates
% [x, y] = geohammer (lon, lat)
% Input: lon, lat are geographic positions in degrees
% Output: x, y are the map positions using a Hammer projection

global georadius
D = sqrt (2 ./ (1 + cosd (lat) .* cosd (lon ./ 2)));
x = 2 .* georadius .* D .* cosd (lat) .* sind (lon ./ 2);
y = georadius .* D .* sind (lat);

function [x, y] = geosinusoidal (lon, lat)
% GEOSINUSOIDAL  geographic positions to Sinusoidal map coordinates
% [x, y] = geosinusoidal (lon, lat)
% Input: lon, lat are geographic positions in degrees
% Output: x, y are the map positions using a Sinusoidal projection

global georadius
x = deg2rad (lon) .* georadius .* cosd (lat);
y = deg2rad (lat) .* georadius;

% Convenience trig functions that take degrees instead of radians

function y = sind (x)
y = sin (x .* pi ./ 180);

function y = cosd (x)
y = cos (x .* pi ./ 180);