% Calculate sensitivities for fitting by finite difference 
%
% arFiniteDifferences(dp)    
%   dp:     parameter variation             [1e-3]


function arFiniteDifferences(dp)

global ar

if(isempty(ar))
    error('please initialize by arInit')
end

if(~exist('dp', 'var'))
    dp = 1e-3;
end

pRef = ar.p;
arChi2(false,[]);

ar.sresFD = ar.res' * ones(1,length(pRef)) + 0;

arWaitbar(0);
for jp=1:length(ar.pLabel)
    arWaitbar(jp, length(ar.pLabel));
    for jm=1:length(ar.model)
        for jc=1:length(ar.model(jm).condition)
            if(jp==1)
                ar.model(jm).condition(jc).suExpSimuFD = zeros(size(ar.model(jm).condition(jc).suExpSimu));
                ar.model(jm).condition(jc).sxExpSimuFD = zeros(size(ar.model(jm).condition(jc).sxExpSimu));
                ar.model(jm).condition(jc).szExpSimuFD = zeros(size(ar.model(jm).condition(jc).szExpSimu));
            end
            qp = ismember(ar.model(jm).condition(jc).p, ar.pLabel{jp}); %R2013a compatible
            if(sum(qp)>0)
                ar.model(jm).condition(jc).suExpSimuFD(:,:,qp) = ar.model(jm).condition(jc).uExpSimu;
                ar.model(jm).condition(jc).sxExpSimuFD(:,:,qp) = ar.model(jm).condition(jc).xExpSimu;
                ar.model(jm).condition(jc).szExpSimuFD(:,:,qp) = ar.model(jm).condition(jc).zExpSimu;
            end
        end
        for jd=1:length(ar.model(jm).data)
            if(jp==1)
                ar.model(jm).data(jd).syExpSimuFD = zeros(size(ar.model(jm).data(jd).syExpSimu));
                ar.model(jm).data(jd).systdExpSimuFD = zeros(size(ar.model(jm).data(jd).systdExpSimu));
                ar.model(jm).data(jd).sresFD = zeros(size(ar.model(jm).data(jd).sres));
                if(ar.config.fiterrors == 1)
                    ar.model(jm).data(jd).sreserrFD = zeros(size(ar.model(jm).data(jd).sres));
                end
            end
            qp = ismember(ar.model(jm).data(jd).p, ar.pLabel{jp}); %R2013a compatible
            if(sum(qp)>0)
                ar.model(jm).data(jd).syExpSimuFD(:,:,qp) = ar.model(jm).data(jd).yExpSimu;
                ar.model(jm).data(jd).systdExpSimuFD(:,:,qp) = ar.model(jm).data(jd).ystdExpSimu;
                ar.model(jm).data(jd).sresFD(:,:,qp) = ar.model(jm).data(jd).res;
                if(ar.config.fiterrors == 1)
                    ar.model(jm).data(jd).sreserrFD(:,:,qp) = ar.model(jm).data(jd).reserr;
                end
            end
        end
    end
    
    % disturb p(jp)
    ar.p(jp) = ar.p(jp) + dp;
    arChi2(false,[]);
    
    for jm=1:length(ar.model)
        for jc=1:length(ar.model(jm).condition)
            qp = ismember(ar.model(jm).condition(jc).p, ar.pLabel{jp}); %R2013a compatible
            if(sum(qp)>0)
                ar.model(jm).condition(jc).suExpSimuFD(:,:,qp) = ...
                    (ar.model(jm).condition(jc).uExpSimu - ar.model(jm).condition(jc).suExpSimuFD(:,:,qp)) / dp;
                ar.model(jm).condition(jc).sxExpSimuFD(:,:,qp) = ...
                    (ar.model(jm).condition(jc).xExpSimu - ar.model(jm).condition(jc).sxExpSimuFD(:,:,qp)) / dp;
                ar.model(jm).condition(jc).szExpSimuFD(:,:,qp) = ...
                    (ar.model(jm).condition(jc).zExpSimu - ar.model(jm).condition(jc).szExpSimuFD(:,:,qp)) / dp;
            end
        end
        for jd=1:length(ar.model(jm).data)
            qp = ismember(ar.model(jm).data(jd).p, ar.pLabel{jp}); %R2013a compatible
            if(sum(qp)>0)
                ar.model(jm).data(jd).syExpSimuFD(:,:,qp) = ...
                    (ar.model(jm).data(jd).yExpSimu - ar.model(jm).data(jd).syExpSimuFD(:,:,qp)) / dp;
                ar.model(jm).data(jd).systdExpSimuFD(:,:,qp) = ...
                    (ar.model(jm).data(jd).ystdExpSimu - ar.model(jm).data(jd).systdExpSimuFD(:,:,qp)) / dp;
                ar.model(jm).data(jd).sresFD(:,:,qp) = ...
                    (ar.model(jm).data(jd).res - ar.model(jm).data(jd).sresFD(:,:,qp)) / dp;
                if(ar.config.fiterrors == 1)
                    ar.model(jm).data(jd).sreserrFD(:,:,qp) = ...
                        (ar.model(jm).data(jd).reserr - ar.model(jm).data(jd).sreserrFD(:,:,qp)) / dp;
                end
            end
        end
    end
    
    ar.sresFD(:,jp) = (ar.res' - ar.sresFD(:,jp)) / dp + 0;
    
    % reset p
    ar.p = pRef;
    arChi2(false,[]);
end
arWaitbar(-1);