%% Visualizing Geoid Height for Earth Geopotential Model 1996
% This example shows how to calculate the Earth's Geoid height using the EGM96 Geopotential 
% Model of the Aerospace Toolbox(TM) software. It also shows how to visualize 
% the results with contour maps overlaid on maps of the Earth.  
% The Mapping Toolbox(TM) and Simulink(R) 3D Animation(TM) are required to
% generate the visualizations.

% Copyright 2010-2020 The MathWorks, Inc.

%% Generating Values for Earth Geopotential Model 1996
% Calculate values for the Earth's Geopotential using the |geoidheight| 
% function to implement the EGM96 Geopotential Model.

%%
% The following code can be used to generate 260281 data points for
% calculating values of the Earth's Geoid height using |geoidheight|. To
% reduce computational overhead, this example includes a mat-file that
% contains this data.

% % Set amount of increment between degrees
% gridDegInc = 0.5;                   %degrees
% 
% % Longitude value in degrees to use for latitude sweep.
% lon = -180:gridDegInc:180;      %degrees
% 
% % Geodetic Latitude values to sweep.
% geod_lat = -90:gridDegInc:90;       %degrees
% 
% % Loop through longitude values for each array of latitudes -90:90.
% for lonIdx = size(lon,2):-1:1
%     
%     % Longitude must be the same dimension as the latitude array
%     lon = lon(lonIdx)*ones(1,numLatitude);      % degrees
%     geoidResults(1:end,lonIdx) = geoidheight(geod_lat,lon,'None');
%   
% end

%% Loading Geoid Data File and Coastal Data
geoidFileName = 'GeoidResults_05deg_180.mat';
load(geoidFileName);
coast = load('coastlines.mat');

%% Plot 2-D View of Geoid Height

% Create 2-D plot using |meshm|
h2D = figure;
set(h2D,'Position',[20 75 700 600],'Toolbar','figure');

% Reference matrix for mapping geoid heights to lat/lon on globe.
RRR = makerefmat('RasterSize',size(geoidResults), ...
    'Latlim', [-90 90], 'Lonlim', [-180 180] );

ast2DGeoidPlot(RRR,geoidResults,coast,gridDegInc)

% Viewing Geoid height using VR canvas
www2D = vrworld('astGeoidHeights.wrl');
open(www2D)

% Actual geoid heights for reference
geoidGrid = vrnode(www2D,'EGM96_Grid');
actualHeights = getfield(geoidGrid,'height'); %#ok<GFLD>

% Initialize heights to 0 for slider control
geoidGrid.height = 0*actualHeights;

% Size canvas for plotting and set parameters
geoidcanvas2D = vr.canvas(www2D,'Parent',h2D,...
    'Antialiasing', 'on','NavSpeed','veryslow',...
    'NavMode','Examine','Units', 'normalized',...
    'Viewpoint','Perspective','Position',[.15 .04 .7 .42]);

% Create slider
slid=astGeoidSlider(geoidcanvas2D);

%% Plot 3-D View of Geoid Height
 
h3D = figure;
set(h3D,'Position',[20 75 700 600]);

% Set up axes
hmapaxis = axesm ('globe','Grid', 'on');
set(hmapaxis,'Position',[.1 .5 .8 .4])

view(85,0)
axis off

% Plot data on 3-D globe
meshm(geoidResults,RRR)

% Plot land mass outline
plotm(coast.coastlat,coast.coastlon,'Color','k')
colormap('jet');

% Plot Title
title({'EGM96 Geoid Heights';['Grid Increment: ' ,num2str(gridDegInc), ' Degrees; Height Units: Meters']})

colorbar;

% 3-D Globe: Geoid Height Using VR Canvas
www3D = vrworld('astGeoidSphere.wrl');
open(www3D)

% Position canvas
geoidcanvas3D = vr.canvas(www3D,'Parent',h3D,...
    'Antialiasing', 'on','NavSpeed','veryslow',...
    'NavMode','Examine','Units', 'normalized',...
    'Position',[.15 .04 .7 .4]);
vrdrawnow;

%% Clean Up

close(h2D,h3D)
close(www2D);close(www3D);
delete(www2D);delete(www3D);

displayEndOfDemoMessage(mfilename)
