function batch_processor_partitioned_rrho_tickonov(Mday,nods,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state,timeKBR,obs_kbr,x0x)
% Batch processing algorithm for GRACE using range-rate observations with Tickonov regularization.
% 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_rrho_tickonov.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);
% 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;
%new Tickonov regularisation by Florian
x_d = 1e-6;
v_d = 0.1;
Mxx = diag([x_d x_d x_d v_d v_d v_d x_d x_d x_d v_d v_d v_d]);
Mxc = zeros(12,noc);
Mcc = zeros(noc);
Nx = zeros(12,1);
Nc = zeros(noc,1);
% alloc. vars.
i_call = 0;
dXdt=zeros(abm_order, lengthX0);
%% The if loop limits saving of'ykbr_save','Hx_save','Hc_save' just for last day
if Mday==nods
ykbr_save = zeros(mKBR,1);
Hx_save = zeros(mKBR,12);
Hc_save = zeros(mKBR,noc);
%Run through each observation:
for ii = 1:mKBR
% Read observation:
Ykbr = obs_kbr(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
% 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_kbr]= hmat(X,noc);
% Time update:
H = Htilde*phibuff;
Hx = H(:,1:12);
Hc = H(:,13:end);
% Calculate the observation deviation:
y= Ykbr - Y_star_kbr(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:
ykbr_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','ykbr_save','Hx_save','Hc_save','-v7.3');
else
%Run through each observation:
for ii = 1:mKBR
% Read observation:
Ykbr = obs_kbr(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
% 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_kbr]= hmat(X,noc);
% Time update:
H = Htilde*phibuff;
Hx = H(:,1:12);
Hc = H(:,13:end);
% Calculate the observation deviation:
y = Ykbr - Y_star_kbr(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;
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