
function hpiloc_smallHpiArray( n_coils, r_ring, meg_array, hpi_moment, noise, time, f_smpl , path )
%hpiloc_varyNcoil Summary of this function goes here
%   Detailed explanation goes here

%% -- Definitions --
load([path, 'grads.mat']);
load([path, 'headSurface.mat']);
load([path, 'exclude.mat']);
j=0;
include = [];
for i=1:102
    if ~any(i == exclude)      
        j = j+1;
        include(j) = i;
    end
end
hpigrad.tra = diag(ones(1,10));

if strcmp(meg_array,'elekta')
    meggrad = meggrad_elekta;
elseif strcmp(meg_array,'onscalp')
    meggrad = meggrad_onscalp;
else 
    disp('ERROR: wrong MEG array')
    meggrad = [];
end

headmodel = [];
headmodel.unit = 'm';
headmodel.type = 'infinite';

%hpicoilfreq = [ 370 420 470 520 570 620 670 720 770 820 ];
hpicoilfreq = [ 218 225 232 239 246 253 260 267 274 281 ];
hpi_freq = hpicoilfreq(1:length(hpigrad.label));

%% -- Calculations --
meg_error = zeros(102,1);

for k = 1:length(n_coils)     % Loop through noise levels
    for l = 1:length(include)-44
        m = include(l);
        %% HPI array
        hpi_array = [];
        tmp = sensorRing(r_ring, n_coils(k), meggrad_onscalp.coilpos(m,:), meggrad_onscalp.coilori(m,:));
        for i = 1:size(tmp,1)
            [hpi_array.pos(i,:), hpi_array.ori(i,:)] = snap2surf(tmp(i,:), bnd.pos, bnd.tri, bnd.nn);
        end
        hpigrad_array = [];
        hpigrad_array.unit = 'm';
        hpigrad_array.tra = zeros(n_coils(k));
        for i = 1:n_coils(k)
            hpigrad_array.label{1,i} = hpigrad.label{1,i};
            hpigrad_array.chanpos(i,:) = hpi_array.pos(i,:);
            hpigrad_array.coilpos(i,:) = hpi_array.pos(i,:);
            hpigrad_array.coilori(i,:) = hpi_array.ori(i,:);
            hpigrad_array.tra(i,i) = hpigrad.tra(i,i);
        end
   
       %% - Generate data -       
        amp = hpi_moment;
        t = (1:time)/f_smpl;
        data = [];
        data.time{1} = t;
        data.trial{1} = zeros(102,time);
        data.label = meggrad.label;
        for i = 1:n_coils(k)
            signal = amp*cos(2*pi*hpi_freq(i)*t);
            lf = ft_compute_leadfield(hpigrad_array.coilpos(i,:),meggrad,headmodel);
            lf = lf *  hpigrad_array.coilori(i,:)';
            data.trial{1} = data.trial{1}+lf*signal;
        end
        data.trial{1} = data.trial{1}+noise.*randn(102,time); % Add noise     

        %% - FFT -
        cfg = [];
        cfg.method = 'mtmfft';
        cfg.output = 'fourier';
        cfg.taper = 'boxcar'; 
        cfg.foi = hpi_freq(1:n_coils(k))';
        freq = ft_freqanalysis(cfg, data);

        %% - Freq reverse -
        freqrev = [];
        freqrev.label = hpigrad_array.label;
        freqrev.dimord = freq.dimord;
        freqrev.fourierspctrm = nan(1, n_coils(k), 102);
        freqrev.freq = (1:102)'; % all frequencies combined into one reverse HPI topography per channel
        freqrev.fourierspctrm = permute(freq.fourierspctrm(:,:,:), [1 3 2]);

        %% - Fit mags -
        tmp_coilpos = meggrad.coilpos;
        j = 0;
        for i=1:102
            if norm(meggrad.coilpos(i,:) - meggrad.coilpos(m,:)) < r_ring && j<5
                j= j+1;
                cfg = [];
                cfg.channel = 'all';
                cfg.headmodel = headmodel;
                cfg.gridsearch = 'no';
                cfg.dip.pos = tmp_coilpos(i,:);
                
                cfg.dip.unit = 'm';
                cfg.model = 'regional';
                cfg.frequency = i;
                cfg.grad = hpigrad_array; 
                megfit_individual = ft_dipolefitting(cfg, freqrev);
                meg_error(j)=norm(tmp_coilpos(i,:)-megfit_individual.dip.pos);
                clc;
                [k m]
            end
        end
        locError_mean(l,k) = mean(meg_error); 
        locError_max(l,k) = max(meg_error);
        locError_min(l,k) = min(meg_error);
        locError_std(l,k) = std(meg_error);
        locError_median(l,k) = median(meg_error); 
    end
end

%% - Temporarily save results -
res = [];
res.n_coils = n_coils;
res.r_ring = r_ring;
res.c_ring = meggrad.coilpos;
res.noise = noise;
res.meg_array = meg_array;
res.hpi_moment = hpi_moment;
res.locError_mean = locError_mean;
res.locError_max = locError_max;
res.locError_min = locError_min;
res.locError_std = locError_std;
res.locError_median = locError_median;

%% -- Save results -- 
save([path 'smallHpiArray_' meg_array '.mat'],'res');

%% -- Plot --
f = figure();
hold on
plot(n_coils,mean(locError_mean,1).*1e3,'*-')
hold off
xlabel('Number of HPI coils')
ylabel('Localization error [mm]')
title(['Localization error vs. number of HPI coils in small ring-array (r =' num2str(r_ring) ' m)'])
saveas(f, [path 'Loc_smallHpiArray/smallHpiArray_' meg_array '_mean'], 'fig');

f = figure();
hold on
plot(n_coils,mean(locError_median,1).*1e3,'*-')
hold off
xlabel('Number of HPI coils')
ylabel('Localization error [mm]')
title(['Localization error vs. number of HPI coils in small ring-array (r =' num2str(r_ring) ' m)'])
saveas(f, [path 'Loc_smallHpiArray/smallHpiArray_' meg_array '_median'], 'fig');
end

