MATLAB Animation Code

% By Jesper G. Pedersen

% Netcdf file containing simulated data:
ncfile = 'tapes.nc';

% Read simulated values of qc(x,y,z) at time step 13:
qc = ncread(ncfile,'qc',[1 1 1 13],[inf inf inf 1]);

% Number of grid points in x-,y-,z-directions:
[nx,ny,nz]=size(qc);

% Grid resolution:
dx = 10;
dy = 10;
dz = 5;

% Grid points:
x = 0:dx:dx*(nx-1);
y = 0:dy:dy*(ny-1);
z = 0:dz:dz*(nz-1);

% Create meshgrid-matrices X and Y of size (nx,ny) for interpolation:
[Y,X] = meshgrid(y,x);
% X contains x-values increasing row-wise and
% Y contains y-values increasing column-wise.

% Simulated qc is interpolated to points separated by 5 m in the horizontal
% directions to match the vertical resolution and make the animation look nice:
xi=x(1):5:x(end);
yi=y(1):5:y(end);
zi=z;

[Yi,Xi] = meshgrid(yi,xi);

qci = zeros(length(xi),length(yi),length(zi));
for k=1:nz
    qci(:,:,k)=interp2(Y,X,qc(:,:,k),Yi,Xi,'spline');    
end
clear qc

% Create custom colormap composed of "polarmap" (from the Matlab Central File Exchange)
% and the built-in "gray" colormap:
cmap1 = polarmap(101);      %Dark blue to red
cmap1 = cmap1(15:51,:);     %Light blue to white
cmap2 = flipud(gray(90));   %White to black
cmap2 = cmap2(1:60,:);      %White to gray
cmapqc = [cmap1;cmap2];     %Light blue to white to gray


% Plot interpolated qc as a function of xi and zi at each yi-value (except at
% yi(end) where qci = qci(:,1,:) due to periodic boundary conditions):

% All text and numbers are displayed using the latex interpreter:
in = 'interpreter';
la = 'latex';

%Font sizes:
fs1 = 13;
fs2 = 11;

% Loop through yi-indices:
for j=1:length(yi)-1
    
    f=figure('Visible','Off');
    
    %Plot qci':    
    imagesc(xi/1000,zi/1000,squeeze(1000*qci(:,j,:))')
    colormap(cmapqc)
    % qci(:,j,:) is transposed to have x (rows) varying along the horizontal axis,
    % and z (columns) along the vertical axis.
    
    % imagesc sets the direction of the vertical axis to 'reverse'.
    % Here it is set back to 'normal' and the two axis are set 'equal':
    set(gca,'ydir','normal')
    axis equal
    
    % Insert x-, and y-labels and set axis limits and ticklabel font size:
    xlabel('x [km]',in,la,'fontsize',fs1)
    ylabel('z [km]',in,la,'fontsize',fs1)
    axis([0 xi(end)/1000 0 1])
    set(gca,'fontsize',fs2)
    
    % Range of qci-values to be represented by the colormap:
    caxis([0 0.9])
    
    % Title displaying the y-position using 4 digits including 3 decimals:
    title(['$y=' num2str(yi(j)/1000,'%4.3f\n') '$~km'],in,la,'fontsize',fs1);
    
    % Insert colorbar with label and specify ticklabel font size:
    cb=colorbar;
    ylabel(cb,'$q_c \: \mathrm{[g\,kg^{-1}]}$',in,la,'fontsize',fs1)
    set(cb,'fontsize',fs2)
    
    % Set the size of the figure in inches:
    sx = 11;
    sy = 4;
    set(f, 'PaperpositionMode','manual')
    set(f, 'PaperUnits', 'inches');
    set(f, 'Papersize',[sx sy])
    set(f, 'Paperposition',[0. 0.0 sx sy])
    
    % Figure is named according to y-position.
    if j < 10
        fname=['qc' '_y00' num2str(j)];
    elseif j<100
        fname=['qc' '_y0' num2str(j)];
    else
        fname=['qc' '_y' num2str(j)];
    end
        
    % Figure is saved as an eps file with 100 dots per inch:
    print(f,'-depsc','tempfile.eps','-r100')
    
    % Font of ticklabels is set to match the font used elsewhere in the figure:
    fin = fopen('tempfile.eps');
    fout = fopen([fname '.eps'],'w');
    while ~feof(fin)
        s = fgetl(fin);
        s = strrep(s,'mwa_cmr10','Helvetica');
        fprintf(fout,'%s\n',s);      
    end
    fclose(fin);
    fclose(fout);
    system('rm tempfile.eps');
end

% When all .eps files have been created, they are made into pdf-files (with epstopdf), and
% an infinitely looping gif animation is created with imagemagick.
%
% In a Linux terminal:
%  "find . -name "*.eps" -exec epstopdf {} \;"
%  and
%  "convert -delay 5 -density 100 -loop 0 *.pdf qc.gif"