%% 2 steps, get data, make a plot
%% Step 1, get data
% TODO: Change this to your own path
f = '~/Documents/checkouts/baart_f/documents/presentations/probabilisticcourse/jb_ibn_srn_gtn_giy.nc';

% From file f, get variable height, starting at t=0, lat=0,lon=0, get 1
% timestep, and all lat and lon
h = nc_varget(f, 'height', [0,0,0], [1,-1,-1]);
% Get the whole time vector
t = nc_varget(f, 'time');
% Transform them into matlab dates
dates = datenum('1990-1-1') + double(t);
% Get the lats
lat = nc_varget(f, 'lat');
% Get the lons
lon = nc_varget(f, 'lon');
% Make them into two matrices (matlab doesn't get vectors for
% coordinates....)
[Lat, Lon] = meshgrid(lat,lon);

% Load some topography data (from matlab)
load topo

% Use a globe
grs80 = almanac('earth','grs80','km');
domeRadius =  3000;  % km
domeLat =  39;       % degrees
domeLon = -77;       % degrees
domeAlt = 0;         % km


%% Part 2 the plotting
% Use 3d rendering
figure('Renderer','opengl')
% Set the globe as an axis
ax = axesm('globe','Geoid',grs80,'Grid','on', ...
    'GLineWidth',1,'GLineStyle','-',...
    'Gcolor',[0.9 0.9 0.9],'Galtitude',100);
% Not sure what this does....
set(ax,'Position',[0 0 1 1]);
% Hide the axes.
axis equal off
% Look from the default 3d view
view(3)
% I want a black background, If I ask for black I get dark grey, at least
% it is an improvement....
colordef black

% Show the areas
geoshow('landareas.shp', 'FaceColor', 'black');
% Put the plot here on the lower left, 800x600px
set(gcf, 'Position', [0 0 800 600]) 
% Not sure what this is, but it shines...
material ([.7, .9, .8])
% Let's keep looking at the same spot so remember where we look at.
obj = camtarget;
% All these things are set to auto by default, which gives an animation
% which requires in epileptic risk warning.
campos('manual')
camva('manual')
camtarget('manual')
camzoom(1.5)

%% This is the animation part
nt = (length(t)-1);
% Let's scale it up so we can make it look more 2012 like....
scale = 4;
for i=0:nt
    % get data for the next timestep
    h = nc_varget(f, 'height', [i,0,0], [1,-1,-1]);
    % Plot a new surface
    handle = geoshow(Lat, Lon, h'*scale,'DisplayType','surface');
    % clim([-300,300]*scale);
    set(gca, 'Clim', [-300,300]*scale);
    % Let's show a colorbar because we scaled things up a bit
    cbh = colorbar('location', 'southoutside', 'Position' ,[0.2 0.1 0.6 0.0595]);
    % Downscale the labels
    set(cbh, 'XTickLabel', num2str(str2num(get(cbh, 'XTickLabel'))/scale))
    % Move around the earth a bit
    camorbit(5, 0)
    % Show some light from the right.
    camh = camlight('right');
    % This makes it shine more
    lighting phong; 
    % Show the date at the north Pole
    texth = text(0,0,8000,datestr(dates(i+1)), 'Color', 'white');
    % Store this frame
    F(i+1) = getframe(gcf);
    % Delete all the parts, for the next timestep
    delete(cbh);
    delete(handle);
    delete(texth);
    delete(camh);
    
end
%% Now play or store the movie.... 
% movie(F, 5)
% movie2avi(F,'mlanimation.avi', 'compression', 'None')