From 51ae64264854c2769602c95a2949f27e0fd7050c Mon Sep 17 00:00:00 2001
From: Axel Schnitger <axel.schnitger@aei.mpg.de>
Date: Fri, 13 Sep 2019 12:58:32 +0200
Subject: [PATCH] Add proper path handling.
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.
---
.gitignore | 2 +
README.md | 141 ++++++-
function_eba/calculatePnm.m | 41 +++
function_eba/designmatrixn.m | 41 +++
function_eba/potential_EBAs.m | 41 +++
function_eba/sortCoefficients.m | 38 ++
function_eba/threeDmulti_rot.m | 45 +++
.../batch_processor_partitioned_pos.m | 184 ++++++++++
.../batch_processor_partitioned_pos_rrho.m | 201 ++++++++++
.../batch_processor_partitioned_rrho.m | 179 +++++++++
...atch_processor_partitioned_rrho_tickonov.m | 182 +++++++++
function_gfr/hmat_pos.m | 42 +++
function_gfr/save_vars_for_plotting.m | 17 +
get_input_data.sh | 31 +-
.../continue_gfr_par_from_iteration.m | 32 +-
gfr_parallel.m => release1/gfr_parallel.m | 61 +--
.../save_vars2continue_itr.m | 4 +-
release2/gfr_eba.m | 151 ++++++++
release2/gfr_pos.m | 333 +++++++++++++++++
release2/gfr_pos_rrho.m | 347 ++++++++++++++++++
release2/gfr_rrho.m | 317 ++++++++++++++++
release2/gfr_rrho_tickonov.m | 322 ++++++++++++++++
release2/results/.gitkeep | 0
23 files changed, 2674 insertions(+), 78 deletions(-)
create mode 100644 function_eba/calculatePnm.m
create mode 100644 function_eba/designmatrixn.m
create mode 100644 function_eba/potential_EBAs.m
create mode 100644 function_eba/sortCoefficients.m
create mode 100644 function_eba/threeDmulti_rot.m
create mode 100644 function_gfr/batch_processor_partitioned_pos.m
create mode 100644 function_gfr/batch_processor_partitioned_pos_rrho.m
create mode 100644 function_gfr/batch_processor_partitioned_rrho.m
create mode 100644 function_gfr/batch_processor_partitioned_rrho_tickonov.m
create mode 100644 function_gfr/hmat_pos.m
create mode 100644 function_gfr/save_vars_for_plotting.m
rename continue_gfr_par_from_iteration.m => release1/continue_gfr_par_from_iteration.m (92%)
rename gfr_parallel.m => release1/gfr_parallel.m (85%)
rename save_vars2continue_itr.m => release1/save_vars2continue_itr.m (70%)
create mode 100644 release2/gfr_eba.m
create mode 100644 release2/gfr_pos.m
create mode 100644 release2/gfr_pos_rrho.m
create mode 100644 release2/gfr_rrho.m
create mode 100644 release2/gfr_rrho_tickonov.m
create mode 100644 release2/results/.gitkeep
diff --git a/.gitignore b/.gitignore
index 050fd26..b5a48cf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -228,3 +228,5 @@ TSWLatexianTemp*
# End of https://www.gitignore.io/api/latex
variational.pdf
+input_data/
+results/
diff --git a/README.md b/README.md
index 927fcc1..63f9514 100644
--- a/README.md
+++ b/README.md
@@ -3,7 +3,8 @@
## Overview
GRACETOOLS is a collection of MATLAB scripts that can be used for gravity field
-recovery using GRACE type satellite observations.
+recovery using GRACE type satellite observations. There are two releases so
+far.
**Features**
@@ -19,10 +20,19 @@ recovery using GRACE type satellite observations.
<!-- vim-markdown-toc GitLab -->
- [Requirements](#requirements)
-- [Installation](#installation)
+- [Download](#download)
+ - [Get the Scripts](#get-the-scripts)
+ - [Get the Data](#get-the-data)
- [Usage](#usage)
-- [Performance](#performance)
+ - [Release 1](#release-1)
+ - [Performance](#performance)
+ - [Reference](#reference)
+ - [Release 2](#release-2)
+ - [Performance](#performance-1)
+ - [Reference](#reference-1)
- [Files](#files)
+ - [From `function_gfr`](#from-function_gfr)
+ - [From `function_eba`](#from-function_eba)
- [Contributors](#contributors)
<!-- vim-markdown-toc -->
@@ -38,7 +48,9 @@ Please reference this article if you use GRACETOOLS:
MATLAB and the [Parallel Computing
Toolbox](https://de.mathworks.com/help/distcomp/).
-## Installation
+## Download
+
+### Get the Scripts
Download GRACETOOLS as a zip file and extract it somewhere on your computer.
@@ -50,12 +62,7 @@ Or clone the repository with `git clone`:
$ git clone git@gitlab.aei.uni-hannover.de:geoq/gracetools.git
```
-## Usage
-
-The main m-files to run are:
-
-- `gfr_parallel.m`
-- `continue_gfr_par_from_iteration.m`
+### Get the Data
You need an example dataset to run these scripts. If you're on Linux or MacOS
use the `get_input_data.sh` script to get the dataset:
@@ -73,8 +80,8 @@ now unnecessary `input_data.zip` file. For completeness here are the terminal
commands:
```
-$ mv ~/Downloads/input_data.zip input_data/
-$ cd input_data/
+$ mv ~/Downloads/input_data.zip path/to/GRACETOOLS/input_data/
+$ cd path/to/GRACETOOLS/input_data/
$ unzip input_data.zip
$ rm -f input_data.zip
```
@@ -83,6 +90,8 @@ The directory structure of `input_data/` should look like this:
```
input_data
+├── EGM96.gfc
+├── GGM05S.gfc
├── MockData_do10_nod4
│ ├── 2005-05-01
│ │ ├── GNV1B_2005-05-01_A_02.asc
@@ -100,17 +109,110 @@ input_data
│ ├── GNV1B_2005-05-04_A_02.asc
│ ├── GNV1B_2005-05-04_B_02.asc
│ └── KBR1B_2005-05-04_X_02.asc
-├── EGM96.gfc
-└── GGM05S.gfc
+└── MockData_do20_nod12
+ ├── 2005-05-01
+ │ ├── GNV1B_2005-05-01_A_02.asc
+ │ ├── GNV1B_2005-05-01_B_02.asc
+ │ └── KBR1B_2005-05-01_X_02.asc
+ ├── 2005-05-02
+ │ ├── GNV1B_2005-05-02_A_02.asc
+ │ ├── GNV1B_2005-05-02_B_02.asc
+ │ └── KBR1B_2005-05-02_X_02.asc
+ ├── 2005-05-03
+ │ ├── GNV1B_2005-05-03_A_02.asc
+ │ ├── GNV1B_2005-05-03_B_02.asc
+ │ └── KBR1B_2005-05-03_X_02.asc
+ ├── 2005-05-04
+ │ ├── GNV1B_2005-05-04_A_02.asc
+ │ ├── GNV1B_2005-05-04_B_02.asc
+ │ └── KBR1B_2005-05-04_X_02.asc
+ ├── 2005-05-05
+ │ ├── GNV1B_2005-05-05_A_02.asc
+ │ ├── GNV1B_2005-05-05_B_02.asc
+ │ └── KBR1B_2005-05-05_X_02.asc
+ ├── 2005-05-06
+ │ ├── GNV1B_2005-05-06_A_02.asc
+ │ ├── GNV1B_2005-05-06_B_02.asc
+ │ └── KBR1B_2005-05-06_X_02.asc
+ ├── 2005-05-07
+ │ ├── GNV1B_2005-05-07_A_02.asc
+ │ ├── GNV1B_2005-05-07_B_02.asc
+ │ └── KBR1B_2005-05-07_X_02.asc
+ ├── 2005-05-08
+ │ ├── GNV1B_2005-05-08_A_02.asc
+ │ ├── GNV1B_2005-05-08_B_02.asc
+ │ └── KBR1B_2005-05-08_X_02.asc
+ ├── 2005-05-09
+ │ ├── GNV1B_2005-05-09_A_02.asc
+ │ ├── GNV1B_2005-05-09_B_02.asc
+ │ └── KBR1B_2005-05-09_X_02.asc
+ ├── 2005-05-10
+ │ ├── GNV1B_2005-05-10_A_02.asc
+ │ ├── GNV1B_2005-05-10_B_02.asc
+ │ └── KBR1B_2005-05-10_X_02.asc
+ ├── 2005-05-11
+ │ ├── GNV1B_2005-05-11_A_02.asc
+ │ ├── GNV1B_2005-05-11_B_02.asc
+ │ └── KBR1B_2005-05-11_X_02.asc
+ ├── 2005-05-12
+ │ ├── GNV1B_2005-05-12_A_02.asc
+ │ ├── GNV1B_2005-05-12_B_02.asc
+ │ └── KBR1B_2005-05-12_X_02.asc
+ └── states0_perfect.mat
+
```
-## Performance
+## Usage
+
+### Release 1
+
+The main m-files in the first release (`release1/`) are:
+
+- `gfr_parallel.m`
+- `save_vars2continue_itr.m`
+- `continue_gfr_par_from_iteration.m`
+
+First run `gfr_parallel.m`. Use `save_vars2continue_itr.m` to save all necessary
+variables in the `results/` subdirectory. You can continue the estimation from
+the last iteration with `continue_gfr_par_from_iteration.m`.
+
+#### Performance
With an older Intel i5 (4 cores) a test case with degree and order 10 with 4
-days of observation, each iteration takes about 10min.
+days of observation, each iteration takes about 10 min.
+
+#### Reference
+
+- Darbeheshti N., Wöske F., Weigelt M., McCullough C., Wu H. (2018) GRACETOOLS - GRACE gravity field recovery tools, Geosciences 2018, 8(9), 350; https://doi.org/10.3390/geosciences8090350.
+
+### Release 2
+
+The main m-files in the second release (`release2/`) are:
+
+- `gfr_eba.m` - GFR with energy balance approach
+- `gfr_pos.m` - GFR with variational equations using positions
+- `gfr_pos_rrho.m` - GFR with variational equations using positions and range rates
+- `gfr_rrho.m` - GFR with variational equations using range rates
+- `gfr_rrho_tickonov.m`- GFR with variational equations using range rates and Tickonov regularization
+
+Navigate to the `release2/` subdirectory and run the scripts. Results will be
+stored in the `results` subdirectory.
+
+#### Performance
+
+With an older Intel i5 (4 cores) a test case with degree and order 20 with 12
+days of observation, each iteration takes about 20 min.
+
+#### Reference
+
+- Darbeheshti N., Wöske F., Weigelt M., Wu H., Mccullough C. (2020) Comparison of Spacewise and Timewise Methods for GRACE Gravity Field Recovery. In: Montillet JP., Bos M. (eds) Geodetic Time Series Analysis in Earth Sciences. Springer Geophysics. Springer, Cham, https://doi.org/10.1007/978-3-030-21718-1_10
+
+- [IUGG Talk Montreal 2019](https://seafile.projekt.uni-hannover.de/f/2668157cf326417485c7/)
## Files
+### From `function_gfr`
+
- `batch_processor_partitioned.m` - batch processing algorithm for GRACE
range-rate observations
- `cs2sc.m`- converts the square containing spherical harmonics coefficients
@@ -162,6 +264,13 @@ days of observation, each iteration takes about 10min.
- `vec2cs.m`- rearranges a vector shaped set of spherical harmonic coefficients
into cs-format
+### From `function_eba`
+
+- `calculatePnm.m` - Calculation of normalized Legendre Functions using stable recursion formulas
+- `designmatrixn.m` - Calculation of normalized Legendre Functions using stable recursion formulas
+- `potential_EBAs.m` - Calculates the potential along the GRACE simulated orbit
+- `sortCoefficients.m` - Sorts spherical harmonic coefficients
+
## Contributors
- [Neda Darbeheshti](mailto:neda.darbeheshti@aei.mpg.de)
diff --git a/function_eba/calculatePnm.m b/function_eba/calculatePnm.m
new file mode 100644
index 0000000..8c85758
--- /dev/null
+++ b/function_eba/calculatePnm.m
@@ -0,0 +1,41 @@
+function Pnm = calculatePnm(theta, n_max)
+%
+% Calculation of normalized Legendre Functions using stable recursion
+% formulas.
+%
+% --Input--
+% theta Co-Latitude (1x1)
+% n_max Maximum degree (1x1)
+%
+% --Output--
+% Pnm Legendre functions (n_max x n_max)
+%
+% Author: Anne Springer
+% Date: 2015-08-24
+%
+
+t = cos(theta);
+
+Pnm = zeros(n_max+1, n_max+1);
+Pnm(1,1) = 1;
+an = sqrt(3);
+
+for n = 1:n_max
+ if n>1
+ an = sqrt((2*n+1)/(2*n));
+ end
+
+ % P_{n,n}
+ Pnm(n+1, n+1) = an * sqrt(1-t^2)*Pnm(n,n);
+
+ % P_{n,n-1}
+ bn = sqrt((2*n + 1)*(2*n - 1)/((n+(n-1))*(n-(n-1))));
+ Pnm(n+1, n) = bn*Pnm(n,n)*t;
+
+ % P_{n,m}
+ for m = 0:n-2
+ bnm = sqrt((2*n + 1)*(2*n - 1)/((n+m)*(n-m)));
+ cnm = sqrt((2*n+1)*(n-m-1)*(n+m-1)/((2*n-3)*(n+m)*(n-m))) ;
+ Pnm(n+1, m+1) = bnm*t*Pnm(n, m+1) - cnm*Pnm(n-1, m+1);
+ end
+end
\ No newline at end of file
diff --git a/function_eba/designmatrixn.m b/function_eba/designmatrixn.m
new file mode 100644
index 0000000..b669957
--- /dev/null
+++ b/function_eba/designmatrixn.m
@@ -0,0 +1,41 @@
+function A = designmatrixn(longitude,latitude,height,n_max,GM,R)
+%
+% Calculation of normalized Legendre Functions using stable recursion
+% formulas.
+%
+% --Input--
+% longitude Longitudes (n x 1) in Radian
+% latitudes Latitudes (n x 1) in Radian
+% height Heights (n x 1)
+% n_max Maximum degree (1 x 1)
+% GM Gravitational constant times mass of the Earth
+% R Radius of the Earth
+%
+% --Output--
+% A Design matrix (n x n_max)
+%
+% Author: Anne Springer
+% Date: 2015-08-24
+%
+
+% Degree and order
+idx = tril(ones(n_max+1))~=0;
+degree = repmat((0:n_max)',1,n_max+1);
+degree = degree(idx);
+order = repmat(0:n_max,n_max+1,1);
+order = order(idx);
+
+% Design matrix
+A = zeros(length(longitude),length(degree)*2-(n_max+1));
+for i = 1:length(A)
+ Pnm = calculatePnm((pi/2-latitude(i)), n_max);
+ Pnm = Pnm(idx);
+ factor = GM/R*(R/height(i)).^(degree+1);
+ cos_ = cos(order*longitude(i));
+ sin_ = sin(order*longitude(i));
+ A(i,:) = [factor.*cos_.*Pnm; factor(n_max+2:end).*sin_(n_max+2:end).*Pnm(n_max+2:end)]';
+end
+
+
+
+
diff --git a/function_eba/potential_EBAs.m b/function_eba/potential_EBAs.m
new file mode 100644
index 0000000..65ae2ff
--- /dev/null
+++ b/function_eba/potential_EBAs.m
@@ -0,0 +1,41 @@
+function v12sst = potential_EBAs(rho_dot,s_ea,s_eb,we)
+% POTENTIAL_EBAS.m calculates the potential along the GRACE simulated orbit.
+
+% rho_dot is the range rate
+% s_ea and s_eb are the state vectors, position and velocity in itrf
+% we is the mean Earth rotation rate
+%
+% Written by Neda Darbeheshti, AEI Hannover, 2018-12.
+
+global e12 en er
+
+xdot12=s_eb(:,4:6)-s_ea(:,4:6);
+sumxdot=s_eb(:,4:6)+s_ea(:,4:6);
+
+%unit vectors
+
+e12=s_eb(:,1:3)-s_ea(:,1:3);
+e12 = bsxfun(@rdivide,e12,sqrt(sum(e12.^2,2)));
+
+en = [s_ea(:,2).*s_eb(:,3)-s_ea(:,3).*s_eb(:,2) s_ea(:,3).*s_eb(:,1)-s_ea(:,1).*s_eb(:,3) s_ea(:,1).*s_eb(:,2)-s_ea(:,2).*s_eb(:,1)];
+en = bsxfun(@rdivide,en,sqrt(sum(en.^2,2)));
+
+er = [e12(:,2).*en(:,3)-e12(:,3).*en(:,2) e12(:,3).*en(:,1)-e12(:,1).*en(:,3) e12(:,1).*en(:,2)-e12(:,2).*en(:,1)];
+er = bsxfun(@rdivide,er,sqrt(sum(er.^2,2)));
+
+e3=[0 0 1];
+
+% calculate potential along each vector
+for i=1:length(rho_dot)
+
+v12sst_rho(i)=.5*rho_dot(i)*e12(i,:)*sumxdot(i,:)';
+v12sst_n(i)=.5*en(i,:)*xdot12(i,:)'*en(i,:)*sumxdot(i,:)';
+v12sst_r(i)=.5*er(i,:)*xdot12(i,:)'*er(i,:)*sumxdot(i,:)';
+
+v12sst_rota(i)=-we^2/2*(norm(cross(e3,s_eb(i,1:3)))^2-norm(cross(e3,s_ea(i,1:3)))^2);
+
+end
+
+% sum all potentials
+v12sst=v12sst_n+v12sst_r+v12sst_rho+v12sst_rota;
+end
diff --git a/function_eba/sortCoefficients.m b/function_eba/sortCoefficients.m
new file mode 100644
index 0000000..ffa9714
--- /dev/null
+++ b/function_eba/sortCoefficients.m
@@ -0,0 +1,38 @@
+function [cnm, snm] = sortCoefficients(x, n_max)
+
+%
+% SORTCOEFFICIENTS.m sorts spherical harmonic coefficients
+%
+% --Input--
+% x Parameter vector containing the estimated spherical
+% harmonic coefficients (n x 1)
+% n_max Maximum degree (1x1)
+%
+% --Output--
+% cnm, snm Potential coefficients divided into sine component and
+% cosine component
+%
+% Author: Anne Springer
+% Date: 2015-08-24
+%
+
+% Degree and order
+idx = tril(ones(n_max+1))~=0;
+degree = repmat((0:n_max)',1,n_max+1);
+degree = degree(idx);
+order = repmat(0:n_max,n_max+1,1);
+order = order(idx);
+
+% Cosine coefficients
+cnm = zeros(n_max+1);
+for i = 1:length(degree)
+ cnm(degree(i)+1, order(i)+1) = x(i);
+end
+
+% Sine coefficients
+snm = zeros(n_max+1);
+for j = 1:length(degree(n_max+2:end))
+ snm(degree(j+n_max+1)+1, order(j+n_max+1)+1) = x(i+j);
+end
+
+end
diff --git a/function_eba/threeDmulti_rot.m b/function_eba/threeDmulti_rot.m
new file mode 100644
index 0000000..b6d4ee3
--- /dev/null
+++ b/function_eba/threeDmulti_rot.m
@@ -0,0 +1,45 @@
+function C = threeDmulti_rot(A,B,type)
+
+% This function multiplies 3 dimensional rotation matrices in form A(3*3*n)
+% with position vectors in form B(n*3). This is mainly used in coordinate
+% transformation between intertial and earth-fixed frames.
+%
+% Input: R = rotation matrices [3x3xN]
+% B = vectors to be rotated [Nx3]
+% type = 'forward' or 'backward' ('backward' for inverse rotation);
+% default is 'forward'
+%
+% Output: C = rotated vectors [Nx3]
+%
+% Example: R = RotmatICRF2SF(SCA1B_A(:,2:5));
+% v = acc_ICRF;
+% --> acc_SF = threeDmulti_rot(R,v);
+%
+% Written by Majid Naeimi, 7 April 2015, IfE, Hanover, Germany
+% Last updated by Henry Wegener, AEI Hannover, 2016-07-08.
+
+
+if ~( size(A,1) == 3 && size(A,2) == 3 && size(B,2) == 3 )
+ error('Dimensions of input arguments are incorrect. Rotation matrices (first argument) must be of size 3x3xN \n Input vectors (second argument) must be of size Nx3 !')
+elseif ( size(A,3) > 1 && size(B,1) == 1 )
+ B = repmat(B,size(A,3),1);
+elseif ~( size(A,3) == size(B,1) )
+ error('Length of vector (second argument) must be the same as for the matrices OR it must be 1 !');
+end
+
+if ( nargin == 2 )
+ type = 'forward';
+end
+
+Ap = reshape(A,9,[]);
+switch type
+ case 'forward'
+ C = [Ap(1,:)'.*B(:,1)+Ap(4,:)'.*B(:,2)+Ap(7,:)'.*B(:,3) Ap(2,:)'.*B(:,1)+Ap(5,:)'.*B(:,2)+Ap(8,:)'.*B(:,3) Ap(3,:)'.*B(:,1)+Ap(6,:)'.*B(:,2)+Ap(9,:)'.*B(:,3)];
+ case 'backward'
+ C = [Ap(1,:)'.*B(:,1)+Ap(2,:)'.*B(:,2)+Ap(3,:)'.*B(:,3) Ap(4,:)'.*B(:,1)+Ap(5,:)'.*B(:,2)+Ap(6,:)'.*B(:,3) Ap(7,:)'.*B(:,1)+Ap(8,:)'.*B(:,2)+Ap(9,:)'.*B(:,3)];
+ otherwise
+ error(['Unknown type: ',type,'. type must be either "forward" (for forward rotation) \n or "backward" (for inverse rotation)!'])
+end
+
+
+end
\ No newline at end of file
diff --git a/function_gfr/batch_processor_partitioned_pos.m b/function_gfr/batch_processor_partitioned_pos.m
new file mode 100644
index 0000000..7bfeb0e
--- /dev/null
+++ b/function_gfr/batch_processor_partitioned_pos.m
@@ -0,0 +1,184 @@
+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
+
+
diff --git a/function_gfr/batch_processor_partitioned_pos_rrho.m b/function_gfr/batch_processor_partitioned_pos_rrho.m
new file mode 100644
index 0000000..28a9941
--- /dev/null
+++ b/function_gfr/batch_processor_partitioned_pos_rrho.m
@@ -0,0 +1,201 @@
+function batch_processor_partitioned_pos_rrho(Mday,nods,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state,timeKBR,obs_gnv,obs_kbr,WMAT)
+% Batch processing algorithm for GRACE using positions and range-rate 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_rrho.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
+% numbers of observations for each time step
+no = size(obs_gnv,2); %(3 or 6 positions, positios + velocities)
+
+% 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;
+ Mxx=zeros(12);
+ 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'ygnv_save','ykbr_save','Hx_save','Hc_save' just for last day
+if Mday==nods
+ ykbr_save = zeros(mKBR,1);
+ ygnv_save = zeros(mKBR,no);
+ Hx_save = zeros(mKBR,no+1,12);
+ Hc_save = zeros(mKBR,no+1,noc);
+ %Run through each observation:
+for ii = 1:mKBR
+ % Read observation:
+ Ykbr = obs_kbr(ii,:)';
+ Ygnv = obs_gnv(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_kbr,Y_star_kbr]= hmat(X,noc);
+
+ % Get the Htilde Matrix:
+ [Htilde_gnv,Y_star_gnv]= hmat_pos(X,noc,no);
+
+ %builld Htilde
+ Htilde=[Htilde_kbr; Htilde_gnv];
+
+ % Time update:
+ H = Htilde*phibuff;
+ Hx = H(:,1:12);
+ Hc = H(:,13:end);
+
+ % Calculate the observation deviation:
+ ykbr = Ykbr - Y_star_kbr(2,:);
+ ygnv = Ygnv - Y_star_gnv;
+
+ %builld y
+ y=[ykbr;ygnv];
+
+ % accumulate Lambda:
+ %L = L + H'*WMAT*H;
+ Mxx=Mxx + Hx'*WMAT*Hx;
+ Mxc=Mxc + Hx'*WMAT*Hc;
+ Mcc=Mcc + Hc'*WMAT*Hc;
+ % accumulate N:
+ %N = N + H'*WMAT*y;
+ Nx = Nx + Hx'*WMAT*y;
+ Nc = Nc + Hc'*WMAT*y;
+ % Store data:
+ ykbr_save(ii) = ykbr;
+ ygnv_save(ii,:) = ygnv;
+ 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','ygnv_save','Hx_save','Hc_save','-v7.3');
+else
+ %Run through each observation:
+for ii = 1:mKBR
+ % Read observation:
+ Ykbr = obs_kbr(ii,:)';
+ Ygnv = obs_gnv(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_kbr,Y_star_kbr]= hmat(X,noc);
+
+ % Get the Htilde Matrix:
+ [Htilde_gnv,Y_star_gnv]= hmat_pos(X,noc,no);
+
+ %builld Htilde
+ Htilde=[Htilde_kbr; Htilde_gnv];
+
+ % Time update:
+ H = Htilde*phibuff;
+ Hx = H(:,1:12);
+ Hc = H(:,13:end);
+
+ % Calculate the observation deviation:
+ ykbr = Ykbr - Y_star_kbr(2,:);
+ ygnv = Ygnv - Y_star_gnv;
+
+ %builld y
+ y=[ykbr;ygnv];
+
+ % accumulate Lambda:
+ %L = L + H'*WMAT*H;
+ Mxx=Mxx + Hx'*WMAT*Hx;
+ Mxc=Mxc + Hx'*WMAT*Hc;
+ Mcc=Mcc + Hc'*WMAT*Hc;
+ % accumulate N:
+ %N = N + H'*WMAT*y;
+ Nx = Nx + Hx'*WMAT*y;
+ Nc = Nc + Hc'*WMAT*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
+
+
diff --git a/function_gfr/batch_processor_partitioned_rrho.m b/function_gfr/batch_processor_partitioned_rrho.m
new file mode 100644
index 0000000..c464d0f
--- /dev/null
+++ b/function_gfr/batch_processor_partitioned_rrho.m
@@ -0,0 +1,179 @@
+function batch_processor_partitioned_rrho(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.
+% 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.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;
+ %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'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
\ No newline at end of file
diff --git a/function_gfr/batch_processor_partitioned_rrho_tickonov.m b/function_gfr/batch_processor_partitioned_rrho_tickonov.m
new file mode 100644
index 0000000..c7127da
--- /dev/null
+++ b/function_gfr/batch_processor_partitioned_rrho_tickonov.m
@@ -0,0 +1,182 @@
+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
\ No newline at end of file
diff --git a/function_gfr/hmat_pos.m b/function_gfr/hmat_pos.m
new file mode 100644
index 0000000..021ca35
--- /dev/null
+++ b/function_gfr/hmat_pos.m
@@ -0,0 +1,42 @@
+function [Hi_tilda, Y_star] = hmat_pos(X0,NOCs,no)
+
+% This function makes the Htilde matrix.
+% References:-McCullough CM (2017) Gravity Field Estimation for Next Generation Satellite
+% Missions, page 97, equation (4.4)
+% Liu(2008), page 28, equation (2.51)
+% -Statistical Orbit Determination (Tapley et al., 2004)
+%
+% Input: X0 = state vector including [RA; RdotA; PHI]
+% NOCs = total number of coefficients to be estimated
+% no = number of observations:3 (for positions) or 6 for (positios + velocities)
+%
+% Outputs: Hi_tilda = [dRA/dRA; dRA/dRdotA; dRA/dRB; dRA/dRdotB; dRA/dKlm]
+%
+%
+% Written by Florian Wöske, ZARM Uni Bremen,2018-11
+
+%% calculate gnv partials w.r.t SHC
+ % Extract the current state transition matrix:
+ phi = reshape(X0(12+1:end), 12, 12+NOCs);
+ % Extract the partials from the current state transition matrix:
+ % if no == 3 always positions
+ dr1_dklm=phi(1:3,13:end);
+ if no == 6
+ dr2_dklm=phi(7:9,13:end);
+ end
+
+ I = eye(3,3);
+ Z = zeros(3,3);
+
+ % fill Hi_tilda matrix
+ if no == 3
+ Hi_tilda=[I, Z, Z, Z, dr1_dklm];
+ % observation is directly x,y,z
+ Y_star = X0(1:3);
+ elseif no == 6
+ Hi_tilda=[I, Z, Z, Z, dr1_dklm; Z, Z, I, Z, dr2_dklm];
+ % observation is directly x,y,z
+ Y_star =[X0(1:3); X0(7:9)];
+ end
+
+end
diff --git a/function_gfr/save_vars_for_plotting.m b/function_gfr/save_vars_for_plotting.m
new file mode 100644
index 0000000..0c367c9
--- /dev/null
+++ b/function_gfr/save_vars_for_plotting.m
@@ -0,0 +1,17 @@
+%% Save all important variables of a gravity estimation
+% this can be used to save all important results, as well as to continue
+% the estimation from the last iteration by eg. using
+% continue_gfr_par_from_iteration.m .
+% All variable are saved in one .mat file, and can be directly loaded in
+% continue_gfr_par_from_iteration.m.
+%%
+% Filename and path of file to save
+FileName= ['results/','GFR_itr_',num2str(ItrNo_max), '_do_', num2str(lmaxcs), '_nods_', num2str(nods), '_','tno_', num2str(tno),'_', datestr(now, 'dd-mmm-yyyy'),'.mat'];
+
+% save all necessary variables in one file
+save(FileName,'xhat_save','chat_save','eps_save','chat_dv_save','dv_save','state0','x0x','field', 'field_cs', 'field_vec0', 'mKBR', 'GM', 'ae', 't0', 'nods', 'lmaxcs', 'lmaxf', 'data_folder', 'state00', 't_itr_nod','tno');
+
+% actual time
+c = clock;
+str_display_out = ['results saved as: ', FileName,', Time: ',num2str(c(1)),'-',num2str(c(2)),'-',num2str(c(3)),', ',num2str(c(4)),':',num2str(c(5)),':',num2str(c(6)) ];
+disp(str_display_out);
diff --git a/get_input_data.sh b/get_input_data.sh
index 1d637cc..edbc057 100755
--- a/get_input_data.sh
+++ b/get_input_data.sh
@@ -1,16 +1,20 @@
#!/bin/bash
+#
+# This script downloads the input data. You need wget und unzip to run it.
+#
+# ./get_input_data.sh
+#
URL="https://seafile.projekt.uni-hannover.de/f/3ee8fef6c6cb485489b4/?dl=1"
INPUT_DATA=input_data.zip
-downloadData(){
+download_data(){
echo "Downloading data from $URL..."
wget --content-disposition "url?raw=1" -U mozilla "$URL" -O $INPUT_DATA
echo "Done."
}
-
-setUp(){
+set_up(){
echo "Start extracting data file..."
unzip $INPUT_DATA -d input_data/
echo "Done."
@@ -19,7 +23,7 @@ setUp(){
echo "Done."
}
-cleanUp(){
+clean_up(){
echo "Removing subdirectories from input directory."
rm -rf input_data/*
echo "Done."
@@ -28,7 +32,7 @@ cleanUp(){
echo "Done."
}
-getData(){
+get_data(){
if [[ -d input_data/ ]]; then
echo "Directory 'input_data/' already exists."
if [[ -f input_data/EGM96.gfc ]]; then
@@ -36,9 +40,9 @@ getData(){
read -r -p "Cleanup and download again? [y/N]" response
case $response in
[yY][eE][sS]|[yY])
- cleanUp
- downloadData
- setUp
+ clean_up
+ download_data
+ set_up
;;
*)
echo "Aborting."
@@ -46,18 +50,17 @@ getData(){
;;
esac
else
- downloadData
- setUp
+ download_data
+ set_up
fi
else
- downloadData
- setUp
+ download_data
+ set_up
fi
}
-
main(){
- getData
+ get_data
}
main
diff --git a/continue_gfr_par_from_iteration.m b/release1/continue_gfr_par_from_iteration.m
similarity index 92%
rename from continue_gfr_par_from_iteration.m
rename to release1/continue_gfr_par_from_iteration.m
index 0d559f5..7dc1e3d 100644
--- a/continue_gfr_par_from_iteration.m
+++ b/release1/continue_gfr_par_from_iteration.m
@@ -6,7 +6,6 @@
%
% Written by Florian Wöske, ZARM Uni Bremen, 2018-07
-
%%
% variables that are continued, and thus need to be known:
% xhat_save (12, nods, ItrNo_max+1)
@@ -27,31 +26,32 @@
% additional needed variables:
% mKBR, GM, ae, t0, nods, lmaxcs, lmaxf
-% filename of file with saved variables:
-FileName = 'GFR_Iteration_12_do_10_nods_4_2005-05-04_011.mat';
+% Result file from gfr_parallel.m:
+GFR_Iteration_Result = 'results/GFR_Iteration_12_do_10_nods_4_2005-05-04_011.mat';
%% variable inputs
+
% add folder with all functions
-FolderName = pwd;
-addpath(genpath(FolderName));
+Release1_Path = pwd;
+addpath(genpath(Release1_Path));
% folder to save intermediate results
mkdir('results/OutputArcParall');
mkdir('results/OutputNormalEqu');
+% add folder name containing functions and input data (observations, ref.
+% field, initial values...)
+Input_Data_Path = fullfile(Release1_Path, '/../input_data');
+addpath(genpath(Input_Data_Path));
+
% Define the whole number of iterations for batch processor (previous + new
% iterations)
ItrNo_max = 12;
-% add folder name containing functions and input data (observations, ref.
-% field, initial values...)
-DataFolder='input_data';
-addpath(genpath(DataFolder));
-
%% load data from previous iteration
% load data
-load(FileName);
+load(GFR_Iteration_Result);
% last iteration:
itr_prev = size(xhat_save,3)-1;
@@ -68,7 +68,6 @@ n_back= ((lmaxf+1)^2 + lmaxf+1)/2;
ggm05s = importGravityField('GGM05S.gfc');
ggm05s = ggm05s(1:n_back,:); % cut to desired length
-
% extend variables for additional iterations:
xhat_save_prev = xhat_save;
chat_save_prev = chat_save;
@@ -111,7 +110,7 @@ for Mday = 1:nods
% load KBR data
data_KBR = ['KBR1B_', date, '_X_02.asc'];
- data_file = fullfile(data_folder, date, data_KBR);
+ data_file = fullfile(Input_Observation_Path, date, data_KBR);
KBRL1Bi = readKBR(data_file);
KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
@@ -119,13 +118,12 @@ end
%% iteration loop
-
for ItrNo = (itr_prev+1):ItrNo_max
% batch over all days
parfor Mday = 1:nods
tic
- batch_processor_partitioned(Mday,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
+ batch_processor_partitioned(Mday,Release1_Path,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
toc
end
@@ -144,7 +142,7 @@ for ItrNo = (itr_prev+1):ItrNo_max
irhodot_deviation = zeros(mKBR,nods);
for Mday = 1:nods
FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
- File = fullfile(FolderName,'results/OutputArcParall', FileName);
+ File = fullfile(Release1_Path,'results/OutputArcParall', FileName);
arcstru=load(File);
invMxx(:,:,Mday) = arcstru.invMxx;
iMxc(:,:,Mday) = arcstru.Mxc;
@@ -215,7 +213,7 @@ for ItrNo = (itr_prev+1):ItrNo_max
% % save final monthly gravity solution in a file
% FileName= ['gfr_NODs',num2str(nods,'%.2d'),'_DO',num2str(lmaxcs,'%.2d'),'_ItrNo',num2str(ItrNo,'%.2d'),'.gfc'];
-% File = fullfile(FolderName,'results/OutputNormalEqu', FileName);
+% File = fullfile(Release1_Path,'results/OutputNormalEqu', FileName);
% fileID = fopen(File,'w');
% for i=1:length(ordering2)
% fprintf(fileID,'%4i %4i %17.16e %17.16e\n',ordering2(i,:));
diff --git a/gfr_parallel.m b/release1/gfr_parallel.m
similarity index 85%
rename from gfr_parallel.m
rename to release1/gfr_parallel.m
index fe055ce..b270dbd 100644
--- a/gfr_parallel.m
+++ b/release1/gfr_parallel.m
@@ -18,13 +18,22 @@ format longg;
%% variable inputs
% add folder with all functions
-FolderName = pwd;
-addpath(genpath(FolderName));
+Release1_Path = pwd;
+addpath(genpath(Release1_Path));
+
+% EBA_Path = fullfile(Release1_Path,'/../function_eba/');
+% addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release1_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
% folder to save intermediate results
mkdir('results/OutputArcParall');
mkdir('results/OutputNormalEqu');
+% folder where observation data are stored
+Input_Observation_Path = fullfile(Release1_Path,'/../input_data','MockData_do10_nod4');
+addpath(genpath(Input_Observation_Path));
+
% Define the number of iterations for batch processor,
ItrNo_max = 5;
@@ -37,10 +46,6 @@ lmaxf = 10;
% number of days
nods = 4;
-% folder where observation data are stored
-data_folder = fullfile(FolderName,'input_data','MockData_do10_nod4');
-addpath(genpath(data_folder));
-
% number of observations in each arc (one day, 17280 is the standard)
mKBR = 17280;
@@ -55,14 +60,14 @@ GM = 0.3986004415E+15;
ae = 0.6378136300E+07;
% reference for comparison
-ggm05s = readGFC('GGM05S.gfc');
+ggm05s = readGFC('../input_data/GGM05S.gfc');
ggm05s = ggm05s(1:n_back,:); % cut to desired length
ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
% d_vec_i = ggm05s(:,3:4)-gravityf(:,3:4);
% d_vec_i = [ggm05s(:,1:2) d_vec_i];
-%Load a-priori/ reference gravity field
-egm96 = readGFC('EGM96.gfc');
+% Load a-priori/ reference gravity field
+egm96 = readGFC('../input_data/EGM96.gfc');
egm96 = egm96(1:n_back,:); % cut to desired length
egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
@@ -101,26 +106,25 @@ for Mday = 1:nods
% load KBR data
data_KBR = ['KBR1B_', date, '_X_02.asc'];
- data_file = fullfile(data_folder, date, data_KBR);
+ data_file = fullfile(Input_Observation_Path, date, data_KBR);
KBRL1Bi = readKBR(data_file);
KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
% load GNV data GRACE A
data_GNV = ['GNV1B_', date, '_A_02.asc'];
- data_file = fullfile(data_folder, date, data_GNV);
+ data_file = fullfile(Input_Observation_Path, date, data_GNV);
GNVi = readGNV(data_file);
state0(1:6,Mday) = GNVi(1,2:7);
% load GNV data GRACE B
data_GNV = ['GNV1B_', date, '_B_02.asc'];
- data_file = fullfile(data_folder, date, data_GNV);
+ data_file = fullfile(Input_Observation_Path, date, data_GNV);
GNVi = readGNV(data_file);
state0(7:12,Mday) = GNVi(1,2:7);
-
end
- % save initial state vector for comparisons
- state00 = state0;
+% save initial state vector for comparisons
+state00 = state0;
%% iteration
@@ -140,7 +144,7 @@ for ItrNo = 0:ItrNo_max
% batch over all days
parfor Mday = 1:nods
tic
- batch_processor_partitioned(Mday,FolderName,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
+ batch_processor_partitioned(Mday,Release1_Path,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
toc
end
@@ -159,7 +163,7 @@ for ItrNo = 0:ItrNo_max
irhodot_deviation = zeros(mKBR,nods);
for Mday = 1:nods
FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
- File = fullfile(FolderName,'results/OutputArcParall', FileName);
+ File = fullfile(Release1_Path,'results/OutputArcParall', FileName);
arcstru=load(File);
invMxx(:,:,Mday) = arcstru.invMxx;
iMxc(:,:,Mday) = arcstru.Mxc;
@@ -200,8 +204,6 @@ for ItrNo = 0:ItrNo_max
% title(['Range-rate residuals, Itr',num2str(ItrNo,'%.2d'),',Day',num2str(Mday,'%.2d')])
end
-
-
% save xhat and chat of iterations
xhat_save(:,:,ItrNo+1) = xhat;
chat_save(:,ItrNo+1) = chat;
@@ -228,16 +230,16 @@ for ItrNo = 0:ItrNo_max
chat_vec = [nm1',outC1',outS1'];
chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
-% dv_geoidn(chat_vec,lmaxcs);
+ % dv_geoidn(chat_vec,lmaxcs);
-% % save final monthly gravity solution in a file
-% FileName= ['gfr_NODs',num2str(nods,'%.2d'),'_DO',num2str(lmaxcs,'%.2d'),'_ItrNo',num2str(ItrNo,'%.2d'),'.gfc'];
-% File = fullfile(FolderName,'results/OutputNormalEqu', FileName);
-% fileID = fopen(File,'w');
-% for i=1:length(ordering2)
-% fprintf(fileID,'%4i %4i %17.16e %17.16e\n',ordering2(i,:));
-% end
-% fclose(fileID);
+ % % save final monthly gravity solution in a file
+ % FileName= ['gfr_NODs',num2str(nods,'%.2d'),'_DO',num2str(lmaxcs,'%.2d'),'_ItrNo',num2str(ItrNo,'%.2d'),'.gfc'];
+ % File = fullfile(Release1_Path,'results/OutputNormalEqu', FileName);
+ % fileID = fopen(File,'w');
+ % for i=1:length(ordering2)
+ % fprintf(fileID,'%4i %4i %17.16e %17.16e\n',ordering2(i,:));
+ % end
+ % fclose(fileID);
% set up values for next iteration
state0 = state0 + xhat;
@@ -258,7 +260,6 @@ for ItrNo = 0:ItrNo_max
dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
fprintf('done with iteration %d \n', ItrNo)
-
end
% plot chat from all iterations:
@@ -267,6 +268,7 @@ for i = 1:ItrNo_max+1
semilogy(2:lmaxcs,chat_dv_save(i,:),'.-','LineWidth',2,'MarkerSize',20);
hold on
end
+
fs = 12;
set(gcf,'PaperPositionMode','auto')
set(gca,'FontSize',fs);
@@ -291,6 +293,7 @@ legend('ggm05s', 'estimated field', 'ggm05 - estimated field')
% plot gravity field error from each iteration additional to initial fields and
% initial error
+
figure;
dv_geoidn(ggm05s,lmaxcs);
hold on
diff --git a/save_vars2continue_itr.m b/release1/save_vars2continue_itr.m
similarity index 70%
rename from save_vars2continue_itr.m
rename to release1/save_vars2continue_itr.m
index 37aee8b..4899ba8 100644
--- a/save_vars2continue_itr.m
+++ b/release1/save_vars2continue_itr.m
@@ -6,7 +6,7 @@
%%
% Filename and path of file to save
-FileName= ['GFR_Iteration_',num2str(ItrNo_max), '_do_', num2str(lmaxcs), '_nods_', num2str(nods), '_', date, '_011.mat'];
+FileName= ['results/','GFR_Iteration_',num2str(ItrNo_max), '_do_', num2str(lmaxcs), '_nods_', num2str(nods), '_', date, '_011.mat'];
% save all necessary variables in one file
-save(FileName,'xhat_save','chat_save','eps_save','chat_dv_save','dv_save','state0','x0x','field', 'field_cs', 'err_vec', 'field_vec0', 'mKBR', 'GM', 'ae', 't0', 'nods', 'lmaxcs', 'lmaxf', 'data_folder','state00');
+save(FileName,'xhat_save','chat_save','eps_save','chat_dv_save','dv_save','state0','x0x','field', 'field_cs', 'err_vec', 'field_vec0', 'mKBR', 'GM', 'ae', 't0', 'nods', 'lmaxcs', 'lmaxf', 'Input_Observation_Path','state00');
diff --git a/release2/gfr_eba.m b/release2/gfr_eba.m
new file mode 100644
index 0000000..184f4f3
--- /dev/null
+++ b/release2/gfr_eba.m
@@ -0,0 +1,151 @@
+%% GRACE gravity field recovery (GFR) using Energy Balance Approach (EBA)
+% This is the main m-file preforming the gravity field estimation.
+% All necessary parameters and data are specified in the following.
+%
+% Written by Neda Darbeheshti, AEI, 2019-06.
+
+clear all
+format longg;
+
+%% variable inputs
+global e12 en er
+
+% add folder with all functions
+Release2_Path = pwd;
+addpath(genpath(Release2_Path));
+EBA_Path = fullfile(Release2_Path,'/../function_eba/');
+addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release2_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
+
+% folder where observation data are stored
+% Use twelve days of data
+Input_Observation_Path = fullfile(Release2_Path,'/../input_data','MockData_do20_nod12');
+% Use four days of data (optional)
+% Input_Observation_Path = fullfile(Release2_Path,'/../input_data','MockData_do10_nod4');
+addpath(genpath(Input_Observation_Path));
+
+% order of spherical harmonic coefficients to be estimated
+n_max = 20;
+
+% number of days
+nod = 12;
+
+% number of observations in each arc (one day, 17280 is the standard)
+mKBR = 17280;
+
+%% fixed inputs
+
+% GM and Earth radius of gravity field
+a = 6378136.3; %m
+GM = 3.9860044150e+14;
+we0=7.29211514670698*10^(-5); % rad/s
+
+% size of gravity fields in vec format
+length_N=(n_max/2+1)*(n_max+1)*2-(n_max+1);
+
+% initialize for normal equations
+N = zeros(length_N+1,length_N+1);
+n0 = zeros(length_N+1,1);
+
+% initialze one column for energy constant
+A0 = ones(3600*24/5,1);
+
+%% load observations data for all days from files
+
+t0 = [2005 05 01 0 0 2];
+t0i = t0;
+tic
+figure
+
+for day=1:nod
+ if day >1
+ t0i(3)=t0i(3)+1;
+ end
+ date = time2str(t0i);
+
+ % load GNV data GRACE A
+ data_GNV = ['GNV1B_', date, '_A_02.asc'];
+ data_file = fullfile(Input_Observation_Path, date, data_GNV);
+ GNVA = readGNV(data_file);
+
+ % load GNV data GRACE B
+ data_GNV = ['GNV1B_', date, '_B_02.asc'];
+ data_file = fullfile(Input_Observation_Path, date, data_GNV);
+ GNVB = readGNV(data_file);
+
+ % load KBR data
+ data_KBR = ['KBR1B_', date, '_X_02.asc'];
+ data_file = fullfile(Input_Observation_Path, date, data_KBR);
+ KBRB = readKBR(data_file);
+
+ % transformation to the Earth_fixed frame
+ RotMat= Ri2e(GNVA(:,1));
+ RotMat_dot= Ri2e_dot(GNVA(:,1));
+ rA_e = threeDmulti_rot(RotMat,GNVA(:,2:4),'forward');
+ vA_e = threeDmulti_rot(RotMat,GNVA(:,5:7),'forward')+threeDmulti_rot(RotMat_dot,GNVA(:,2:4),'forward');
+ rB_e = threeDmulti_rot(RotMat,GNVB(:,2:4),'forward');
+ vB_e = threeDmulti_rot(RotMat,GNVB(:,5:7),'forward')+threeDmulti_rot(RotMat_dot,GNVB(:,2:4),'forward');
+
+ % check if the transformation is correct
+ plot3(vA_e(:,1),vA_e(:,2),vA_e(:,3))
+ hold on
+
+% state vector
+sa = [rA_e vA_e];
+sb = [rB_e vB_e];
+
+% calculate potential observable
+v12sst0g = potential_EBAs(KBRB(:,3),sa,sb,we0);
+
+[lon1, lat1, h1] = cart2sph(sa(:,1),sa(:,2),sa(:,3));
+[lon2, lat2, h2] = cart2sph(sb(:,1),sb(:,2),sb(:,3));
+
+%% Computation of design matrix A
+A2g = designmatrixn(lon2,lat2,h2,n_max,GM,a);
+A1g = designmatrixn(lon1,lat1,h1,n_max,GM,a);
+
+A = A2g-A1g;
+A = [A, A0];
+
+% Update N:normal equation matrix
+N = N + A'*A;
+
+% Update n0:normal equation vector
+n0 = n0+ A'*v12sst0g';
+
+fprintf('done with day %d \n', day)
+end
+
+%% Cholesky decomposition
+L = chol(N, 'lower');
+
+z0 = L\n0;
+x0 = L'\z0;
+
+%% estimate residuals
+yhat = A*x0;
+eps = v12sst0g'-yhat;
+
+toc
+
+%% writing degree variances
+[cnm, snm] = sortCoefficients(x0,n_max);
+
+kk = 0;
+
+for i=1:length(cnm)
+ for j=1:length(snm)
+ if j<=i
+ kk=kk+1;
+ xegm(kk,1)=i-1;
+ xegm(kk,2)=j-1;
+ xegm(kk,3)=cnm(i,j);
+ xegm(kk,4)=snm(i,j);
+ end
+ end
+end
+
+%% save results in a file
+FileName=[Release2_Path,'/results/EBA','_do_', num2str(n_max), '_nods_', num2str(nod),'_',datestr(now, 'dd-mmm-yyyy'),'.mat'];
+save(FileName,'xegm','eps','GM','a','n_max','e12', 'en', 'er','rA_e');
diff --git a/release2/gfr_pos.m b/release2/gfr_pos.m
new file mode 100644
index 0000000..57fa4b9
--- /dev/null
+++ b/release2/gfr_pos.m
@@ -0,0 +1,333 @@
+%% GRACE gravity field recovery (GFR) using positions
+% This is the main m-file preforming the gravity field estimation. It is
+% optimised for parallel computing using local and global parameters for
+% different arcs. All necessary parameters and data are specified in the
+% following. You have to set a number of iterations you want to perform.
+% save_vars_for_plotting.m saves all the necessary variables for plotting.
+%
+% Written by Florian Wöske, ZARM Uni Bremen, 2018-12.
+% Adapted by Neda Darbeheshti, AEI, 2019-06.
+
+clear all
+% close all;
+% clc;
+format longg;
+
+nods = 1;
+
+% check if Matlab toolbox licence is available, if not end job
+[status,errmsg] = license('checkout', 'Distrib_Computing_Toolbox');
+if status ~= 0 || nods == 1
+
+%% variable inputs
+
+% add folder name containing functions and input data (observations, ref.
+% field, initial values...)
+Release2_Path = pwd;
+addpath(genpath(Release2_Path));
+EBA_Path = fullfile(Release2_Path,'/../function_eba/');
+addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release2_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
+
+Input_Data_Path = [Release2_Path,'/../input_data'];
+addpath(genpath(Input_Data_Path));
+
+% folder where observation data are stored
+Input_Observation_Path = [Input_Data_Path, '/MockData_do20_nod12'];
+
+% Define the number of iterations for batch processor,
+ItrNo_max = 1;
+
+% order of spherical harmonic coefficients to be estimated
+lmaxcs = 20;
+
+% parameters of background gravity field
+lmaxf = 20;
+
+% number of days
+nods = 12;
+
+% number of GNV observations per time step(eg. 3A positions + 3B positions)
+no = 6;
+
+% total number of observatrion per time step
+tno=6;
+
+% Name of Run/ and simulation extention of results file name
+name_ext_str = 'test_new';
+
+% GNV Observation in ECI or ECEF (0, 1)
+obs_frame = 0;
+
+% number of observations in each arc (one day, 17280 is the standard)
+mKBR = 17280;
+% mKBR = 20;
+
+% a guide to relationship between order of spherical harmonic coefficients and number of days
+% Degree 10 1 day Degree 70 30 days
+% Degree 50 7 days Degree 120 16 days(from Gunter 2000)
+
+%% fixed inputs
+
+% n: size of gravity fields in vec format, estimated and background
+n_est = ((lmaxcs+1)^2 + lmaxcs+1)/2;
+n_back = ((lmaxf+1)^2 + lmaxf+1)/2;
+
+GM = 0.3986004415E+15;
+ae = 0.6378136300E+07;
+
+% reference for comparison
+ggm05s = importGravityField('GGM05S.gfc');
+ggm05s = ggm05s(1:n_back,:); % cut to desired length
+ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
+
+% Load a-priori/ reference gravity field
+egm96 = importGravityField('EGM96.gfc');
+egm96 = egm96(1:n_back,:); % cut to desired length
+egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
+% reformatting spherical harmonics coefficients
+field_cs = egm96_cs;
+
+field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+field = field_sc(:)';
+[outC1,outS1,nm1] = cs2vec(field_cs,false);
+field_vec0 = [nm1',outC1',outS1'];
+
+err_vec = egm96;
+err_vec(:,3:4) = ggm05s(:,3:4) - egm96(:,3:4);
+
+% initialzes Legendre polynomial calculation
+data_plm = initplm(lmaxf,2);
+
+% folder to save intermediate results
+tc = clock;
+tcs = num2str(round(tc(6)));
+Temp_Directory = [Release2_Path,'/jn', tcs, '_' name_ext_str,'/OutputArcParall'];
+mkdir(Temp_Directory); % add time tag to counter same name_ext_str
+
+%% load observations data and initial states for all days from files
+
+t0 = [2005 05 01 0 0 2];
+date0 = time2str(t0);
+
+if no == 3 % just position observations
+ GNVL1B = zeros(mKBR,4,nods);
+else % position and velocity observations
+ GNVL1B = zeros(mKBR,7,nods);
+end
+
+state0 = zeros(12,nods);
+
+% State deviation:
+x0x = zeros(12,nods);
+t0i = t0;
+
+for Mday = 1:nods
+
+ if Mday >1
+ t0i(3)=t0i(3)+1;
+ end
+
+ date = time2str(t0i);
+
+ % load GNV data GRACE A
+ data_GNV = [date, '/GNV1B_', date, '_A_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(1:6,Mday) = [r_eci_1; v_eci_1];
+
+ else
+ state0(1:6,Mday) = GNVi(1,2:7);
+ GNVL1B(:,1:4,Mday) = GNVi(1:mKBR,1:4);
+ end
+
+ % load GNV data GRACE B
+ data_GNV = [date, '/GNV1B_', date, '_B_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(7:12,Mday) = [r_eci_1; v_eci_1];
+ else
+ state0(7:12,Mday) = GNVi(1,2:7);
+ GNVL1B(:,5:7,Mday) = GNVi(1:mKBR,2:4);
+ end
+
+end
+
+% save initial state vector for comparisons
+state00 = state0;
+
+%% iteration
+
+% number of coefficients
+% 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
+
+xhat_save = zeros (12, nods, ItrNo_max+1);
+chat_save = zeros (noc, ItrNo_max+1);
+chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+
+% save tic toc times:
+t_itr_nod = zeros(ItrNo_max+1,nods+6); % +6 for clock
+t_nod = zeros(1,nods);
+
+for ItrNo = 0:ItrNo_max
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['start iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+ % batch over all days
+ parfor Mday = 1:nods
+ % for Mday = 1:nods % for breakpoints and debugging
+ tic
+ batch_processor_partitioned_pos(Mday,nods,Temp_Directory,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),GNVL1B(:,1,Mday),GNVL1B(:,:,Mday),x0x(:,Mday))
+ t_nod(Mday) = toc;
+ end
+
+ t_itr_nod(ItrNo+1,1:nods) = t_nod;
+ t_itr_nod(ItrNo+1,nods+1:end) = c; % save time when each iteration was started
+
+ % solve normal equations
+
+ % read data from arc parallel
+ invMxx = zeros(12,12,nods);
+ iMxc = zeros(12,noc,nods);
+ iMcc = zeros(noc,noc,nods);
+ iNx = zeros(12,1,nods);
+ iNc = zeros(noc,1,nods);
+ iL = zeros(noc,noc,nods);
+ iN = zeros(noc,1,nods);
+ for Mday = 1:nods
+ FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
+ File = fullfile(Temp_Directory, FileName);
+ arcstru=load(File);
+ invMxx(:,:,Mday) = arcstru.invMxx;
+ iMxc(:,:,Mday) = arcstru.Mxc;
+ iMcc(:,:,Mday) = arcstru.Mcc;
+ iNx(:,:,Mday) = arcstru.Nx;
+ iNc(:,:,Mday) = arcstru.Nc;
+ McxTinvMxx = iMxc(:,:,Mday)'*invMxx(:,:,Mday);
+ iL(:,:,Mday) = McxTinvMxx*iMxc(:,:,Mday);
+ iN(:,:,Mday) = McxTinvMxx*iNx(:,:,Mday);
+ end
+
+ % sum over the number of days
+ sumiMcc = sum(iMcc,3);
+ sumiNc = sum(iNc,3);
+ sumiL = sum(iL,3);
+ sumiN = sum(iN,3);
+ L = sumiMcc-sumiL;
+ N = sumiNc-sumiN;
+
+ % estimate global parameters
+ chat = L \ N;
+
+ % estimate local parameters
+ xhat = zeros(12,nods);
+
+ for Mday = 1:nods
+ xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
+ end
+
+ % save xhat and chat of iterations
+ xhat_save(:,:,ItrNo+1) = xhat;
+ chat_save(:,ItrNo+1) = chat;
+
+ % put coefficients in right order for plotting and output
+ ko=1;
+ for i=1:lmaxcs+1
+ for j=1:lmaxcs+1
+ if i<=j
+ ordering1(ko,1)=j-1;
+ ordering1(ko,2)=i-1;
+ ko=ko+1;
+ end
+ end
+ end
+ nocC=(noc+lmaxcs-1)/2;
+ ordering1(3:end,3)=[chat(1:lmaxcs-1)', 0 ,chat(lmaxcs:nocC)']; % put C coeffs.
+ ordering1(3+lmaxcs:end,4)=chat(nocC+1:end); % put S coeffs.
+ ordering1_cs = vec2cs([ordering1(:,1) ordering1(:,2)]',ordering1(:,3),ordering1(:,4));
+
+ chat_cs = ordering1_cs;
+ % chat_vec = ordering1; % orderwise ordering
+ [outC1,outS1,nm1] = cs2vec(chat_cs,false); %degreewise
+ chat_vec = [nm1',outC1',outS1'];
+
+ chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
+
+ % set up values for next iteration
+ state0 = state0 + xhat;
+ x0x = x0x - xhat;
+
+ % Add CS coeffs in cs format to reference field (just up to the order
+ % that was estimated)
+ field_cs(1:lmaxcs+1,1:lmaxcs+1) = field_cs(1:lmaxcs+1,1:lmaxcs+1) + chat_cs;
+ field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+ field = field_sc(:)';
+
+ [outC1,outS1,nm1] = cs2vec(field_cs,false);
+ field_vec = [nm1',outC1',outS1'];
+
+ % save dv of true gravity field minus estimated from each iteration
+ d_vec = ggm05s(1:n_back,3:4)-field_vec(:,3:4);
+ d_vec = [field_vec(:,1:2) d_vec];
+ dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['done with iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+end
+
+% estimate residuals for the last day and last iteration
+iHx = arcstru.Hx_save;
+iHc = arcstru.Hc_save;
+iygnv_save = arcstru.ygnv_save;
+
+% estimate GNVA residuals
+yhatXA=squeeze(iHx(:,1,:))*xhat(:,nods)+squeeze(iHc(:,1,:))*chat;
+eps_save(:,1)= iygnv_save(:,1)-yhatXA;
+
+yhatYA=squeeze(iHx(:,2,:))*xhat(:,nods)+squeeze(iHc(:,2,:))*chat;
+eps_save(:,2) = iygnv_save(:,2)-yhatYA;
+
+yhatZA=squeeze(iHx(:,3,:))*xhat(:,nods)+squeeze(iHc(:,3,:))*chat;
+eps_save(:,3) = iygnv_save(:,3)-yhatZA;
+
+% estimate GNVB residuals
+yhatXB=squeeze(iHx(:,4,:))*xhat(:,nods)+squeeze(iHc(:,4,:))*chat;
+eps_save(:,4)= iygnv_save(:,4)-yhatXB;
+
+yhatYB=squeeze(iHx(:,5,:))*xhat(:,nods)+squeeze(iHc(:,5,:))*chat;
+eps_save(:,5) = iygnv_save(:,5)-yhatYB;
+
+yhatZB=squeeze(iHx(:,6,:))*xhat(:,nods)+squeeze(iHc(:,6,:))*chat;
+eps_save(:,6) = iygnv_save(:,6)-yhatZB;
+
+% save results in .mat file:
+save_vars_for_plotting
+
+% remove intermediate folder
+FolderName_dd = Temp_Directory(1:end-16); % delete folder and subfolders
+rmdir(FolderName_dd,'s');
+
+else
+ % error message and end computation/ job
+ fprintf(['licence for MATLAB parallel computation toolbox not available\nComputation terminated\n'])
+ exit
+ return
+end
diff --git a/release2/gfr_pos_rrho.m b/release2/gfr_pos_rrho.m
new file mode 100644
index 0000000..1c6c8da
--- /dev/null
+++ b/release2/gfr_pos_rrho.m
@@ -0,0 +1,347 @@
+%% GRACE gravity field recovery (GFR) using range-rates and positions
+% This is the main m-file preforming the gravity field estimation. It is
+% optimised for parallel computing using local and global parameters for
+% different arcs. All necessary parameters and data are specified in the
+% following. You have to set a number of iterations you want to perform.
+% save_vars_for_plotting.m saves all the necessary variables for plotting.
+%
+% Written by Florian Wöske, ZARM Uni Bremen, 2018-12.
+% Adapted by Neda Darbeheshti, AEI, 2019-06.
+
+clear all
+% close all;
+% clc;
+format longg;
+
+nods = 1;
+
+% check if Matlab toolbox licence is available, if not end job
+[status,errmsg] = license('checkout', 'Distrib_Computing_Toolbox');
+if status ~= 0 || nods == 1
+
+%% variable inputs
+
+% add folder name containing functions and input data (observations, ref.
+% field, initial values...)
+Release2_Path = pwd;
+addpath(genpath(Release2_Path));
+EBA_Path = fullfile(Release2_Path,'/../function_eba/');
+addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release2_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
+Input_Data_Path = [Release2_Path,'/../input_data'];
+addpath(genpath(Input_Data_Path));
+
+% folder where observation data are stored
+Input_Observation_Path = [Input_Data_Path, '/MockData_do20_nod12'];
+
+% Define the number of iterations for batch processor,
+ItrNo_max = 1;
+
+% order of spherical harmonic coefficients to be estimated
+lmaxcs = 20;
+
+% parameters of background gravity field
+lmaxf = 20;
+
+% number of days
+nods = 12;
+
+% number of GNV observations per time step(eg. 3A positions + 3B positions)
+no = 6;
+
+% total number of observatrion per time step
+tno=no+1;
+
+% Name of Run/ and simulation extention of results file name
+name_ext_str = 'test_new';
+
+% GNV Observation in ECI or ECEF (0, 1)
+obs_frame = 0;
+
+% number of observations in each arc (one day, 17280 is the standard)
+mKBR = 17280;
+% mKBR = 20;
+
+% a guide to relationship between order of spherical harmonic coefficients and number of days
+% Degree 10 1 day Degree 70 30 days
+% Degree 50 7 days Degree 120 16 days(from Gunter 2000)
+
+%% fixed inputs
+
+% n: size of gravity fields in vec format, estimated and background
+n_est = ((lmaxcs+1)^2 + lmaxcs+1)/2;
+n_back = ((lmaxf+1)^2 + lmaxf+1)/2;
+
+GM = 0.3986004415E+15;
+ae = 0.6378136300E+07;
+
+% reference for comparison
+ggm05s = importGravityField('GGM05S.gfc');
+ggm05s = ggm05s(1:n_back,:); % cut to desired length
+ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
+
+% Load a-priori/ reference gravity field
+egm96 = importGravityField('EGM96.gfc');
+egm96 = egm96(1:n_back,:); % cut to desired length
+egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
+
+% reformatting spherical harmonics coefficients
+field_cs = egm96_cs;
+
+field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+field = field_sc(:)';
+[outC1,outS1,nm1] = cs2vec(field_cs,false);
+field_vec0 = [nm1',outC1',outS1'];
+
+err_vec = egm96;
+err_vec(:,3:4) = ggm05s(:,3:4) - egm96(:,3:4);
+
+% initialzes Legendre polynomial calculation
+data_plm = initplm(lmaxf,2);
+
+% folder to save intermediate results
+tc = clock;
+tcs = num2str(round(tc(6)));
+
+Temp_Directory = [Release2_Path,'/jn', tcs, '_' name_ext_str,'/OutputArcParall'];
+mkdir(Temp_Directory); % add time tag to counter same name_ext_str
+
+%% load observations data and initial states for all days from files
+
+t0 = [2005 05 01 0 0 2];
+date0 = time2str(t0);
+
+if no == 3 % just position observations
+ GNVL1B = zeros(mKBR,3,nods);
+else % position and velocity observations
+ GNVL1B = zeros(mKBR,6,nods);
+end
+
+KBRL1B = zeros(mKBR,4,nods);
+state0 = zeros(12,nods);
+
+% State deviation:
+x0x = zeros(12,nods);
+t0i = t0;
+
+for Mday = 1:nods
+
+ if Mday >1
+ t0i(3)=t0i(3)+1;
+ end
+
+ date = time2str(t0i);
+
+ % load KBR data
+ data_KBR = [date, '/KBR1B_', date, '_X_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_KBR);
+ KBRL1Bi = readKBR(data_file);
+ KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
+
+ % load GNV data GRACE A
+ data_GNV = [date, '/GNV1B_', date, '_A_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(1:6,Mday) = [r_eci_1; v_eci_1];
+
+ else
+ state0(1:6,Mday) = GNVi(1,2:7);
+ GNVL1B(:,1:3,Mday) = GNVi(1:mKBR,2:4);
+ end
+
+ % load GNV data GRACE B
+ data_GNV = [date, '/GNV1B_', date, '_B_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(7:12,Mday) = [r_eci_1; v_eci_1];
+ else
+ state0(7:12,Mday) = GNVi(1,2:7);
+ GNVL1B(:,4:6,Mday) = GNVi(1:mKBR,2:4);
+ end
+
+end
+
+% save initial state vector for comparisons
+state00 = state0;
+
+%% iteration
+
+% number of coefficients
+% 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
+
+xhat_save = zeros (12, nods, ItrNo_max+1);
+chat_save = zeros (noc, ItrNo_max+1);
+chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+
+% save tic toc times:
+t_itr_nod = zeros(ItrNo_max+1,nods+6); % +6 for clock
+t_nod = zeros(1,nods);
+
+% define diagonal weighting matrix
+WMAT0=10^(2*2)*eye(no+1);
+WMAT0(1,1)=10^(7*2);
+
+for ItrNo = 0:ItrNo_max
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['start iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+ % batch over all days
+ parfor Mday = 1:nods
+ %for Mday = 1:nods % for breakpoints and debugging
+ tic
+ batch_processor_partitioned_pos_rrho(Mday,nods,Temp_Directory,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),GNVL1B(:,:,Mday),KBRL1B(:,3,Mday),WMAT0)
+ t_nod(Mday) = toc;
+ end
+
+ t_itr_nod(ItrNo+1,1:nods) = t_nod;
+ t_itr_nod(ItrNo+1,nods+1:end) = c; % save time when each iteration was started
+
+ %% solve normal equations
+
+ % read data from arc parallel
+ invMxx = zeros(12,12,nods);
+ iMxc = zeros(12,noc,nods);
+ iMcc = zeros(noc,noc,nods);
+ iNx = zeros(12,1,nods);
+ iNc = zeros(noc,1,nods);
+ iL = zeros(noc,noc,nods);
+ iN = zeros(noc,1,nods);
+ for Mday = 1:nods
+ FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
+ File = fullfile(Temp_Directory, FileName);
+ arcstru=load(File);
+ invMxx(:,:,Mday) = arcstru.invMxx;
+ iMxc(:,:,Mday) = arcstru.Mxc;
+ iMcc(:,:,Mday) = arcstru.Mcc;
+ iNx(:,:,Mday) = arcstru.Nx;
+ iNc(:,:,Mday) = arcstru.Nc;
+ McxTinvMxx = iMxc(:,:,Mday)'*invMxx(:,:,Mday);
+ iL(:,:,Mday) = McxTinvMxx*iMxc(:,:,Mday);
+ iN(:,:,Mday) = McxTinvMxx*iNx(:,:,Mday);
+ end
+
+ % sum over the number of days
+ sumiMcc = sum(iMcc,3);
+ sumiNc = sum(iNc,3);
+ sumiL = sum(iL,3);
+ sumiN = sum(iN,3);
+ L = sumiMcc-sumiL;
+ N = sumiNc-sumiN;
+
+ % estimate global parameters
+ chat = L \ N;
+
+ % estimate local parameters
+ xhat = zeros(12,nods);
+
+ for Mday = 1:nods
+ xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
+ end
+
+ % save xhat and chat of iterations
+ xhat_save(:,:,ItrNo+1) = xhat;
+ chat_save(:,ItrNo+1) = chat;
+
+ % put coefficients in right order for plotting and output
+ ko=1;
+ for i=1:lmaxcs+1
+ for j=1:lmaxcs+1
+ if i<=j
+ ordering1(ko,1)=j-1;
+ ordering1(ko,2)=i-1;
+ ko=ko+1;
+ end
+ end
+ end
+ nocC=(noc+lmaxcs-1)/2;
+ ordering1(3:end,3)=[chat(1:lmaxcs-1)', 0 ,chat(lmaxcs:nocC)']; % put C coeffs.
+ ordering1(3+lmaxcs:end,4)=chat(nocC+1:end); % put S coeffs.
+ ordering1_cs = vec2cs([ordering1(:,1) ordering1(:,2)]',ordering1(:,3),ordering1(:,4));
+
+ chat_cs = ordering1_cs;
+ % chat_vec = ordering1; % orderwise ordering
+ [outC1,outS1,nm1] = cs2vec(chat_cs,false); % degreewise
+ chat_vec = [nm1',outC1',outS1'];
+
+ chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
+
+ % set up values for next iteration
+ state0 = state0 + xhat;
+ x0x = x0x - xhat;
+
+ % Add CS coeffs in cs format to reference field (just up to the order
+ % that was estimated)
+ field_cs(1:lmaxcs+1,1:lmaxcs+1) = field_cs(1:lmaxcs+1,1:lmaxcs+1) + chat_cs;
+ field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+ field = field_sc(:)';
+
+ [outC1,outS1,nm1] = cs2vec(field_cs,false);
+ field_vec = [nm1',outC1',outS1'];
+
+ % save dv of true gravity field minus estimated from each iteration
+ d_vec = ggm05s(1:n_back,3:4)-field_vec(:,3:4);
+ d_vec = [field_vec(:,1:2) d_vec];
+ dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['done with iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+end
+ % estimate residuals for the last day and last iteration
+iHx = arcstru.Hx_save;
+iHc = arcstru.Hc_save;
+iykbr_save = arcstru.ykbr_save;
+iygnv_save = arcstru.ygnv_save;
+
+ % estimate range rate residuals
+ yhatK=squeeze(iHx(:,1,:))*xhat(:,nods)+squeeze(iHc(:,1,:))*chat;
+ eps_save(:,1) = iykbr_save-yhatK;
+ % estimate GNVA residuals
+ yhatXA=squeeze(iHx(:,2,:))*xhat(:,nods)+squeeze(iHc(:,2,:))*chat;
+ eps_save(:,2)= iygnv_save(:,1)-yhatXA;
+
+ yhatYA=squeeze(iHx(:,3,:))*xhat(:,nods)+squeeze(iHc(:,3,:))*chat;
+ eps_save(:,3) = iygnv_save(:,2)-yhatYA;
+
+ yhatZA=squeeze(iHx(:,4,:))*xhat(:,nods)+squeeze(iHc(:,4,:))*chat;
+ eps_save(:,4) = iygnv_save(:,3)-yhatZA;
+ % estimate GNVB residuals
+ yhatXB=squeeze(iHx(:,5,:))*xhat(:,nods)+squeeze(iHc(:,5,:))*chat;
+ eps_save(:,5)= iygnv_save(:,4)-yhatXB;
+
+ yhatYB=squeeze(iHx(:,6,:))*xhat(:,nods)+squeeze(iHc(:,6,:))*chat;
+ eps_save(:,6) = iygnv_save(:,5)-yhatYB;
+
+ yhatZB=squeeze(iHx(:,7,:))*xhat(:,nods)+squeeze(iHc(:,7,:))*chat;
+ eps_save(:,7) = iygnv_save(:,6)-yhatZB;
+
+% save results in .mat file:
+save_vars_for_plotting
+
+% remove intermediate folder
+FolderName_dd = Temp_Directory(1:end-16); % delete folder and subfolders
+rmdir(FolderName_dd,'s');
+
+else
+ % error message and end computation/ job
+ fprintf(['licence for MATLAB parallel computation toolbox not available\nComputation terminated\n'])
+ exit
+ return
+end
diff --git a/release2/gfr_rrho.m b/release2/gfr_rrho.m
new file mode 100644
index 0000000..0137971
--- /dev/null
+++ b/release2/gfr_rrho.m
@@ -0,0 +1,317 @@
+%% GRACE gravity field recovery (GFR) using range-rates
+% This is the main m-file preforming the gravity field estimation. It is
+% optimised for parallel computing using local and global parameters for
+% different arcs. All necessary parameters and data are specified in the
+% following. You have to set a number of iterations you want to perform.
+% save_vars_for_plotting.m saves all the necessary variables for plotting.
+%
+% Written by Florian Wöske, ZARM Uni Bremen, 2018-12.
+% Adapted by Neda Darbeheshti, AEI, 2019-06.
+
+clear all
+% close all;
+% clc;
+format longg;
+
+nods = 1;
+
+% check if Matlab toolbox licence is available, if not end job
+[status,errmsg] = license('checkout', 'Distrib_Computing_Toolbox');
+if status ~= 0 || nods == 1
+
+%% variable inputs
+
+% add folder name containing functions and input data (observations, ref.
+% field, initial values...)
+Release2_Path = pwd;
+addpath(genpath(Release2_Path));
+EBA_Path = fullfile(Release2_Path,'/../function_eba/');
+addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release2_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
+Input_Data_Path = [Release2_Path,'/../input_data'];
+addpath(genpath(Input_Data_Path));
+
+% folder where observation data are stored
+Input_Observation_Path = [Input_Data_Path, '/MockData_do20_nod12'];
+
+% Define the number of iterations for batch processor,
+ItrNo_max = 1;
+
+% order of spherical harmonic coefficients to be estimated
+lmaxcs = 20;
+% parameters of background gravity field
+lmaxf = 20;
+% number of days
+nods = 12;
+% number of GNV observations per time step(eg. 3A positions + 3B positions)
+no = 0;
+% total number of observatrion per time step
+tno=1;
+% Name of Run/ and simulation extention of results file name
+name_ext_str = 'test_new';
+
+% GNV Observation in ECI or ECEF (0, 1)
+obs_frame = 0;
+
+% number of observations in each arc (one day, 17280 is the standard)
+mKBR = 17280;
+%mKBR = 20;
+
+% a guide to relationship between order of spherical harmonic coefficients and number of days
+% Degree 10 1 day Degree 70 30 days
+% Degree 50 7 days Degree 120 16 days(from Gunter 2000)
+
+%% fixed inputs
+
+% n: size of gravity fields in vec format, estimated and background
+n_est = ((lmaxcs+1)^2 + lmaxcs+1)/2;
+n_back = ((lmaxf+1)^2 + lmaxf+1)/2;
+
+GM = 0.3986004415E+15;
+ae = 0.6378136300E+07;
+
+% reference for comparison
+ggm05s = importGravityField('GGM05S.gfc');
+ggm05s = ggm05s(1:n_back,:); % cut to desired length
+ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
+
+% Load a-priori/ reference gravity field
+egm96 = importGravityField('EGM96.gfc');
+egm96 = egm96(1:n_back,:); % cut to desired length
+egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
+% reformatting spherical harmonics coefficients
+field_cs = egm96_cs;
+
+field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+field = field_sc(:)';
+[outC1,outS1,nm1] = cs2vec(field_cs,false);
+field_vec0 = [nm1',outC1',outS1'];
+
+err_vec = egm96;
+err_vec(:,3:4) = ggm05s(:,3:4) - egm96(:,3:4);
+
+% initialzes Legendre polynomial calculation
+data_plm = initplm(lmaxf,2);
+
+% folder to save intermediate results
+tc = clock;
+tcs = num2str(round(tc(6)));
+Temp_Directory = [Release2_Path,'/jn', tcs, '_' name_ext_str,'/OutputArcParall'];
+mkdir(Temp_Directory); % add time tag to counter same name_ext_str
+
+%% load observations data and initial states for all days from files
+
+t0 = [2005 05 01 0 0 2];
+date0 = time2str(t0);
+
+if no == 3 % just position observations
+ GNVL1B = zeros(mKBR,3,nods);
+else % position and velocity observations
+ GNVL1B = zeros(mKBR,6,nods);
+end
+
+KBRL1B = zeros(mKBR,4,nods);
+state0 = zeros(12,nods);
+
+% State deviation:
+x0x = zeros(12,nods);
+t0i = t0;
+
+for Mday = 1:nods
+
+ if Mday >1
+ t0i(3)=t0i(3)+1;
+ end
+
+ date = time2str(t0i);
+
+ % load KBR data
+ data_KBR = [date, '/KBR1B_', date, '_X_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_KBR);
+ KBRL1Bi = readKBR(data_file);
+ KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
+
+ % load GNV data GRACE A
+ data_GNV = [date, '/GNV1B_', date, '_A_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(1:6,Mday) = [r_eci_1; v_eci_1];
+
+ else
+ state0(1:6,Mday) = GNVi(1,2:7);
+ end
+
+ % load GNV data GRACE B
+ data_GNV = [date, '/GNV1B_', date, '_B_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(7:12,Mday) = [r_eci_1; v_eci_1];
+ else
+ state0(7:12,Mday) = GNVi(1,2:7);
+ end
+end
+
+% save initial state vector for comparisons
+state00 = state0;
+
+%% iteration
+
+% number of coefficients
+% 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
+
+% SH deviation:
+% x0c = zeros(noc,1);
+xhat_save = zeros (12, nods, ItrNo_max+1);
+chat_save = zeros (noc, ItrNo_max+1);
+chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+
+% save tic toc times:
+t_itr_nod = zeros(ItrNo_max+1,nods+6); % +6 for clock
+t_nod = zeros(1,nods);
+
+for ItrNo = 0:ItrNo_max
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['start iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+ % batch over all days
+ parfor Mday = 1:nods
+ % for Mday = 1:nods % for breakpoints and debugging
+ tic
+ batch_processor_partitioned_rrho(Mday,nods,Temp_Directory,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
+ t_nod(Mday) = toc;
+ end
+
+ t_itr_nod(ItrNo+1,1:nods) = t_nod;
+ t_itr_nod(ItrNo+1,nods+1:end) = c; % save time when each iteration was started
+
+ %% solve normal equations
+
+ % read data from arc parallel
+ invMxx = zeros(12,12,nods);
+ iMxc = zeros(12,noc,nods);
+ iMcc = zeros(noc,noc,nods);
+ iNx = zeros(12,1,nods);
+ iNc = zeros(noc,1,nods);
+ iL = zeros(noc,noc,nods);
+ iN = zeros(noc,1,nods);
+ for Mday = 1:nods
+ FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
+ File = fullfile(Temp_Directory, FileName);
+ arcstru=load(File);
+ invMxx(:,:,Mday) = arcstru.invMxx;
+ iMxc(:,:,Mday) = arcstru.Mxc;
+ iMcc(:,:,Mday) = arcstru.Mcc;
+ iNx(:,:,Mday) = arcstru.Nx;
+ iNc(:,:,Mday) = arcstru.Nc;
+ McxTinvMxx = iMxc(:,:,Mday)'*invMxx(:,:,Mday);
+ iL(:,:,Mday) = McxTinvMxx*iMxc(:,:,Mday);
+ iN(:,:,Mday) = McxTinvMxx*iNx(:,:,Mday);
+ end
+
+ % sum over the number of days
+ sumiMcc = sum(iMcc,3);
+ sumiNc = sum(iNc,3);
+ sumiL = sum(iL,3);
+ sumiN = sum(iN,3);
+ L = sumiMcc-sumiL;
+ N = sumiNc-sumiN;
+
+ % estimate global parameters
+ chat = L \ N;
+
+ % estimate local parameters
+ xhat = zeros(12,nods);
+
+ for Mday = 1:nods
+ xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
+ end
+
+ % save xhat and chat of iterations
+ xhat_save(:,:,ItrNo+1) = xhat;
+ chat_save(:,ItrNo+1) = chat;
+
+ % put coefficients in right order for plotting and output
+ ko=1;
+ for i=1:lmaxcs+1
+ for j=1:lmaxcs+1
+ if i<=j
+ ordering1(ko,1)=j-1;
+ ordering1(ko,2)=i-1;
+ ko=ko+1;
+ end
+ end
+ end
+ nocC=(noc+lmaxcs-1)/2;
+ ordering1(3:end,3)=[chat(1:lmaxcs-1)', 0 ,chat(lmaxcs:nocC)']; % put C coeffs.
+ ordering1(3+lmaxcs:end,4)=chat(nocC+1:end); % put S coeffs.
+ ordering1_cs = vec2cs([ordering1(:,1) ordering1(:,2)]',ordering1(:,3),ordering1(:,4));
+
+ chat_cs = ordering1_cs;
+ % chat_vec = ordering1; % orderwise ordering
+ [outC1,outS1,nm1] = cs2vec(chat_cs,false); %degreewise
+ chat_vec = [nm1',outC1',outS1'];
+
+ chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
+
+ % set up values for next iteration
+ state0 = state0 + xhat;
+ x0x = x0x - xhat;
+
+ % Add CS coeffs in cs format to reference field (just up to the order
+ % that was estimated)
+ field_cs(1:lmaxcs+1,1:lmaxcs+1) = field_cs(1:lmaxcs+1,1:lmaxcs+1) + chat_cs;
+ field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+ field = field_sc(:)';
+ % x0c=x0c-chat;
+
+ [outC1,outS1,nm1] = cs2vec(field_cs,false);
+ field_vec = [nm1',outC1',outS1'];
+
+ % save dv of true gravity field minus estimated from each iteration
+ d_vec = ggm05s(1:n_back,3:4)-field_vec(:,3:4);
+ d_vec = [field_vec(:,1:2) d_vec];
+ dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['done with iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+end
+ % estimate residuals for the last day and last iteration
+iHx = arcstru.Hx_save;
+iHc = arcstru.Hc_save;
+iykbr_save = arcstru.ykbr_save;
+
+% estimate range rate residuals
+yhatK=iHx*xhat(:,nods)+iHc*chat;
+eps_save(:,1) = iykbr_save-yhatK;
+
+% save results in .mat file:
+save_vars_for_plotting
+
+% remove intermediate folder
+FolderName_dd = Temp_Directory(1:end-16); % delete folder and subfolders
+rmdir(FolderName_dd,'s');
+
+else
+ % error message and end computation/ job
+ fprintf(['licence for MATLAB parallel computation toolbox not available\nComputation terminated\n'])
+ exit
+ return
+end
diff --git a/release2/gfr_rrho_tickonov.m b/release2/gfr_rrho_tickonov.m
new file mode 100644
index 0000000..353263f
--- /dev/null
+++ b/release2/gfr_rrho_tickonov.m
@@ -0,0 +1,322 @@
+%% GRACE gravity field recovery (GFR) using range-rates with Tickonov regularization
+% This is the main m-file preforming the gravity field estimation. It is
+% optimised for parallel computing using local and global parameters for
+% different arcs. All necessary parameters and data are specified in the
+% following. You have to set a number of iterations you want to perform.
+% save_vars_for_plotting.m saves all the necessary variables for plotting.
+%
+% Written by Florian Wöske, ZARM Uni Bremen, 2018-12.
+% Adapted by Neda Darbeheshti, AEI, 2019-06.
+
+clear all
+% close all;
+% clc;
+format longg;
+
+nods = 1;
+
+% check if Matlab toolbox licence is available, if not end job
+[status,errmsg] = license('checkout', 'Distrib_Computing_Toolbox');
+if status ~= 0 || nods == 1
+%% variable inputs
+
+% add folder name containing functions and input data (observations, ref.
+% field, initial values...)
+Release2_Path = pwd;
+addpath(genpath(Release2_Path));
+EBA_Path = fullfile(Release2_Path,'/../function_eba/');
+addpath(genpath(EBA_Path));
+GFR_Path = fullfile(Release2_Path,'/../function_gfr/');
+addpath(genpath(GFR_Path));
+Input_Data_Path = [Release2_Path,'/../input_data'];
+addpath(genpath(Input_Data_Path));
+
+% folder where observation data are stored
+Input_Observation_Path = [Input_Data_Path, '/MockData_do20_nod12'];
+
+% Define the number of iterations for batch processor,
+ItrNo_max = 1;
+
+% order of spherical harmonic coefficients to be estimated
+lmaxcs = 20;
+
+% parameters of background gravity field
+lmaxf = 20;
+
+% number of days
+nods = 12;
+
+% number of GNV observations per time step(eg. 3A positions + 3B positions)
+no = 0;
+
+% total number of observatrion per time step
+tno=1;
+% Name of Run/ and simulation extention of results file name
+name_ext_str = 'test_new';
+
+% GNV Observation in ECI or ECEF (0, 1)
+obs_frame = 0;
+
+% number of observations in each arc (one day, 17280 is the standard)
+mKBR = 17280;
+% mKBR = 20;
+
+% a guide to relationship between order of spherical harmonic coefficients and number of days
+% Degree 10 1 day Degree 70 30 days
+% Degree 50 7 days Degree 120 16 days(from Gunter 2000)
+
+%% fixed inputs
+
+% n: size of gravity fields in vec format, estimated and background
+n_est = ((lmaxcs+1)^2 + lmaxcs+1)/2;
+n_back = ((lmaxf+1)^2 + lmaxf+1)/2;
+
+GM = 0.3986004415E+15;
+ae = 0.6378136300E+07;
+
+% reference for comparison
+ggm05s = importGravityField('GGM05S.gfc');
+ggm05s = ggm05s(1:n_back,:); % cut to desired length
+ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
+
+% Load a-priori/ reference gravity field
+egm96 = importGravityField('EGM96.gfc');
+egm96 = egm96(1:n_back,:); % cut to desired length
+egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));
+
+% reformatting spherical harmonics coefficients
+field_cs = egm96_cs;
+
+field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+field = field_sc(:)';
+[outC1,outS1,nm1] = cs2vec(field_cs,false);
+field_vec0 = [nm1',outC1',outS1'];
+
+err_vec = egm96;
+err_vec(:,3:4) = ggm05s(:,3:4) - egm96(:,3:4);
+
+% initialzes Legendre polynomial calculation
+data_plm = initplm(lmaxf,2);
+
+% folder to save intermediate results
+tc = clock;
+tcs = num2str(round(tc(6)));
+Temp_Directory = [Release2_Path,'/jn', tcs, '_' name_ext_str,'/OutputArcParall'];
+mkdir(Temp_Directory); % add time tag to counter same name_ext_str
+
+%% load observations data and initial states for all days from files
+
+t0 = [2005 05 01 0 0 2];
+date0 = time2str(t0);
+
+if no == 3 % just position observations
+ GNVL1B = zeros(mKBR,3,nods);
+else % position and velocity observations
+ GNVL1B = zeros(mKBR,6,nods);
+end
+
+KBRL1B = zeros(mKBR,4,nods);
+state0 = zeros(12,nods);
+
+% State deviation:
+x0x = zeros(12,nods);
+t0i = t0;
+
+for Mday = 1:nods
+
+ if Mday >1
+ t0i(3)=t0i(3)+1;
+ end
+
+ date = time2str(t0i);
+
+ % load KBR data
+ data_KBR = [date, '/KBR1B_', date, '_X_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_KBR);
+ KBRL1Bi = readKBR(data_file);
+ KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
+
+ % load GNV data GRACE A
+ data_GNV = [date, '/GNV1B_', date, '_A_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(1:6,Mday) = [r_eci_1; v_eci_1];
+
+ else
+ state0(1:6,Mday) = GNVi(1,2:7);
+ end
+
+ % load GNV data GRACE B
+ data_GNV = [date, '/GNV1B_', date, '_B_02.asc'];
+ data_file = fullfile(Input_Observation_Path, data_GNV);
+ GNVi = readGNV(data_file);
+ if obs_frame == 1 % trafo to eci frame
+ RotMat_1 = Ri2e(data_GNV(1,1,Mday));
+ RotMat_dot_1 = Ri2e_dot(data_GNV(1,1,Mday));
+ v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
+ r_eci_1 = RotMat_1'*data_GNV(1,2:4,Mday)';
+ state0(7:12,Mday) = [r_eci_1; v_eci_1];
+ else
+ state0(7:12,Mday) = GNVi(1,2:7);
+ end
+
+end
+
+% save initial state vector for comparisons
+state00 = state0;
+
+%% iteration
+
+% number of coefficients
+% 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
+
+% SH deviation:
+% x0c = zeros(noc,1);
+xhat_save = zeros (12, nods, ItrNo_max+1);
+chat_save = zeros (noc, ItrNo_max+1);
+chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+dv_save = zeros (ItrNo_max+1,lmaxcs-1);
+
+% save tic toc times:
+t_itr_nod = zeros(ItrNo_max+1,nods+6); % +6 for clock
+t_nod = zeros(1,nods);
+
+for ItrNo = 0:ItrNo_max
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['start iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+ % batch over all days
+ parfor Mday = 1:nods
+ %for Mday = 1:nods % for breakpoints and debugging
+ tic
+ batch_processor_partitioned_rrho_tickonov(Mday,nods,Temp_Directory,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
+ t_nod(Mday) = toc;
+ end
+
+ t_itr_nod(ItrNo+1,1:nods) = t_nod;
+ t_itr_nod(ItrNo+1,nods+1:end) = c; % save time when each iteration was started
+
+ %% solve normal equations
+
+ % read data from arc parallel
+ invMxx = zeros(12,12,nods);
+ iMxc = zeros(12,noc,nods);
+ iMcc = zeros(noc,noc,nods);
+ iNx = zeros(12,1,nods);
+ iNc = zeros(noc,1,nods);
+ iL = zeros(noc,noc,nods);
+ iN = zeros(noc,1,nods);
+ for Mday = 1:nods
+ FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
+ File = fullfile(Temp_Directory, FileName);
+ arcstru=load(File);
+ invMxx(:,:,Mday) = arcstru.invMxx;
+ iMxc(:,:,Mday) = arcstru.Mxc;
+ iMcc(:,:,Mday) = arcstru.Mcc;
+ iNx(:,:,Mday) = arcstru.Nx;
+ iNc(:,:,Mday) = arcstru.Nc;
+ McxTinvMxx = iMxc(:,:,Mday)'*invMxx(:,:,Mday);
+ iL(:,:,Mday) = McxTinvMxx*iMxc(:,:,Mday);
+ iN(:,:,Mday) = McxTinvMxx*iNx(:,:,Mday);
+ end
+
+ % sum over the number of days
+ sumiMcc = sum(iMcc,3);
+ sumiNc = sum(iNc,3);
+ sumiL = sum(iL,3);
+ sumiN = sum(iN,3);
+ L = sumiMcc-sumiL;
+ N = sumiNc-sumiN;
+
+ % estimate global parameters
+ chat = L \ N;
+
+ % estimate local parameters
+ xhat = zeros(12,nods);
+
+ for Mday = 1:nods
+ xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
+ end
+
+ % save xhat and chat of iterations
+ xhat_save(:,:,ItrNo+1) = xhat;
+ chat_save(:,ItrNo+1) = chat;
+
+ % put coefficients in right order for plotting and output
+ ko=1;
+ for i=1:lmaxcs+1
+ for j=1:lmaxcs+1
+ if i<=j
+ ordering1(ko,1)=j-1;
+ ordering1(ko,2)=i-1;
+ ko=ko+1;
+ end
+ end
+ end
+ nocC=(noc+lmaxcs-1)/2;
+ ordering1(3:end,3)=[chat(1:lmaxcs-1)', 0 ,chat(lmaxcs:nocC)']; % put C coeffs.
+ ordering1(3+lmaxcs:end,4)=chat(nocC+1:end); % put S coeffs.
+ ordering1_cs = vec2cs([ordering1(:,1) ordering1(:,2)]',ordering1(:,3),ordering1(:,4));
+
+ chat_cs = ordering1_cs;
+ % chat_vec = ordering1; % orderwise ordering
+ [outC1,outS1,nm1] = cs2vec(chat_cs,false); % degreewise
+ chat_vec = [nm1',outC1',outS1'];
+
+ chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
+
+ % set up values for next iteration
+ state0 = state0 + xhat;
+ x0x = x0x - xhat;
+
+ % Add CS coeffs in cs format to reference field (just up to the order
+ % that was estimated)
+ field_cs(1:lmaxcs+1,1:lmaxcs+1) = field_cs(1:lmaxcs+1,1:lmaxcs+1) + chat_cs;
+ field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
+ field = field_sc(:)';
+ % x0c=x0c-chat;
+
+ [outC1,outS1,nm1] = cs2vec(field_cs,false);
+ field_vec = [nm1',outC1',outS1'];
+
+ % save dv of true gravity field minus estimated from each iteration
+ d_vec = ggm05s(1:n_back,3:4)-field_vec(:,3:4);
+ d_vec = [field_vec(:,1:2) d_vec];
+ dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);
+
+ c = clock;
+ [~,cs]=time2str(c);
+ fprintf(['done with iteration ',num2str(ItrNo) ,', time: ', cs,'\n'])
+
+end
+ % estimate residuals for the last day and last iteration
+iHx = arcstru.Hx_save;
+iHc = arcstru.Hc_save;
+iykbr_save = arcstru.ykbr_save;
+
+% estimate range rate residuals
+yhatK=iHx*xhat(:,nods)+iHc*chat;
+eps_save(:,1) = iykbr_save-yhatK;
+
+% save results in .mat file:
+save_vars_for_plotting
+
+% remove intermediate folder
+FolderName_dd = Temp_Directory(1:end-16); % delete folder and subfolders
+rmdir(FolderName_dd,'s');
+
+else
+ % error message and end computation/ job
+ fprintf(['licence for MATLAB parallel computation toolbox not available\nComputation terminated\n'])
+ exit
+ return
+end
diff --git a/release2/results/.gitkeep b/release2/results/.gitkeep
new file mode 100644
index 0000000..e69de29
--
GitLab