function batch_processor_partitioned(Mday,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state,timeKBR,observation,x0x)

% BATCH_PROCESSOR_PARTITIONED provides a batch processing algorithm for GRACE
% range-rate observations. 
%
% Separation between local and global parameters, which are estimated for
% different arcs. E.g. 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
%           parallelization
% 
% Example: see gfr_parallel.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: initializes Legendre polynomial calculation
%
% Written by Neda Darbeheshti, AEI Hannover, 2018-07
% Last updated by Florian Wöske, ZARM Uni Bremen,2018-07


%% 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);

% states, rho and rhodot 
states = zeros(mKBR,12);
GiR_save = zeros(mKBR,2);

rhodot_deviation = zeros(mKBR,1);
Hx_save = zeros(mKBR,12);
Hc_save = zeros(mKBR,noc);

% All parameters needed 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;
%L = eye(12+noc) / P0;
%N = P0 \ 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);

% Run through each observation:
for ii = 1:mKBR
    % Read observation:
    Y = observation(ii,:)';
    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
    
    % save the states
    states(ii,:) = X(1:12)';
    
    % 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]= hmat(X,noc);
    % save rho and rhodot 
    GiR_save(ii,:) = Y_star;
    
    % Time update:
    H = Htilde*phibuff;
    Hx = H(1:12);
    Hc = H(13:end);
    % Calculate the observation deviation:
    y = Y - Y_star(2,:);
    % 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:
    rhodot_deviation(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,'results/OutputArcParall', FileName);
% save(File,'Mxc','Nx','Mcc','Nc','iL','iN','invMxx','rhodot_deviation','Hx_save','Hc_save');
save(File,'invMxx','Mxc','Mcc','Nx','Nc','rhodot_deviation','Hx_save','Hc_save', 'states', 'GiR_save');
