diff --git a/src/CG_secondOrder.m b/src/CG_secondOrder.m
index 0c931c3634785132e50fcca8abc0f49eb0454a3f..d1540276bfc1b6bf7e3923480cae955ca8d7a740 100644
--- a/src/CG_secondOrder.m
+++ b/src/CG_secondOrder.m
@@ -1,155 +1,171 @@
-function [F, nr_of_iterations, Res] = CG_secondOrder(G, sampling, kappa, P_Omega, z, R, kmax, tol)
-
-% initialize some variables
-nd = sampling.nd; % number of incident waves
-nxhat = sampling.nxhat; % number of detectors on the sphere
-
-nr_of_sources = length(R);
-N = ceil(exp(1)/2*kappa*R); % compute cutoff parameters for individual sources
-
-[X,Y] = meshgrid(1:nr_of_sources, 1:nr_of_sources);
-indices = [[0,0]; [reshape(X,nr_of_sources^2,1),reshape(Y,nr_of_sources^2,1)]]; % (J^2+1)*2-array that indicates the finite dimensional subspace where the corresponding block should lie in.
-% (0,0) indicates $V_\Omega$; (j,l) indicates $V_{N_l,N_j}$, i.e.
-% cutoff-parameter N_l in row direction und N_j in column direction
-clear X Y
-
-% right hand side of the LGS that should be solved:
-B = cell(nr_of_sources^2+1, 1);
-B{1} = zeros(nxhat,nd);
-for iterk = 2:nr_of_sources^2+1
- incidentDir = indices(iterk,1);
- detectorPos = indices(iterk,2);
- B{iterk} = projOp(G,sampling,kappa,z(:,incidentDir),z(:,detectorPos),N(incidentDir), N(detectorPos));
- clear incidentDir detectorPos
-end
-
-% disp('CG Algorithm for quadratic approximations starts.')
-% tic
-
-% CG method:
-F = cell(nr_of_sources^2+1, 1); % current iterate
-for iterk = 1:nr_of_sources^2+1
-
- F{iterk} = zeros(nxhat,nd); % zero as initial iterate
-
-end
-
-% R = myminus(B, applyM(F, sampling, nr_of_sources, kappa, z, N), nr_of_sources);
-R = B; % Since F=0.
-D = R;
-
-Res = NaN*ones(kmax,1);
-for iteri = 1:kmax
-
- MD = applyM(D, P_Omega, sampling, nr_of_sources, kappa, z, N, indices);
- DMD = myscalprod(D, MD, nr_of_sources);
- normR2 = mynorm2(R, nr_of_sources);
- alpha = normR2/DMD;
-
- F = myminus(F, myprod(-alpha, D, nr_of_sources), nr_of_sources);
-
- Rneu = myminus(R, myprod(alpha, MD, nr_of_sources), nr_of_sources);
-
- normRneu2 = mynorm2(Rneu, nr_of_sources);
- beta = normRneu2/normR2;
- R = Rneu;
-
- D = myminus(Rneu, myprod(-beta, D, nr_of_sources), nr_of_sources);
-
- Res(iteri, 1) = normR2;
- if normR2 <= tol
- Res = Res(1:iteri);
- kmax = iteri;
- break
- end
-end
-
-% toc
-% disp('CG Algorithm for quadratic approximations finished.')
-
-nr_of_iterations = kmax;
-
-
-function Y = applyM(X, P_Omega, sampling, nr_of_sources, kappa, z, N, indices)
-
-Y = cell(nr_of_sources^2+1, 1);
-
-for iterk = 1:nr_of_sources^2+1 % loop for rows
-
- if iterk == 1
-
- Yhelp = X{1};
- for iterl = 2:nr_of_sources^2+1 % loop for columns
- incidentDir = indices(iterl,1);
- detectorPos = indices(iterl,2);
- Yhelp = Yhelp + P_Omega.*projOp(X{iterl},sampling,kappa,z(:,incidentDir),z(:,detectorPos),N(incidentDir),N(detectorPos));
- clear incidentDir detectorPos
- end
- Y{1} = Yhelp;
- clear Yhelp
- else
- incidentDirk = indices(iterk,1);
- detectorPosk = indices(iterk,2);
- Yhelp = projOp(P_Omega.*X{1},sampling,kappa,z(:,incidentDirk),z(:,detectorPosk),N(incidentDirk),N(detectorPosk));
- for iterl = 2:nr_of_sources^2+1
- if iterl == iterk
- Yhelp = Yhelp + X{iterl};
- else
- incidentDirl = indices(iterl,1);
- detectorPosl = indices(iterl,2);
- Bla = projOp(X{iterl},sampling,kappa,z(:,incidentDirl),z(:,detectorPosl),N(incidentDirl),N(detectorPosl));
- Yhelp = Yhelp + projOp(Bla, sampling,kappa,z(:,incidentDirk),z(:,detectorPosk),N(incidentDirk), N(detectorPosk));
- end
- end
- Y{iterk} = Yhelp;
-
- end
-
-end
-
-
-function t = mynorm2(X, nr_of_sources)
-
-t = 0;
-for iterk = 1:nr_of_sources^2+1
-
- t = t + norm(X{iterk}, 'fro')^2;
-
-end
-
-
-function t = myscalprod(X, Y, nr_of_sources)
-
-t = 0;
-for iterk = 1:nr_of_sources^2+1
-
- A = X{iterk};
- B = Y{iterk};
-
- t = t + trace(A'*B);
-
-end
-
-
-
-function Z = myminus(X, Y, nr_of_sources)
-
-Z = cell(nr_of_sources^2+1, 1);
-for iterk = 1:nr_of_sources^2+1
-
- Z{iterk} = X{iterk} - Y{iterk};
-
-end
-
-
-
-function Z = myprod(alpha, X, nr_of_sources)
-
-Z = cell(nr_of_sources^2+1, 1);
-for iterk = 1:nr_of_sources^2+1
-
- Z{iterk} = alpha*X{iterk};
-
-end
-
-
+function [F, nr_of_iterations, Res] = CG_secondOrder(G, sampling, kappa, P_Omega, z, R, kmax, tol)
+
+% Solves the least squares problem for the splitting and/ or completion problem
+% in Born approximation of order 2 numerically by using the cg method.
+%
+% INPUT: G Observed far field matrix, nxhat*nd-array.
+% sampling Structure containing information about the discretization.
+% kappa Wave number, >0.
+% P_Omega Projection operator corresponding to non-observable part Omega, nxhat*nd-array
+% with only 0 and 1 entries, zero-matrix for splitting only.
+% z Positions of the individual scatterers, 2*nr_of_sources-array,
+% 2*1-array for completion only.
+% R Sizes of the individual scatterers, i.e. radii of balls containing them,
+% vector of length nr_of_sources, scalar for completion only.
+% kmax Maximal number of CG iterations for stopping criterion.
+% tol Tolerance for residuum norm as stopping criterion.
+%
+% OUTPUT: F Structure containing numerical reconstructions for solutions of LS problem.
+% nr_of_iterations Number of performed CG iterations.
+% Res Residuum-norms for all CG iterations, vector of length nr_of_iterations.
+%
+% *********************************************************************************************************
+
+% initialize some variables
+nd = sampling.nd; % number of incident waves
+nxhat = sampling.nxhat; % number of observation directions
+
+nr_of_sources = length(R);
+N = ceil(exp(1)/2*kappa*R); % compute cutoff parameters for individual scatterers
+
+[X,Y] = meshgrid(1:nr_of_sources, 1:nr_of_sources);
+indices = [[0,0]; [reshape(X,nr_of_sources^2,1),reshape(Y,nr_of_sources^2,1)]]; % (J^2+1)*2-array that indicates the finite dimensional subspace where the corresponding block should lie in.
+% (0,0) indicates $V_\Omega$; (j,l) indicates $V_{N_l,N_j}$, i.e.
+% cutoff-parameter N_l in row direction und N_j in column direction
+clear X Y
+
+% right hand side of the LGS that should be solved:
+B = cell(nr_of_sources^2+1, 1);
+B{1} = zeros(nxhat,nd);
+for iterk = 2:nr_of_sources^2+1
+ incidentDir = indices(iterk,1);
+ detectorPos = indices(iterk,2);
+ B{iterk} = projOp(G,sampling,kappa,z(:,incidentDir),z(:,detectorPos),N(incidentDir), N(detectorPos));
+ clear incidentDir detectorPos
+end
+
+% disp('CG Algorithm for quadratic approximations starts.')
+% tic
+
+% CG method:
+F = cell(nr_of_sources^2+1, 1); % current iterate
+for iterk = 1:nr_of_sources^2+1
+
+ F{iterk} = zeros(nxhat,nd); % zero as initial iterate
+
+end
+
+R = B; % Since F=0.
+D = R;
+
+Res = NaN*ones(kmax,1);
+for iteri = 1:kmax
+
+ MD = applyM(D, P_Omega, sampling, nr_of_sources, kappa, z, N, indices);
+ DMD = myscalprod(D, MD, nr_of_sources);
+ normR2 = mynorm2(R, nr_of_sources);
+ alpha = normR2/DMD;
+
+ F = myminus(F, myprod(-alpha, D, nr_of_sources), nr_of_sources);
+
+ Rneu = myminus(R, myprod(alpha, MD, nr_of_sources), nr_of_sources);
+
+ normRneu2 = mynorm2(Rneu, nr_of_sources);
+ beta = normRneu2/normR2;
+ R = Rneu;
+
+ D = myminus(Rneu, myprod(-beta, D, nr_of_sources), nr_of_sources);
+
+ Res(iteri, 1) = normR2;
+ if normR2 <= tol
+ Res = Res(1:iteri);
+ kmax = iteri;
+ break
+ end
+end
+
+% toc
+% disp('CG Algorithm for quadratic approximations finished.')
+
+nr_of_iterations = kmax;
+
+
+function Y = applyM(X, P_Omega, sampling, nr_of_sources, kappa, z, N, indices)
+
+Y = cell(nr_of_sources^2+1, 1);
+
+for iterk = 1:nr_of_sources^2+1 % loop for rows
+
+ if iterk == 1
+
+ Yhelp = X{1};
+ for iterl = 2:nr_of_sources^2+1 % loop for columns
+ incidentDir = indices(iterl,1);
+ detectorPos = indices(iterl,2);
+ Yhelp = Yhelp + P_Omega.*projOp(X{iterl},sampling,kappa,z(:,incidentDir),z(:,detectorPos),N(incidentDir),N(detectorPos));
+ clear incidentDir detectorPos
+ end
+ Y{1} = Yhelp;
+ clear Yhelp
+ else
+ incidentDirk = indices(iterk,1);
+ detectorPosk = indices(iterk,2);
+ Yhelp = projOp(P_Omega.*X{1},sampling,kappa,z(:,incidentDirk),z(:,detectorPosk),N(incidentDirk),N(detectorPosk));
+ for iterl = 2:nr_of_sources^2+1
+ if iterl == iterk
+ Yhelp = Yhelp + X{iterl};
+ else
+ incidentDirl = indices(iterl,1);
+ detectorPosl = indices(iterl,2);
+ Bla = projOp(X{iterl},sampling,kappa,z(:,incidentDirl),z(:,detectorPosl),N(incidentDirl),N(detectorPosl));
+ Yhelp = Yhelp + projOp(Bla, sampling,kappa,z(:,incidentDirk),z(:,detectorPosk),N(incidentDirk), N(detectorPosk));
+ end
+ end
+ Y{iterk} = Yhelp;
+
+ end
+
+end
+
+
+function t = mynorm2(X, nr_of_sources)
+
+t = 0;
+for iterk = 1:nr_of_sources^2+1
+
+ t = t + norm(X{iterk}, 'fro')^2;
+
+end
+
+
+function t = myscalprod(X, Y, nr_of_sources)
+
+t = 0;
+for iterk = 1:nr_of_sources^2+1
+
+ A = X{iterk};
+ B = Y{iterk};
+
+ t = t + trace(A'*B);
+
+end
+
+
+function Z = myminus(X, Y, nr_of_sources)
+
+Z = cell(nr_of_sources^2+1, 1);
+for iterk = 1:nr_of_sources^2+1
+
+ Z{iterk} = X{iterk} - Y{iterk};
+
+end
+
+
+function Z = myprod(alpha, X, nr_of_sources)
+
+Z = cell(nr_of_sources^2+1, 1);
+for iterk = 1:nr_of_sources^2+1
+
+ Z{iterk} = alpha*X{iterk};
+
+end