Select Git revision
batch_processor_partitioned_pos.m
Axel Schnitger authored
m-files in the release directory are now executable from the release diretory. No need to change the path to the repo's root level.
batch_processor_partitioned_pos.m 6.35 KiB
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