function headmodel = ft_headmodel_singleshell(mesh, varargin)

% FT_HEADMODEL_SINGLESHELL creates a volume conduction model of the

% head for MEG based on a realistic shaped surface of the inside of

% the skull.

%

% The method implemented in this function allows for a simple and

% fast method for the MEG forward calculation for one shell of arbitrary

% shape, based on a correction of the lead field for a spherical

% volume conductor by a superposition of basis functions, gradients

% of harmonic functions constructed from spherical harmonics.

%

% This function implements

% G. Nolte, "The magnetic lead field theorem in the quasi-static

% approximation and its use for magnetoencephalography forward calculation

% in realistic volume conductors", Phys Med Biol. 2003 Nov 21;48(22):3637-52.

%

% Use as

% headmodel = ft_headmodel_singleshell(mesh, ...)

%

% Optional input arguments should be specified in key-value pairs and can include

% order = number of iterations in series expansion (default = 10)

%

% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

% Copyright (C) 2012-2018, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL

%

% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org

% for the documentation and details.

%

% FieldTrip is free software: you can redistribute it and/or modify

% it under the terms of the GNU General Public License as published by

% the Free Software Foundation, either version 3 of the License, or

% (at your option) any later version.

%

% FieldTrip 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 General Public License for more details.

%

% You should have received a copy of the GNU General Public License

% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.

%

% $Id$

% order of spherical spherical harmonics

order = ft_getopt(varargin, 'order', 10);

if isnumeric(mesh) && size(mesh,2)==3

% assume that it is a Nx3 array with vertices

% convert it to a structure, this is needed to determine the units further down

mesh = struct('pos', mesh);

elseif isstruct(mesh) && isfield(mesh,'bnd')

% take the triangulated surfaces from the input structure

mesh = mesh.bnd;

end

% replace pnt with pos

mesh = fixpos(mesh);

if ~isstruct(mesh) || numel(mesh)>1 || ~isfield(mesh, 'pos')

ft_error('the input mesh should be a set of points or a single triangulated surface')

end

% represent the mesh in a headmodel strucure

% the computational parameters will be added later on by ft_prepare_vol_sens

headmodel = [];

headmodel.bnd = mesh;

headmodel.type = 'singleshell';

headmodel.order = order;

headmodel = ft_determine_units(headmodel);