function [XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, xlim, ylim)
% function [XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, xlim, ylim)
% 
% This function meshes the domain.
%
% INPUTS:
%       app:    The application data-structure. Should contain the fields
%                Nx and Ny for the number of x and y control volumes. 
%       xlim:   Array giving the limits of the domain in the x-direction.
%                xlim = [min(x), max(x)]
%       ylim:   Array giving the limits of the domain in the y-direction.
%                ylim = [min(y), max(y)]
%       
% OUTPUTS:
%       XP:     X-coordinate matrix for the Pressure/Tracer grid
%       YP:     Y-coordinate matrix for the Pressure/Tracer grid
%       XU:     X-coordinate matrix for the u-velocity grid
%       YU:     Y-coordinate matrix for the u-velocity grid
%       XV:     X-coordinate matrix for the v-velocity grid
%       YV:     Y-coordinate matrix for the v-velocity grid
%       dx:     Control volume size in the x-direction
%       dy:     Control volume size in the y-direction
%
% Written by:   Matt Ueckermann

dx = abs(xlim(2) - xlim(1))./(app.Nx);
dy = abs(ylim(2) - ylim(1))./(app.Ny); 

xp = linspace(xlim(1) + dx / 2, xlim(2) - dx / 2, app.Nx);
yp = fliplr(linspace(ylim(1) + dy / 2, ylim(2) - dy / 2, app.Ny));
xu = linspace(xlim(1) + dx, xlim(2) - dx, app.Nx - 1);
yu = fliplr(linspace(ylim(1) + dy / 2, ylim(2) - dy / 2, app.Ny));
xv = linspace(xlim(1) + dx / 2, xlim(2) - dx / 2, app.Nx);
yv = fliplr(linspace(ylim(1) + dy, ylim(2) - dy, app.Ny - 1));

[XP YP] = meshgrid(xp, yp);
[XU YU] = meshgrid(xu, yu);
[XV YV] = meshgrid(xv, yv);