function batch_processor_partitioned_pos(Mday,nods,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state,timeKBR,obs_gnv,x0x)
% Batch processing algorithm for GRACE using positions observations. 
% seperation between local and global parameters, which are estimated for
% different arcs (eg. initial states daily and spherical harmonics 
% coefficients for the whole time. Using Partitioned Normal Equations
% based on:
% - Gunter's MSc(2000)thesis, page 27
% - Statistical Orbit Determination (Tapley et al., 2004), page 196-197
% 
% Input:    Mday = current arc number/index
%           FolderName = folder/path where intermediate results are saved
%           lmaxcs = order of spherical harmonic coefficients to be estimated
%           mKBR = number of observations in each arc
%           field = gravity field
%           data_plm = Legendre polynomial, see initplm
%           GM = GM of gravity field
%           ae = Earth radius of gravity field
%           lmaxf = parameters of background gravity field
%           state = initial states for arc [12x1]
%           timeKBR = times of observations
%           observation = vector of observations (eg. range rates) 
%           x0x = state deviation [12x1]
% 
% Output:   function has no output, output is saved in file for easy
%           paralleization
% 
% Exampel: see gfr_pos.m
%  
% Main functions
% - deriv: calculates all partials
% - hmat: makes observation-state mapping matrices Hxtilde and HcTilde 
% Other functions
% - vec2cs: reformatting spherical harmonics coefficients 
% - cs2sc: reformatting spherical harmonics coefficients
% - initplm: initialzes Legendre polynomial calculation
%
% Written by Florian Wöske, ZARM Uni Bremen, 2018-12.
% Adapted by Neda Darbeheshti, AEI, 2019-06.
%% set some initial parameters
% Minimum degree and order of coefficients to be estimated
lmincs=2;
noc=lmaxcs^2-lmincs^2+2*lmaxcs+1;%total number of coefficients to be estimated
% Load P0:A priori state variance???covariance matrix
P0x=eye(12);
P0c=eye(noc);

x0c = zeros(noc,1);

% numbers of observations for each time step
no = size(obs_gnv,2)-1; %(3 or 6 positions)

% All parameters mneeded in derivs function are put in a struct, making it
% easy to change or add the parameters of derivs 
deriv_params.field = field;
deriv_params.data_plm = data_plm;
deriv_params.GM = GM;
deriv_params.ae = ae;
deriv_params.lmaxf = lmaxf;
deriv_params.lmaxcs = lmaxcs;
deriv_params.noc = noc;

%% set integration parameters
% integration method (1: Runge-Kutta 4, 2: RK DP8, 3: ABM)
% RK4 needs 4 calls of deriv function, medium integration accuracy, DP8
% needs 13 calls of deriv, high accuracy also with bigger step size. ABM
% needs to be initialized with a single step integration method like RK or
% DP, as for all multistep integrators. After initialization it needs just
% 2 calls of DERIV and has a high accuracy. Multistep methods with high
% orders tend to get unstable when not decreasing step sizes.
int_meth = 3;
% id ABM integrator, order of ABM scheme (implemented: 1 to 12)
abm_order = 8;
lengthX0=12*(12+noc)+12;
% Define the state transition matrix:
phibuff = eye(12+noc);% a replacement for phi as a square matrix
phi = reshape(phibuff(1:12,1:end),12,(12+noc));%phi not to be a square matrix to save some space
% Form the extended state vector:
X0 = [state;reshape(phi,numel(phi),1)];
    % new
    X = X0;
    %Eqs. (6.3.22), (6.3.23), and (6.3.25)from Tapley's book
    Mxx=eye(12) / P0x;
    Mxc=zeros(12,noc);
    Mcc=eye(noc) / P0c;
    Nx=P0x \ x0x;
    Nc=P0c \ x0c;
    % alloc. vars.
    i_call = 0;
    dXdt=zeros(abm_order, lengthX0);
%% The if loop limits saving of'ygnv_save','Hx_save','Hc_save' just for last day
if Mday==nods
    ygnv_save = zeros(mKBR,no);
    Hx_save = zeros(mKBR,no,12);
    Hc_save = zeros(mKBR,no,noc);
    %Run through each observation:
    for ii = 1:mKBR
        % Read observation:
        Ygnv = obs_gnv(ii,2:end)';
        
        if ii > 1
            i_call = i_call+1; %count number of odeint evaluations
            [X, dXdt] = odeint(@deriv, X, dXdt, timeKBR(ii-1), 5, i_call, int_meth, abm_order, deriv_params);
        end
              
        % Extract the state transition matrix:
        phi = reshape(X(12+1:end), 12, 12+noc);
        [a,b] = size(phi);
        phibuff(1:a,1:b) = phi;
                 
        % Get the Htilde Matrix:
        [Htilde,Y_star_gnv]= hmat_pos(X,noc,no);
                       
        % Time update:
        H = Htilde*phibuff;
        Hx = H(:,1:12);
        Hc = H(:,13:end);
        
        % Calculate the observation deviation:
        y= Ygnv - Y_star_gnv;
        
        % accumulate Lambda:
        %L = L + H'*WMAT*H;
        Mxx=Mxx + Hx'*Hx;
        Mxc=Mxc + Hx'*Hc;
        Mcc=Mcc + Hc'*Hc;
        % accumulate N:
        %N = N + H'*WMAT*y;
        Nx = Nx + Hx'*y;
        Nc = Nc + Hc'*y;
        % Store data:
        ygnv_save(ii,:) = y;
        Hx_save(ii,:,:)=Hx;
        Hc_save(ii,:,:)=Hc;
            end
    invMxx=Mxx\eye(size(Mxx));
% save in a file
FileName= ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
File= fullfile(FolderName, FileName);
save(File,'invMxx','Mxc','Mcc','Nx','Nc','ygnv_save','Hx_save','Hc_save','-v7.3');
else
     %Run through each observation:
    for ii = 1:mKBR
        % Read observation:
          Ygnv = obs_gnv(ii,2:end)';
        
        if ii > 1
            i_call = i_call+1; %count number of odeint evaluations
            [X, dXdt] = odeint(@deriv, X, dXdt, timeKBR(ii-1), 5, i_call, int_meth, abm_order, deriv_params);
        end
             
        % Extract the state transition matrix:
        phi = reshape(X(12+1:end), 12, 12+noc);
        [a,b] = size(phi);
        phibuff(1:a,1:b) = phi;
                
        % Get the Htilde Matrix:
        [Htilde,Y_star_gnv]= hmat_pos(X,noc,no);
                   
        % Time update:
        H = Htilde*phibuff;
        Hx = H(:,1:12);
        Hc = H(:,13:end);
        
        % Calculate the observation deviation:
        y= Ygnv - Y_star_gnv;
        
        % accumulate Lambda:
        %L = L + H'*WMAT*H;
        Mxx=Mxx + Hx'*Hx;
        Mxc=Mxc + Hx'*Hc;
        Mcc=Mcc + Hc'*Hc;
        % accumulate N:
        %N = N + H'*WMAT*y;
        Nx = Nx + Hx'*y;
        Nc = Nc + Hc'*y;
              end
    invMxx=Mxx\eye(size(Mxx));
    % save in a file
FileName= ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
File= fullfile(FolderName, FileName);   
save(File,'invMxx','Mxc','Mcc','Nx','Nc');
end


