diff --git a/GFR_Iteration_10_do_70_nods_12_2005-5-12_011.mat b/GFR_Iteration_10_do_70_nods_12_2005-5-12_011.mat
deleted file mode 100644
index 79d71b26c285e56cc4a7f7d2790c9e818aa8046f..0000000000000000000000000000000000000000
Binary files a/GFR_Iteration_10_do_70_nods_12_2005-5-12_011.mat and /dev/null differ
diff --git a/gfr_parallel4BMDC1.m b/gfr_parallel4BMDC1.m
deleted file mode 100644
index 7cab35ee3b024e6452a7c8948aeb2007d0f6eb9d..0000000000000000000000000000000000000000
--- a/gfr_parallel4BMDC1.m
+++ /dev/null
@@ -1,309 +0,0 @@
-%% GRACE gravity field recovery (GFR)
-%
-% GFR_PARALLEL 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. With
-% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
-% iteration. Therefore you need to save all necessary variables, e.g. use the
-% file save_vars2continue_itr.m to save them in one .mat file. In the end some
-% example plots are given you may change or extend.
-%
-% Written by Florian W??ske, ZARM Uni Bremen, and
-
-% Neda Darbeheshti, AEI Hannover, 2018-07.
-
-%%
-format longg;
-
-%% variable inputs
-
-% add folder with all functions
-FolderName = pwd;
-addpath(genpath(FolderName));
-
-% folder to save intermediate results
-mkdir('results/OutputArcParall');
-mkdir('results/OutputNormalEqu');
-
-% Define the number of iterations for batch processor,
-ItrNo_max = 10;
-
-% order of spherical harmonic coefficients to be estimated
-lmaxcs = 70;
-
-% parameters of background gravity field
-lmaxf = 70;
-
-% number of days
-nods = 12;
-
-% folder where observation data are stored
-data_folder = fullfile(FolderName,'input_data','MDC_basic_challenge1');
-addpath(genpath(data_folder));
-
-% number of observations in each arc (one day, 17280 is the standard)
-mKBR = 17280;
-
-%% 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 and Earth radius of gravity field
-GM = 0.3986004415E+15;
-ae = 0.6378136300E+07;
-
-% reference for comparison
-ggm05s = readGFC('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 = 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(egm96_cs(1:lmaxf+1,1:lmaxf+1),0);
-field = field_sc(:)';
-[outC1,outS1,nm1] = cs2vec(egm96_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);
-
-%% load observations data and initial states for all days from files
-
-t0 = [2005 05 01 0 0 2];
-date0 = time2str(t0);
-
-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);
- date(6)=[]; % wrong date format in Nedas MDC mock data
- % load KBR data
- data_KBR = ['KBR1B_', date, '_X','_MDC-BAS1_01.asc'];
- data_file = fullfile(data_folder, data_KBR);
- KBRL1Bi = readKBR(data_file);
- KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
-end
-% Load initial states in inertial frame from POD
-%GNVi=load('POD_iter_4_do_70_nods_30_2005-5-30_GA_BMDC');
-GNVi=load('POD_iter_4_do_70_nods_12_2005-5-12_A');
-state0(1:6,1:Mday)=GNVi.state_est(1:6,1:Mday);
-%GNVi=load('POD_iter_4_do_70_nods_30_2005-5-30_GB_BMDC');
-GNVi=load('POD_iter_4_do_70_nods_12_2005-5-12_B');
-state0(7:12,1:Mday)=GNVi.state_est(1:6,1:Mday);
- % 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);
-eps_save = zeros(mKBR,nods,ItrNo_max+1);
-chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-
-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))
- toc
- end
-
- %% 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);
- iHx = zeros(mKBR,12,nods);
- iHc = zeros(mKBR,noc,nods);
- irhodot_deviation = zeros(mKBR,nods);
- for Mday = 1:nods
- FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
- File = fullfile(FolderName,'results/OutputArcParall', 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);
- iHx(:,:,Mday) = arcstru.Hx_save;
- iHc(:,:,Mday) = arcstru.Hc_save;
- irhodot_deviation(:,Mday) = arcstru.rhodot_deviation;
- 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);
- yhat = zeros(mKBR,nods);
- for Mday = 1:nods
- xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
- % estimate range rate residuals
- yhat(:,Mday)=iHx(:,:,Mday)*xhat(:,Mday)+iHc(:,:,Mday)*chat;
- eps = irhodot_deviation(:,Mday)-yhat(:,Mday);
- eps_save(:,Mday,ItrNo+1) = eps;
-% figure;plot(eps)
-% xlabel('Observation Number');
-% ylabel('residuals [m/s]');
-% 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;
-
- % 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);
-% 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);
-
- % 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);
-
- fprintf('done with iteration %d \n', ItrNo)
-
-end
-
-% plot chat from all iterations:
-figure;
-for i = 1:ItrNo_max+1
- str = ['iter ',num2str(i-1)];
- semilogy(2:lmaxcs,chat_dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
- hold on
-end
-fs = 12;
-set(gcf,'PaperPositionMode','auto')
-set(gca,'FontSize',fs);
-xlabel('Spherical harmonic degree n')
-ylabel('Square root of degree variances')
-title('chat from each iteration')
-grid on
-legend
-
-figure;
-dv_geoidn(ggm05s,lmaxcs);
-hold on
-% dv_geoidn(d_vec_i,lmaxcs);
-
-[outC1,outS1,nm1] = cs2vec(field_cs,false);
-final_field_vec = [nm1',outC1',outS1'];
-dv_geoidn(final_field_vec,lmaxcs);
-
-d_vec = ggm05s(1:n_back,3:4)-final_field_vec(:,3:4);
-d_vec = [final_field_vec(:,1:2) d_vec];
-dv_geoidn(d_vec,lmaxcs);
-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
-dv_geoidn(field_vec0,lmaxcs);
-dv_geoidn(err_vec,lmaxcs);
-legend('ggm05s','ggm05s + error field', 'error field')
-for i = 1:ItrNo_max+1
- str = ['error iter ',num2str(i-1)];
- semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
-end
-
-% plot range rate residuals in time and frequency for last iteration
-figure;
-plot(eps)
-xlabel('Observation Number');
-ylabel('residuals [m/s]');
-title(['Range-rate residuals, Itr',num2str(ItrNo,'%.2d'),',Day',num2str(Mday,'%.2d')])
-ASD_pw(eps,.2);
-%ASD(eps,.2);
-
diff --git a/gfr_parallel4BMDC1_earthframe.m b/gfr_parallel4BMDC1_earthframe.m
deleted file mode 100644
index c73e36c0bda93a0930d3a33235cc3bad3d956de7..0000000000000000000000000000000000000000
--- a/gfr_parallel4BMDC1_earthframe.m
+++ /dev/null
@@ -1,282 +0,0 @@
-%% GRACE gravity field recovery (GFR)
-%
-% GFR_PARALLEL 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. With
-% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
-% iteration. Therefore you need to save all necessary variables, e.g. use the
-% file save_vars2continue_itr.m to save them in one .mat file. In the end some
-% example plots are given you may change or extend.
-%
-% Written by Florian W??ske, ZARM Uni Bremen, and
-
-% Neda Darbeheshti, AEI Hannover, 2018-07.
-
-%%
-format longg;
-
-%% variable inputs
-
-% add folder with all functions
-FolderName = pwd;
-addpath(genpath(FolderName));
-
-% folder to save intermediate results
-mkdir('results/OutputArcParall');
-mkdir('results/OutputNormalEqu');
-
-% Define the number of iterations for batch processor,
-ItrNo_max = 9;
-
-% order of spherical harmonic coefficients to be estimated
-lmaxcs = 70;
-
-% parameters of background gravity field
-lmaxf = 70;
-
-% number of days
-nods = 13;
-
-% folder where observation data are stored
-data_folder = fullfile(FolderName,'input_data','MDC_basic_challenge1');
-addpath(genpath(data_folder));
-
-% number of observations in each arc (one day, 17280 is the standard)
-mKBR = 17280;
-
-%% 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 and Earth radius of gravity field
-GM = 0.3986004415E+15;
-ae = 0.6378136300E+07;
-
-% reference for comparison
-ggm05s = readGFC('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 = 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(egm96_cs(1:lmaxf+1,1:lmaxf+1),0);
-field = field_sc(:)';
-[outC1,outS1,nm1] = cs2vec(egm96_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);
-
-%% load observations data and initial states for all days from files
-
-t0 = [2005 05 01 0 0 2];
-date0 = time2str(t0);
-
-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);
- date(6)=[]; % wrong date format in Nedas MDC mock data
- % load KBR data
- data_KBR = ['KBR1B_', date, '_X','_MDC-BAS1_01.asc'];
- data_file = fullfile(data_folder, data_KBR);
- KBRL1Bi = readKBR(data_file);
- KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
-
- % load GNV data GRACE A
- data_GNV = ['GNV1B_', date, '_A','_MDC-BAS1_01.asc'];
- data_file = fullfile(data_folder, data_GNV);
- GNVi = readGNV(data_file);
- RotMat_1 = Ri2e(GNVi(1,1));
- RotMat_dot_1 = Ri2e_dot(GNVi(1,1));
- r_eci_1 = RotMat_1'*GNVi(1,2:4)';
- v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
-
- state0(1:6,Mday) = [r_eci_1; v_eci_1];
-
- % load GNV data GRACE B
- data_GNV = ['GNV1B_', date, '_B','_MDC-BAS1_01.asc'];
- data_file = fullfile(data_folder, data_GNV);
- GNVi = readGNV(data_file);
- RotMat_1 = Ri2e(GNVi(1,1));
- RotMat_dot_1 = Ri2e_dot(GNVi(1,1));
- r_eci_1 = RotMat_1'*GNVi(1,2:4)';
- v_eci_1 = RotMat_1'*GNVi(1,5:7)' + RotMat_dot_1'*GNVi(1,2:4)';
-
- state0(7:12,Mday) = [r_eci_1; v_eci_1];
-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);
-eps_save = zeros(mKBR,nods,ItrNo_max+1);
-chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-
-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))
- toc
- end
-
- %% 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);
- iHx = zeros(mKBR,12,nods);
- iHc = zeros(mKBR,noc,nods);
- irhodot_deviation = zeros(mKBR,nods);
- for Mday = 1:nods
- FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
- File = fullfile(FolderName,'results/OutputArcParall', 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);
- iHx(:,:,Mday) = arcstru.Hx_save;
- iHc(:,:,Mday) = arcstru.Hc_save;
- irhodot_deviation(:,Mday) = arcstru.rhodot_deviation;
- 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);
- yhat = zeros(mKBR,nods);
- for Mday = 1:nods
- xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
- % estimate range rate residuals
- yhat(:,Mday)=iHx(:,:,Mday)*xhat(:,Mday)+iHc(:,:,Mday)*chat;
- eps = irhodot_deviation(:,Mday)-yhat(:,Mday);
- eps_save(:,Mday,ItrNo+1) = eps;
-% figure;plot(eps)
-% xlabel('Observation Number');
-% ylabel('residuals [m/s]');
-% 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;
-
- % 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);
-% 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);
-
- % 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);
-
- fprintf('done with iteration %d \n', ItrNo)
-
-end
-
-%%
-% Filename and path of file to save
-FileName= ['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');
-
-
diff --git a/gfr_parallel4BMDC1_inertiapod.m b/gfr_parallel4BMDC1_inertiapod.m
deleted file mode 100644
index 7cab35ee3b024e6452a7c8948aeb2007d0f6eb9d..0000000000000000000000000000000000000000
--- a/gfr_parallel4BMDC1_inertiapod.m
+++ /dev/null
@@ -1,309 +0,0 @@
-%% GRACE gravity field recovery (GFR)
-%
-% GFR_PARALLEL 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. With
-% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
-% iteration. Therefore you need to save all necessary variables, e.g. use the
-% file save_vars2continue_itr.m to save them in one .mat file. In the end some
-% example plots are given you may change or extend.
-%
-% Written by Florian W??ske, ZARM Uni Bremen, and
-
-% Neda Darbeheshti, AEI Hannover, 2018-07.
-
-%%
-format longg;
-
-%% variable inputs
-
-% add folder with all functions
-FolderName = pwd;
-addpath(genpath(FolderName));
-
-% folder to save intermediate results
-mkdir('results/OutputArcParall');
-mkdir('results/OutputNormalEqu');
-
-% Define the number of iterations for batch processor,
-ItrNo_max = 10;
-
-% order of spherical harmonic coefficients to be estimated
-lmaxcs = 70;
-
-% parameters of background gravity field
-lmaxf = 70;
-
-% number of days
-nods = 12;
-
-% folder where observation data are stored
-data_folder = fullfile(FolderName,'input_data','MDC_basic_challenge1');
-addpath(genpath(data_folder));
-
-% number of observations in each arc (one day, 17280 is the standard)
-mKBR = 17280;
-
-%% 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 and Earth radius of gravity field
-GM = 0.3986004415E+15;
-ae = 0.6378136300E+07;
-
-% reference for comparison
-ggm05s = readGFC('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 = 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(egm96_cs(1:lmaxf+1,1:lmaxf+1),0);
-field = field_sc(:)';
-[outC1,outS1,nm1] = cs2vec(egm96_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);
-
-%% load observations data and initial states for all days from files
-
-t0 = [2005 05 01 0 0 2];
-date0 = time2str(t0);
-
-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);
- date(6)=[]; % wrong date format in Nedas MDC mock data
- % load KBR data
- data_KBR = ['KBR1B_', date, '_X','_MDC-BAS1_01.asc'];
- data_file = fullfile(data_folder, data_KBR);
- KBRL1Bi = readKBR(data_file);
- KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);
-end
-% Load initial states in inertial frame from POD
-%GNVi=load('POD_iter_4_do_70_nods_30_2005-5-30_GA_BMDC');
-GNVi=load('POD_iter_4_do_70_nods_12_2005-5-12_A');
-state0(1:6,1:Mday)=GNVi.state_est(1:6,1:Mday);
-%GNVi=load('POD_iter_4_do_70_nods_30_2005-5-30_GB_BMDC');
-GNVi=load('POD_iter_4_do_70_nods_12_2005-5-12_B');
-state0(7:12,1:Mday)=GNVi.state_est(1:6,1:Mday);
- % 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);
-eps_save = zeros(mKBR,nods,ItrNo_max+1);
-chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-dv_save = zeros (ItrNo_max+1,lmaxcs-1);
-
-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))
- toc
- end
-
- %% 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);
- iHx = zeros(mKBR,12,nods);
- iHc = zeros(mKBR,noc,nods);
- irhodot_deviation = zeros(mKBR,nods);
- for Mday = 1:nods
- FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
- File = fullfile(FolderName,'results/OutputArcParall', 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);
- iHx(:,:,Mday) = arcstru.Hx_save;
- iHc(:,:,Mday) = arcstru.Hc_save;
- irhodot_deviation(:,Mday) = arcstru.rhodot_deviation;
- 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);
- yhat = zeros(mKBR,nods);
- for Mday = 1:nods
- xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
- % estimate range rate residuals
- yhat(:,Mday)=iHx(:,:,Mday)*xhat(:,Mday)+iHc(:,:,Mday)*chat;
- eps = irhodot_deviation(:,Mday)-yhat(:,Mday);
- eps_save(:,Mday,ItrNo+1) = eps;
-% figure;plot(eps)
-% xlabel('Observation Number');
-% ylabel('residuals [m/s]');
-% 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;
-
- % 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);
-% 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);
-
- % 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);
-
- fprintf('done with iteration %d \n', ItrNo)
-
-end
-
-% plot chat from all iterations:
-figure;
-for i = 1:ItrNo_max+1
- str = ['iter ',num2str(i-1)];
- semilogy(2:lmaxcs,chat_dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
- hold on
-end
-fs = 12;
-set(gcf,'PaperPositionMode','auto')
-set(gca,'FontSize',fs);
-xlabel('Spherical harmonic degree n')
-ylabel('Square root of degree variances')
-title('chat from each iteration')
-grid on
-legend
-
-figure;
-dv_geoidn(ggm05s,lmaxcs);
-hold on
-% dv_geoidn(d_vec_i,lmaxcs);
-
-[outC1,outS1,nm1] = cs2vec(field_cs,false);
-final_field_vec = [nm1',outC1',outS1'];
-dv_geoidn(final_field_vec,lmaxcs);
-
-d_vec = ggm05s(1:n_back,3:4)-final_field_vec(:,3:4);
-d_vec = [final_field_vec(:,1:2) d_vec];
-dv_geoidn(d_vec,lmaxcs);
-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
-dv_geoidn(field_vec0,lmaxcs);
-dv_geoidn(err_vec,lmaxcs);
-legend('ggm05s','ggm05s + error field', 'error field')
-for i = 1:ItrNo_max+1
- str = ['error iter ',num2str(i-1)];
- semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
-end
-
-% plot range rate residuals in time and frequency for last iteration
-figure;
-plot(eps)
-xlabel('Observation Number');
-ylabel('residuals [m/s]');
-title(['Range-rate residuals, Itr',num2str(ItrNo,'%.2d'),',Day',num2str(Mday,'%.2d')])
-ASD_pw(eps,.2);
-%ASD(eps,.2);
-
diff --git a/luis/MyBatchHLSST20h.sh b/luis/MyBatchHLSST20h.sh
deleted file mode 100644
index bb464b4191101c2f00c3c6a4b0800b0db4d3f8cc..0000000000000000000000000000000000000000
--- a/luis/MyBatchHLSST20h.sh
+++ /dev/null
@@ -1,24 +0,0 @@
-#!/bin/bash -login
-#PBS -l nodes=1:ppn=6
-#PBS -l walltime=12:00:00
-#PBS -l mem=16GB
-#PBS -M weigelt@ife.uni-hannover.de
-#PBS -m abe
-#PBS -o /home/nhglmatt/Matlab/HLSST/${PBS_JOBNAME}_${PBS_JOBID}.txt
-#PBS -k oe
-#PBS -j oe
-
-# Call the calculation
-for midx in `seq 1 1 12`
-do
- echo "Starting processing $name for year $year and month $midx ..."
- # load Matlab module
- module load MATLAB/2017a
- # start Matlab
- matlab –nosplash –nodesktop –nodisplay -r "cd('/home/nhglmatt/Matlab/HLSST/');myadd2path;hlSST('$name',[$year $midx $year $midx],true,'$TMPDIR')"
- wait $!
- # release Matlab module
- module del MATLAB/2017a
-done
-
-
diff --git a/luis/hello.m b/luis/hello.m
deleted file mode 100644
index dac544c58dcfbb73329cf28345c56a70ef422b79..0000000000000000000000000000000000000000
--- a/luis/hello.m
+++ /dev/null
@@ -1,4 +0,0 @@
-% Example MATLAB script for Hello World
-disp ’Hello World ’
-exit
-% end of example file
diff --git a/luis/luis_command.txt b/luis/luis_command.txt
deleted file mode 100644
index a1ca9ae38d74126a652ab17d4a04aa9d6c0e8b89..0000000000000000000000000000000000000000
--- a/luis/luis_command.txt
+++ /dev/null
@@ -1,14 +0,0 @@
-login
-
-ssh nhbfneda@login.cluster.uni-hannover.de
-
-scp /home/neda/Documents/code/gracetools-master/gfr_parallel.m nhbfneda@login.cluster.uni-hannover.de:/home/nhbfneda/
-
--r for directory
-
-for mac
-ssh aeissh
-ssh luis
-
-module load MATLAB/2017a
-matlab
diff --git a/luis/matlab-job-serial.sh b/luis/matlab-job-serial.sh
deleted file mode 100644
index 87f92836e4fd6e09af4fdb90853a02af574abc04..0000000000000000000000000000000000000000
--- a/luis/matlab-job-serial.sh
+++ /dev/null
@@ -1,18 +0,0 @@
-#!/bin/bash
-#PBS -N matlab_serial
-#PBS -M my@email.address
-#PBS -m bae
-#PBS -j oe
-#PBS -l nodes =1:ppn=1
-#PBS -l walltime =00:10:00
-#PBS -l mem=4gb
-# Compute node the job ran on
-echo "Job ran on:" $HOSTNAME
-# Load modules
-module load MATLAB /2017a
-# Change to work dir:
-cd $PBS_O_WORKDIR
-# Log file name
-LOGFILE=$(echo $PBS_JOBID | cut -d"." -f1).log
-# The program to run
-matlab -nodesktop -nosplash < hello.m > $LOGFILE 2>&1