clear
close all

shipLat = 36.8020317;
shipLon = -121.7872017;
shipHeading = 10;
shipHeadingXY = -(shipHeading - 90);

load('mb15.mat');
SCCSMCAlat = [(36+51/60) (36+51/60) (36 + 48/60) (36 + 48/60) (36 + 51/60)];
SCCSMCAlon = [-(121+56/60) -(122+3.8/60) -(122+2.88/60) -(121 + 56/60) -(121 + 56/60)];
PLLSMCAlat = [(36+43/60) (36+43/60) (36 + 41/60) (36 + 41/60) (36 + 43/60)];
PLLSMCAlon = [-(121+56/60) -(122+1.3/60) -(122+0.8/60) -(121 + 56/60) -(121 + 56/60)];

zone = utmzone(shipLat, shipLon);
wgs84 = wgs84Ellipsoid;
ellipsoid = [wgs84.SemimajorAxis wgs84.Eccentricity];

axesm utm
setm(gca,'zone', zone, 'geoid', ellipsoid, 'LabelUnits', 'dm', ...
    'grid','on','Frame', 'on', 'FFill', 1000, 'FLineWidth', 0.5,  ...
    'meridianlabel','on','parallellabel','on', ...
    'MLineLocation', 0.016667, 'PLineLocation', 0.016667, ...
    'MLabelLocation', 0.083333, 'PLabelLocation', 0.083333);

[lat1, lon1] = minvtran(mb15.x(1), mb15.y(1));
[lat2, lon2] = minvtran(mb15.x(end), mb15.y(end));
setm(gca,'maplatlimit',[lat1 lat2], 'maplonlimit', [lon1 lon2]);

cm = haxby(256);
colormap(cm);
V = [0:20:100 200:100:4000];
V = -V;
[C, hF] = contourf(mb15.x,mb15.y,mb15.z, V);
set(hF, 'ShowText', 'on', 'LabelSpacing', 1e6);
fillm(SCCSMCAlat, SCCSMCAlon, 'r');
fillm(PLLSMCAlat, PLLSMCAlon, 'r');
hA = gca;
set(hA,'DataAspectRatio',[1 1 1]);