
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 'wa'
     latlim=[-35 -25];
     lonlim=[105 120];
     etopomag=1;
   case 'mauritius'
     latlim=[-25 -15];
     lonlim=[55 65];
     etopomag=1;
   otherwise
     disp('Unknown area.')
end

% load K18 track

k18=load('k18');

% load Beagle track

beagle=[ 76420    220   -32.0515  115.7433   -27.53    0.34 ; ...
         37942    177   -31.8650  115.5128   -27.66    0.63 ; ...
         18068    131   -30.2733  113.0898   -31.34    0.75 ; ...
         30782    138   -30.0742  112.8548   -31.39    0.77 ; ...
         36856    188   -29.9296  112.6448   -31.39    0.81 ; ...
         53397     37   -29.4986  111.9224   -32.75    0.80 ];

% 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')

    % Draw the beagletrack

    % Compute the elevation profile, defaulting the track type to great circle and the interpolation type to bilinear:
    [z2,rng2,lat2,lon2] = mapprofile(geoid.geoid,geoid.geoidlegend,beagle(:,3),beagle(:,4))

    % Draw the great-circle track
    geoshow(lat2, lon2,'Color','red')
    geoshow(beagle(:,3), beagle(:,4),'DisplayType','Point','Color','red')

    % Plot the geoid height along the great circle track

    value = ltln2val(geoid.geoid, geoid.geoidrefvec,beagle(:,3),beagle(:,4))
    beagle(:,5)-5.55
    beagle(:,5)-5.55-value
    
    rng3(1)=0;
    for ii=2:size(beagle,1)
      [rngii, az] = distance(beagle(ii-1,3),beagle(ii-1,4),beagle(ii,3),beagle(ii,4));
       rng3(ii)=rng3(ii-1)+rngii;
    end

    err=3*beagle(:,6)./sqrt(beagle(:,2));
    
    figure; 
    plot(deg2km(rng2),z2,'r')
    hold on
    %plot(deg2km(rng3),value,'r+')
    plot(deg2km(rng3),beagle(:,5)-5.55,'r*')
    errorbar(deg2km(rng3),beagle(:,5)-5.55,-err,err,'r*')
    a=axis;
    a(1)=-10;a(4)=a(4)+0.5;
    axis(a);
    title('Geoid Heights Perth-Mauritius.')
    ylabel('Geoid height [m]')
    xlabel('Distance along track [km]')
    legend('EGM96','Beagle II')
    % Plot the geoid height along the beagle 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')

    % Compute the elevation profile, defaulting the track type to great circle and the interpolation type to bilinear:
    [z2,rng2,lat2,lon2] = mapprofile(topo.topo,topo.topolegend,beagle(:,3),beagle(:,4))

    % Draw the great-circle track
    geoshow(lat2, lon2,'Color','red')
    geoshow(beagle(:,3), beagle(:,4),'DisplayType','Point','Color','red')

    % 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