function [OUT,mStatesInn,mStatesNOKF,sodaPar] = sodaEnKF(sodaPar,parVec)

mConstants = evalin('base','mConstants');
mStatesNOKFInit = evalin('base','mStatesNOKFInit');


for k=1:numel(parVec)
    mPars{k,1} = sodaPar.parNames{k};
    mPars{k,2} = parVec(k);
end
clear k

measOperator = 1;

nIterObs = numel(sodaPar.obsPoints{1,2});
nStatesKF = sodaPar.nStatesKF;
nStatesNOKF = sodaPar.nStatesNOKF;
statesViz = sodaPar.statesViz;
nStatesViz = numel(statesViz);
nMembers = sodaPar.nEnsembleMembers;
% nOutputs = sodaPar.nOutputs;
% nOptPars = sodaPar.nOptPars;

if strncmp(sodaPar.modeStr,'soda',4)
    TMP=strcmp(sodaPar.stochForcePars,'optimized');
    stochForceParsFlag = any(TMP(:,2));

    if stochForceParsFlag
        Ix = find(TMP(:,2));
        sodaPar.stochForcePars{Ix,2} = parVec(sodaPar.stochForceParsIx);
    end
    clear TMP
end

% mOutputs = cell([nIterObs,nOutputs,nMembers]);

defaultNaNArray = repmat(NaN,[nIterObs,nStatesKF,nMembers]);

mStatesPrior = defaultNaNArray;
mStatesPert = defaultNaNArray;
mStatesPost = defaultNaNArray;
mStatesInn = defaultNaNArray;

mStatesNOKF = repmat(NaN,[nIterObs,nStatesNOKF,nMembers]);

iIter = 1;

for iMember=1:nMembers
    mStatesNOKF(iIter,1:nStatesNOKF,iMember)=mStatesNOKFInit;
end


if strncmp(sodaPar.modeStr,'soda',4)
    stochForce = defaultNaNArray;
    eval(sodaPar.errModelCallStr)
    kalmanGain = repmat(NaN,[nIterObs,nStatesKF,nStatesKF]);
end
perturbedObs = defaultNaNArray;

if ~isfield('sodaPar','isMeas')
    nAux = size(mConstants,1);
    sodaPar.isMeas = repmat(false,[nAux,1]);
    for k=1:nAux
        if any(strcmp(mConstants{k,1},sodaPar.measNames))
            sodaPar.isMeas(k)=true;
        end
    end
end




switch sodaPar.initState
    case 'unifrand'
        mStatesPost(1,1:nStates,1:nMembers) = sodaUnifRandDraw(sodaPar,'stateSpace');
    case 'reference'
        for iMember=1:nMembers
        
            if strncmp(sodaPar.modeStr,'soda',4)
                
                mStatesPrior(1,1:nStatesKF,iMember) = sodaPar.initStateRef;
                
                together = [mStatesPrior(1,1:nStatesKF,iMember);sodaPar.stateSpaceLoBound];
                mStatesPrior(1,1:nStatesKF,iMember) = max(together,[],1);

                together = [mStatesPrior(1,1:nStatesKF,iMember);sodaPar.stateSpaceHiBound];
                mStatesPrior(1,1:nStatesKF,iMember) = min(together,[],1);

                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %             

                mStatesPost(1,1:nStatesKF,iMember) =...
                sodaPar.initStateRef + randn(1,nStatesKF)*sodaPar.initStateErr;

                together = [mStatesPost(1,1:nStatesKF,iMember);sodaPar.stateSpaceLoBound];
                mStatesPost(1,1:nStatesKF,iMember) = max(together,[],1);

                together = [mStatesPost(1,1:nStatesKF,iMember);sodaPar.stateSpaceHiBound];
                mStatesPost(1,1:nStatesKF,iMember) = min(together,[],1);
                
            else

                mStatesPrior(1,1:nStatesKF,iMember) = sodaPar.initStateRef;
                mStatesPost = mStatesPrior;
                
            end
        end
end

for iIter = 1:nIterObs-1
    
    % Propagate all ensemble members one time step...
    % further using the specified model: 
    [mStatesPrior,mStatesNOKF] = sodaPropStates(sodaPar,iIter,...
       mConstants,mStatesPrior,mStatesPost,mPars,mStatesNOKF);

    switch sodaPar.modeStr
        case 'scemua'
            
            mStatesPost = mStatesPrior;
            mStatesPert = mStatesPrior;
            
        case 'soda'
            % determine the stochastic forcing value:
            eval(sodaPar.errModelCallStr)

            % Add a model error to the simulated state values:
            mStatesPert = sodaPerturbModel(sodaPar,mStatesPrior,stochForce,mStatesPert,iIter);

            % Compute the ensemble covariance matrix, Pe (See...
            % Equation 47 of Evensen et al. 2002):
            ensembleCov = sodaCalcEnsembleCov(mStatesPert,iIter);

            % For all ensemble members, perturb the current observation...
            % using a randn distribution:
            [perturbedObs,sodaPar] = sodaPerturbObs(sodaPar,mConstants,perturbedObs,iIter);

            % Compute the measurement error covariance matrix, Re (See...
            % Equation 51 of Evensen et al. 2002):
            measErrCov = sodaCalcMeasErrCov(sodaPar,perturbedObs,iIter);


            % Caclulate the classical Kalman Gain (See...
            % Equation 52 of Evensen et al. 2002):
            kalmanGain = sodaCalcKalmanGain(ensembleCov,measOperator,measErrCov,kalmanGain,iIter);


            % Calculate the state innovations:
            mStatesInn = sodaCalcInnovations(sodaPar,kalmanGain,perturbedObs,...
                measOperator,mStatesPert,mStatesInn,iIter);

            % Add innovation to ensemble of states (See Equation 54 of ...
            % Evensen et al. 2002):
            mStatesPost = sodaUpdateStates(sodaPar,mStatesPert,...
                mStatesPost,mStatesInn,iIter);

        case 'reset'

            [perturbedObs,sodaPar] = sodaPerturbObs(sodaPar,mConstants,perturbedObs,iIter);            
            mStatesPert(iIter+1,:) = mStatesPrior(iIter+1,:);
            
            mStatesInn(iIter+1,:) = perturbedObs(iIter+1,:) - mStatesPrior(iIter+1,:);
            mStatesPost(iIter+1,:) = perturbedObs(iIter+1,:);
            
    end
end

meanEnPredPrior = measOperator*mean(mStatesPrior,3);
meanEnPredPost = measOperator*mean(mStatesPost,3);
OUT = measOperator*mean(mStatesPert,3);

if sodaPar.plotEnsemble && any(strcmp(sodaPar.modeStr,{'soda','reset'}))

    trueObs = mConstants{sodaPar.isMeas,2};

    figure(873)
    clf
    set(gcf,'numbertitle','off','name','ensemble results')

    for k = [1:nStatesViz,nStatesViz+1]
        
        
        if k==nStatesViz+1
            iStateViz = statesViz(end);
        else
            iStateViz = statesViz(k);
        end
        
        if nStatesViz>1
            subplot(ceil(sqrt(nStatesViz+1)),ceil(sqrt(nStatesViz+1)),k)
        else
            subplot(1,2,k)
        end
        
%         for iMember=1:nMembers
%             xOffset = (1:nIterObs-1)-0.125+(((iMember-1)/(nMembers-1))*0.25);
%             plotHandles(:,iMember) = plot(xOffset,mStatesPrior(1:nIterObs-1,iStateViz,iMember)','o',...
%                 xOffset,mStatesPert(1:nIterObs-1,iStateViz,iMember)','+',...
%                 xOffset,mStatesPost(1:nIterObs-1,iStateViz,iMember)','s',...
%                 xOffset,perturbedObs(1:nIterObs-1,iStateViz,iMember)','d');
%             hold on
%         end
%         plotHandles2 = plot(1:nIterObs,trueObs(:,iStateViz),'k*',...
%              1:nIterObs-1,meanEnPredPrior(1:nIterObs-1,iStateViz),...
%              1:nIterObs-1,meanEnPredPost(1:nIterObs-1,iStateViz));

        markerSizeArray = [ 6 , 6 ; 6 , 6 ; 6       ,6];
        colorArray=       {'k','c';[1,0.5,0],'m';[7,244,1]/255,[0,0.5,0]};
        markerArray =     {'s','o';'s','o';'s','o'};
        

        hres1 = 0.2*min(sodaPar.obsPoints{1,2}(2:end)-sodaPar.obsPoints{1,2}(1:end-1))/2;
        hres2 = hres1/(nMembers*2);
        for iMember=1:nMembers
            xOffset = -0.5*hres1+((iMember-1)/nMembers)*hres1+hres2;
            plotHandles{1,1}(:,iMember) = plot(sodaPar.obsPoints{1,2}+xOffset,mStatesPrior(:,iStateViz,iMember)','ks');
            hold on
            plotHandles{2,1}(:,iMember) = plot(sodaPar.obsPoints{1,2}+xOffset,mStatesPert(:,iStateViz,iMember)','ms');
            plotHandles{3,1}(:,iMember) = plot(sodaPar.obsPoints{1,2}+xOffset,mStatesPost(:,iStateViz,iMember)','ys');
            plotHandles{2,2}(:,iMember) = plot(sodaPar.obsPoints{1,2}+xOffset,perturbedObs(:,iStateViz,iMember)','mo');
            
        end

        plotHandles{1,2} = plot(sodaPar.obsPoints{1,2},trueObs(:,iStateViz),'ko');
        
        plotHandles{4,1} = plot(sodaPar.obsPoints{1,2},meanEnPredPrior(:,iStateViz),'-k',...
            'marker','none','linewidth',2);
        plotHandles{5,1} = plot(sodaPar.obsPoints{1,2},meanEnPredPost(:,iStateViz),'-',...
            'color',colorArray{3,1},'marker','none','linewidth',2);

%         set(plotHandles,'markersize',4)
        
        
        for r=1:size(colorArray,1)
            for c=1:size(colorArray,2)
                set(plotHandles{r,c}(:),'markerfacecolor',colorArray{r,c},...
                                        'markeredgecolor','k',...
                                        'marker',markerArray{r,c},...
                                        'markersize',markerSizeArray(r,c))
            end
        end
        clear r c


%         set(plotHandles2(2),'color',[0,0,1])
%         set(plotHandles2(3),'color',[0,0.5,0])
        set(gca,'ylim',[sodaPar.stateSpaceLoBound(iStateViz),...
            sodaPar.stateSpaceHiBound(iStateViz)],'xlim',[1,nIterObs])

%         xlabel('observations')
        ylabel(sodaPar.stateNamesKF{iStateViz})
        
        if k==nStatesViz+1
            axis tight
            set(gca,'xlim',min(get(gca,'xlim'))+[-2,-1])
%             hLegend = legend([plotHandles(1:4,1);plotHandles2],...
%                 '\eta(\Psi_e^a,X,\theta)',...
%                 '\Psi_e^f',...
%                 '\Psi_e^a',...
%                 'H(\Psi_t)+\epsilon_t',...
%                 'true obs',...
%                 'mean(\Psi_e^f)',...
%                 'mean(\Psi_e^a)');
%             set(hLegend,'interpreter','tex');
            hLegend = legend([plotHandles{1,1}(1),plotHandles{2,1}(1),...
                plotHandles{3,1}(1),plotHandles{1,2}(1),plotHandles{2,2}(1),...
                plotHandles{4,1}(1),plotHandles{5,1}(1)],...
                'model, prior',...
                'model, perturbed',...
                'model, posterior',...
                'observation, true',...
                'observation, perturbed',...
                ['mean(',char(39),'model, prior',char(39),')'],...
                ['mean(',char(39),'model, posterior',char(39),')']);
            set(hLegend,'interpreter','none');

            axis off
            
        end
    end

    
    for iSubplot=1:nStatesViz
        if nStatesViz>1
            subplot(ceil(sqrt(nStatesViz+1)),ceil(sqrt(nStatesViz+1)),iSubplot)
        else
            subplot(1,2,iSubplot)
        end
        set(gca,'xlim',[sodaPar.vizWindow])
    end
    
    
    drawnow
end

% (perturbedObs<mStatesPost & mStatesPost<mStatesPert) |...
% (mStatesPert<mStatesPost & mStatesPost<perturbedObs)


if strncmp(sodaPar.modeStr,'soda',4) && stochForceParsFlag
    
    sodaPar.stochForcePars{Ix,2}='optimized';
    
end
