
function beaglemap(select,area)

if nargin <2 
  area='indian-ocean';
end

% set map limits

area=lower(area);
switch area
   case 'indian-ocean'
     latlim=[-45 0];
     lonlim=[44.99 120];
     etopomag=2;
   case 'mauritius'
     latlim=[-25 -15];
     lonlim=[55 65];
     etopomag=1;
   otherwise
     disp('Unknown area.')
end

% load K18 track

k18=load('k18');

% set waypoints

%            Perth                    Mauritius
plat = [  -32.0475  -19.9510  -20.1467 -20.1524 ];
plon = [  115.7442   57.5791   57.4756  57.4888 ];

% Open figure and initialize axis

figure;
worldmap(latlim,lonlim)

% select map type

select=lower(select);
switch select

  case {'geoid','hrgeoid','lrgeoid','egm96'}

    % Display EGM96 geoid (1x1 degree bins)
    switch select
      case 'lrgeoid'
        geoid=load('geoid');
      otherwise
        [geoid.geoid,geoid.geoidrefvec]=egm96geoid(1,latlim,lonlim);
        geoid.geoidlegend=geoid.geoidrefvec;
    end
    geoshow(geoid.geoid, geoid.geoidrefvec, 'DisplayType', 'texturemap');
    hcb = colorbar('horiz');
    set(get(hcb,'Xlabel'),'String','EGM96 Geoid Heights in Meters.')

    % Mask out all the land.
    geoshow('landareas.shp', 'FaceColor', 'black');
    %coast=load('coast');
    %geoshow(coast.lat, coast.long)

    % Compute the elevation profile, defaulting the track type to great circle and the interpolation type to bilinear:
    [z,rng,lat,lon] = mapprofile(geoid.geoid,geoid.geoidlegend,plat,plon);

    % Draw the great-circle track
    geoshow(lat, lon)

    % Draw the K18 track
    geoshow(k18.lat, k18.lon,'DisplayType','Point','Color','Black')

    % Plot the geoid height along the track

    figure; 
    plot(deg2km(rng),z,'r')
    title('EGM Geoid Heights Perth-Mauritius.')
    ylabel('Geoid height [m]')
    xlabel('Distance along track [km]')

  case {'topo','etopo','hrtopo','lrtopo'}

    % Display topography and bathymetry
    switch select
      case 'lrtopo'
        topo=load('topo');
      case 'hrtopo'
        [topo.topo, topo.topolegend]=etopo('etopo5.dat',1,latlim,lonlim);
      otherwise
        [topo.topo, topo.topolegend]=etopo('etopo5.dat',etopomag,latlim,lonlim);
    end
    geoshow(topo.topo, topo.topolegend, 'DisplayType', 'texturemap');
    demcmap(topo.topo)
    hcb = colorbar('horiz');
    set(get(hcb,'Xlabel'),'String','Topography and Bathymetry in Meters.')

    % Compute the elevation profile, defaulting the track type to great circle and the interpolation type to bilinear:
    [z,rng,lat,lon] = mapprofile(topo.topo,topo.topolegend,plat,plon);

    % Draw the great-circle track
    geoshow(lat, lon,'Color','Black')

    % Plot the topography along the track

    figure; 
    plot(deg2km(rng),z,'r')
    title('ETOPO5 Ocean Bathymetry Perth-Mauritius')
    ylabel('Depth [m]')
    xlabel('Distance along track [km]')

  otherwise
    disp('Unknown method.')
end

return