gfr for bmdc1 was added

gfr_parallel4BMDC1.m

+%% 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);
+