function hpiloc_varyPos( pos_error_range, meg_array, hpi_moment, noise, time, f_smpl ,n_reps, path )
%hpiloc_varyPos Summary of this function goes here
%   Detailed explanation goes here

%% -- Definitions --
load([path, 'grads.mat']);
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));

locError_mean = zeros(n_reps,length(pos_error_range));
locError_median = zeros(n_reps,length(pos_error_range));

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


%% - 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:length(hpi_freq)
    signal = amp*cos(2*pi*hpi_freq(i)*t);
    lf = ft_compute_leadfield(hpigrad.coilpos(i,:),meggrad,headmodel);
    lf = lf *  hpigrad.coilori(i,:)';
    data.trial{1} = data.trial{1}+lf*signal;
end
data.trial{1} = data.trial{1}+noise.*randn(102,time); % Add noise

for k = 1:length(pos_error_range)     % Loop through noise levels
    for m = 1:n_reps
        pos_error = pos_error_range(k)*2*(rand(length(hpi_freq),3)-1); % +/- pos_error_range
        hpigrad_false = hpigrad;
        for i = 1:length(hpi_freq)
            hpigrad_false.coilpos(i,:) = hpigrad.coilpos(i,:)+pos_error(i,:);
            hpigrad_false.chanpos(i,:) = hpigrad.coilpos(i,:)+pos_error(i,:);
        end
        %% - Plot FFT -
        % cfg = [];
        % cfg.method = 'mtmfft';
        % cfg.output = 'pow';
        % cfg.taper = 'hanning';
        % freq = ft_freqanalysis(cfg, data);
        % cfg = [];
        % cfg.operation = 'log10';
        % cfg.parameter = 'powspctrm';
        % freq = ft_math(cfg, freq);
        % figure
        % plot(freq.freq, freq.powspctrm);

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

        %% - Freq reverse -
        freqrev = [];
        freqrev.label = hpigrad.label;
        freqrev.dimord = freq.dimord;
        freqrev.fourierspctrm = nan(1, 10, 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;
        parfor i=1:102
            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_false; 
            megfit_individual = ft_dipolefitting(cfg, freqrev);
            meg_error(i)=norm(tmp_coilpos(i,:)-megfit_individual.dip.pos);
            clc;
            [k m]
        end
        locError_mean(m,k) = mean(meg_error);
        locError_max(m,k) = max(meg_error);
        locError_min(m,k) = min(meg_error);
        locError_std(m,k) = std(meg_error);
        locError_median(m,k) = median(meg_error);

        %% - Plot Mag fit -
        % figure
        % for i=1:102
        %   ft_plot_dipole(megfit_individual{i}.dip.pos, megfit_individual{i}.dip.mom, 'color', 'r', 'unit', 'm');
        % end
        % ft_plot_sens(meggrad, 'coil', true, 'coilorientation', true, 'coildiameter', 0.01);        
    end
end

%% - Temporarily save results -
res = [];
res.pos_error_range = pos_error_range;
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;
res.n_reps = n_reps;

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

%% -- Plot --
f = figure();
hold on
plot(pos_error_range.*1e3,mean(locError_mean,1).*1e3,'*-')
hold off
xlabel('Position error [+/- mm]')
ylabel('Localization error [mm]')
title('Localization error vs. HPI position error')
saveas(f, [path 'Loc_pos/pos_' meg_array '_mean'], 'fig');

f = figure();
hold on
plot(pos_error_range.*1e3,mean(locError_median,1).*1e3,'*-')
hold off
xlabel('Position error [+/- mm]')
ylabel('Localization error [mm]')
title('Localization error vs. HPI position error')
saveas(f, [path 'Loc_pos/pos_' meg_array '_median'], 'fig');

end

