Skip to content
Snippets Groups Projects
Commit 51ae6426 authored by Axel Schnitger's avatar Axel Schnitger
Browse files

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.
parent 138f3ee4
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 2035 additions and 78 deletions
.gitignore
+ 2
− 0
View file @ 51ae6426
......@@ -228,3 +228,5 @@ TSWLatexianTemp*
# End of https://www.gitignore.io/api/latex
variational.pdf
input_data/
results/
README.md
+ 125
− 16
View file @ 51ae6426
......@@ -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)
- [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 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)
......
function_eba/calculatePnm.m 0 → 100644
+ 41
− 0
View file @ 51ae6426
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
function_eba/designmatrixn.m 0 → 100644
+ 41
− 0
View file @ 51ae6426
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
function_eba/potential_EBAs.m 0 → 100644
+ 41
− 0
View file @ 51ae6426
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
function_eba/sortCoefficients.m 0 → 100644
+ 38
− 0
View file @ 51ae6426
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
function_eba/threeDmulti_rot.m 0 → 100644
+ 45
− 0
View file @ 51ae6426
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
function_gfr/batch_processor_partitioned_pos.m 0 → 100644
+ 184
− 0
View file @ 51ae6426
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
function_gfr/batch_processor_partitioned_pos_rrho.m 0 → 100644
+ 201
− 0
View file @ 51ae6426
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
function_gfr/batch_processor_partitioned_rrho.m 0 → 100644
+ 179
− 0
View file @ 51ae6426
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
function_gfr/batch_processor_partitioned_rrho_tickonov.m 0 → 100644
+ 182
− 0
View file @ 51ae6426
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
function_gfr/hmat_pos.m 0 → 100644
+ 42
− 0
View file @ 51ae6426
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
function_gfr/save_vars_for_plotting.m 0 → 100644
+ 17
− 0
View file @ 51ae6426
%% 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);
get_input_data.sh
+ 17
− 14
View file @ 51ae6426
#!/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
continue_gfr_par_from_iteration.m → release1/continue_gfr_par_from_iteration.m
+ 15
− 17
View file @ 51ae6426
......@@ -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,:));
......
gfr_parallel.m → release1/gfr_parallel.m
+ 32
− 29
View file @ 51ae6426
......@@ -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');
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,22 +106,21 @@ 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
......@@ -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;
......@@ -232,7 +234,7 @@ for ItrNo = 0: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,:));
......@@ -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
......
save_vars2continue_itr.m → release1/save_vars2continue_itr.m
+ 2
− 2
View file @ 51ae6426
......@@ -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');
release2/gfr_eba.m 0 → 100644
+ 151
− 0
View file @ 51ae6426
%% 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');
release2/gfr_pos.m 0 → 100644
+ 333
− 0
View file @ 51ae6426
%% 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
release2/gfr_pos_rrho.m 0 → 100644
+ 347
− 0
View file @ 51ae6426
%% 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint