function kgrid = makeGrid(varargin)
%MAKEGRID   Create k-space grid structure.
%
% DESCRIPTION:
%       makeGrid creates a MATLAB structure containing the grid coordinates
%       and wavenumber matrices for use in k-space simulations and
%       reconstructions. Note, the grid structures are indexed as: (x) in
%       1D; (z, x) in 2D; (z, x, y) in 3D.
%
% USAGE:
%       kgrid = makeGrid(Nx, dx)
%       kgrid = makeGrid(Nx, dx, Nz, dz)
%       kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz)
%
% INPUTS:
%       Nx, Ny, Nz  - number of pixels in each Cartesian direction
%       dx, dy, dz  - pixel size in each Cartesian direction [m]
%
% OUTPUTS:
%       kgrid       - k-space grid structure used by the simulation and
%                     reconstructions functions within k-Wave
%
%       The structure kgrid has the following fields:
%
%       kgrid.k       - ND grid of the scalar wavenumber
%       kgrid.k_max   - maximum spatial frequency supported by the grid
%       kgrid.t_array - evenly spaced array of time values (set to 'auto')
%
%       And for each spatial dimension x, y, z:
%
%       kgrid.Nx    - number of pixels 
%       kgrid.dx    - size of an individual pixel [m]
%       kgrid.x     - plaid ND grid of the x coordinates centered about 0 [m]  
%       kgrid.x_off - plaid ND grid of the x coordinates beginning at 0 [m]
%       kgrid.x_size- length of grid dimension [m]
%       kgrid.kx    - plaid ND grid of the wavenumber components centered
%                     about 0
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 12th March 2009
%       last update - 4th December 2009
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010 Bradley Treeby and Ben Cox
%
% See also cart2grid, interpCartData, kspaceFirstOrder1D,
% kspaceFirstOrder2D, kspaceFirstOrder3D, kspaceSecondOrder, ndgrid,
% makeTime, smooth

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% assign the input values to the grid object
if nargin == 6
    kgrid.Ny = varargin{3};
    kgrid.dy = varargin{4};
    kgrid.Nz = varargin{5};
    kgrid.dz = varargin{6};
elseif nargin == 4
    kgrid.Nz = varargin{3};
    kgrid.dz = varargin{4};
elseif nargin ~= 2
    error('Incorrect number of input arguments');
end
kgrid.Nx = varargin{1};
kgrid.dx = varargin{2};

switch nargin
    case 2
        % assign the grid parameters for the x spatial direction
        [nx, kgrid.kx, kgrid.x_size] = makeDim(kgrid.Nx, kgrid.dx);
        
        % define a spatial grid that is centered about 0
        kgrid.x = nx*kgrid.x_size; 

        % define a spatial grid that begins at 0
        kgrid.x_off = (nx + 0.5)*kgrid.x_size;  

        % define the scalar wavenumber based on the wavenumber components
        kgrid.k = abs(kgrid.kx);
        
        % define maximum supported frequency
        kgrid.k_max = max(abs(kgrid.kx(:)));
    case 4
        % assign the grid parameters for the x and z spatial directions
        [nx, kx_vec, kgrid.x_size] = makeDim(kgrid.Nx, kgrid.dx);
        [nz, kz_vec, kgrid.z_size] = makeDim(kgrid.Nz, kgrid.dz);
        
        % define a spatial grid that is centered about 0
        [kgrid.z, kgrid.x] = ndgrid(nz*kgrid.z_size, nx*kgrid.x_size); 

        % define a spatial grid that begins at 0
        [kgrid.z_off, kgrid.x_off] = ndgrid((nz + 0.5)*kgrid.z_size, (nx + 0.5)*kgrid.x_size);  

        % define plaid grids of the wavenumber components centered about 0
        [kgrid.kz, kgrid.kx] = ndgrid(kz_vec, kx_vec);

        % define the scalar wavenumber based on the wavenumber components
        kgrid.k = sqrt(kgrid.kx.^2 + kgrid.kz.^2);
        
        % define maximum supported frequency
        kgrid.k_max = min([max(abs(kgrid.kz(:))), max(abs(kgrid.kx(:)))]);        
    case 6
        % assign the grid parameters for the x ,y and z spatial directions
        [nx, kx_vec, kgrid.x_size] = makeDim(kgrid.Nx, kgrid.dx);
        [nz, kz_vec, kgrid.z_size] = makeDim(kgrid.Nz, kgrid.dz);
        [ny, ky_vec, kgrid.y_size] = makeDim(kgrid.Ny, kgrid.dy);

        % define a spatial grid that is centered about 0
        [kgrid.z, kgrid.x, kgrid.y] = ndgrid(nz*kgrid.z_size, nx*kgrid.x_size, ny*kgrid.y_size);

        % define a spatial grid that begins at 0
        [kgrid.z_off, kgrid.x_off, kgrid.y_off] = ndgrid((nz + 0.5)*kgrid.z_size, (nx + 0.5)*kgrid.x_size, (ny + 0.5)*kgrid.y_size);

        % define plaid grids of the wavenumber components centered about 0
        [kgrid.kz, kgrid.kx, kgrid.ky] = ndgrid(kz_vec, kx_vec, ky_vec);

        % define the scalar wavenumber based on the wavenumber components
        kgrid.k = sqrt(kgrid.kx.^2 + kgrid.ky.^2 + kgrid.kz.^2);
        
        % define maximum supported frequency
        kgrid.k_max = min([max(abs(kgrid.kz(:))), max(abs(kgrid.kx(:))), max(abs(kgrid.ky(:)))]);         
end

% set t_array to 'auto' by default
kgrid.t_array = 'auto';

% subfunction to create the grid parameters for a single spatial direction
function [nx, kx_vec, x_size] = makeDim(Nx, dx)

% define the discretisation of the spatial dimension such that there is
% always a DC component
if rem(Nx, 2) == 0
    % grid dimension has an even number of points
    nx = (-0.5:1/Nx:0.5-1/Nx).'; %HACK PUT BACK
%     nx = (-0.5:1/(Nx-1):0.5).'; %DELETE THIS BIT
else
    % grid dimension has an odd number of points
    nx = (-0.5:1/(Nx-1):0.5).';%HACK PUT BACK
%     nx = (-0.5:1/Nx:0.5-1/Nx).';%DELETE THIS BIT
end

% assign the size parameter
x_size = dx*Nx;

% define the wavenumber vector components
kx_vec = (2*pi/dx).*nx;