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

Merge branch 'release2' into 'master' 

Release2

See merge request geoq/gracetools!1
parents 138f3ee4 51ae6426
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
release2/gfr_rrho.m 0 → 100644
+ 317
− 0
View file @ 84c33e44
%% 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
release2/gfr_rrho_tickonov.m 0 → 100644
+ 322
− 0
View file @ 84c33e44
%% 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
release2/results/.gitkeep 0 → 100644
+ 0
− 0
View file @ 84c33e44
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint