Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
1 result

batch_processor_partitioned_rrho_tickonov.m

Blame
  • batch_processor_partitioned_rrho_tickonov.m 6.35 KiB
    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