You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: spm_opm_plotScalpData.m
    path of file: ./repos/spm/toolbox/MEEGtools
    
    the code of the file until where you have to start completion: function [f] = spm_opm_plotScalpData(S)
% Display M/EEG interpolated sensor data on a scalp image
% FORMAT D = spm_opm_amm(S)
%   S               - input structure
%  fields of S:
%   S.D             - SPM MEEG object                                - Default: no Default
%   S.T             - time point to initalise to                    - Default: first sample  
%   S.display       - string to deermine what is plotted   -Default: 'RADIAL'
% OUTPUT:
%   f          - the handle of the figure which displays the interpolated
%                data
%__________________________________________________________________________
%
% This function creates a figure whose purpose is to display an
% interpolation of the sensor data on the scalp (as an image).
%__________________________________________________________________________

% Jean Daunizeau
% Copyrigh
    
function [C,R,out] = fitSphere(x,y,z)
% fitSphere  Fit sphere.
%       A = fitSphere(x,y,z) returns the parameters of the best-fit
%       [C,R,out] = fitSphere(x,y,z) returns the center and radius
%       sphere to data points in vectors (x,y,z) using Taubin's method.
% IN:
%   - x/y/z: 3D cartesian coordinates
% OUT:
%   - C: the center of sphere coordinates
%   - R: the radius of the sphere
%   - out: an output structure devoted to graphical display of the best fit
%   sphere

% Make sugary one and zero vectors
l = ones(length(x),1);
O = zeros(length(x),1);

% Make design mx
D = [(x.*x + y.*y + z.*z) x y z l];

Dx = [2*x l O O O];
Dy = [2*y O l O O];
Dz = [2*z O O l O];

% Create scatter matrices
M = D'*D;
N = Dx'*Dx + Dy'*Dy + Dz'*Dz;

% Extract eigensystem
[v, evalues] = eig(M);
evalues = diag(evalues);
Mrank = sum(evalues > eps*5*norm(M));

if (Mrank == 5)
    % Full rank -- min ev corresponds to solution
    Minverse = v'*diag(1./evalues)*v;
    [v,evalues] = eig(inv(M)*N);
    [dmin,dminindex] = max(diag(evalues));
    pvec = v(:,dminindex(1))';
else
    % Rank deficient -- just extract nullspace of M
    pvec = null(M)';
    [m,n] = size(pvec);
    if m > 1
        pvec = pvec(1,:);
    end