diff --git a/src/FISTA_secondOrder.m b/src/FISTA_secondOrder.m
index 8501cc74bd972b90bb1f3a0de2ed6f05ec584c55..7162c415d26a4e0b1b7d3f6fdbbfbb0f4e8cf14d 100644
--- a/src/FISTA_secondOrder.m
+++ b/src/FISTA_secondOrder.m
@@ -1,176 +1,170 @@
-function [F, nr_of_iterations, Resnorm] = FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)
-%% L1SCHEME:
-% INPUT: farfield - far field matrix
-% kappa - wave number
-% Omega - part of [0,2pi] where the far field pattern is unkown
-% z - 2xm-vector containing the centers of the circles describing
-% the a priori information on the source locations
-% r - 1xm-vector containing the radii of the circles describing
-% the a priori information on the source locations
-% OUTPUT: f - cell array containing the m source components
-% EXAMPLE USAGE:
-
-nr_of_sources = size(z,2);
-
-
-% split and complete far field using iterated soft shrinkage
-landweberparam.kmax = kmax; % max nr of iterations of iterated soft shrinkage algorithm
-landweberparam.tol = 1e-3; % tolerance to be met in the Landweber iteration (stopping criterion)
-landweberparam.omega = 1/(nr_of_sources+1); % relaxation parameter in the Landweber itertation
-softshrinkparam.mu = 1e-3; % threshold parameter in the soft shrinkage operatore
-softshrinkparam.type = 'complex'; % type of soft-shrinkage to be used ('real'/'complex')
-
-softshrinkparam.weights = landweberparam.omega * softshrinkparam.mu;
-
-% fast iteraitve thresholding algorithm:
-
-[F, nr_of_iterations, Resnorm] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam);
-
-indOfInterest = sub2ind([nr_of_sources,nr_of_sources],1:nr_of_sources,1:nr_of_sources);
-
-F = {mysum(F, nr_of_sources^2), F{indOfInterest}};
-
-end
-
-function [F, nr_of_iterations, Resnorm] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam)
-%% FISTA:
-% INPUT: umeas - far field pattern
-% kappa - wave number
-% Omega -
-% z - 2xnr_of_sources vector containing presumed source positions
-% landweberparam.Nmax -- max nr of iterations
-% landweberparam.tol -- tolerance to be met in the Landweber iteration
-% landweberparam.omega -- relaxation parameter in the Landweber itertation
-% softshrinkparam.mu -- threshold parameter in the soft shrinkage operatore
-% softshrinkparam.type -- type of soft-shrinkage to be used ('real'/'complex')
-% OUTPUT: f - cell array containing split far field components
-%
-% EXAMPLE USAGE:
-
-
-% initialize some variables
-nxhat = sampling.nxhat;
-nd = sampling.nd;
-
-nr_of_sources = size(z,2);
-
-[X,Y] = meshgrid(1:nr_of_sources, 1:nr_of_sources);
-indices = [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
-
-% initialize some more variables
-X = cell(nr_of_sources^2,1); % columns of B correspond to the vectors b_j in the manuscript
-for iters=1:nr_of_sources^2
- X{iters} = zeros(nxhat,nd);
-end
-KstarRes = cell(nr_of_sources^2,1);
-iteri = 1; % iteration counter
-Resnorm = zeros(landweberparam.kmax, 1); % norm of the residual in each step of the iteration
-Resnorm(1) = landweberparam.tol+1; % just needed to get iteration started; will be removed later
-
-Y = X; Xalt = X;
-Yhelp = Y;
-talt = 1;
-
-% disp('FISTA Algorithm starts.')
-% tic
-
-% fast iterated soft shrinkage
-while (Resnorm(iteri) > landweberparam.tol) && (iteri < landweberparam.kmax)
-
- % apply operator K
- for iters = 1 : nr_of_sources^2
- incidentDir = indices(iters,1);
- detectorPos = indices(iters,2);
- Yhelp{iters} = fft2(translOp(ifft2(Y{iters})*(nxhat*nd),sampling,-z(:,incidentDir),-z(:,detectorPos), kappa))/(nxhat*nd);
- clear incidentDir detectorPos
- end
- KB = ifft2(mysum(Yhelp,nr_of_sources^2))*(nxhat*nd);
- KB = KB .* P_OmegaC; % apply P_OmegaC
-
- % evaluate residual
- Res = G - KB;
-
- % apply operator K^*
- Res = Res .* P_OmegaC; % apply P_OmegaC
- Res = fft2(Res)/(nxhat*nd);
- for iters = 1 : nr_of_sources^2
- incidentDir = indices(iters,1);
- detectorPos = indices(iters,2);
- KstarRes{iters} = fft2(translOp(ifft2(Res)*(nxhat*nd),sampling,z(:,incidentDir),z(:,detectorPos),kappa))/(nxhat*nd);
- clear incidentDir detectorPos
- end
-
- X = myplus(Y, myprod(landweberparam.omega,KstarRes,nr_of_sources^2),nr_of_sources^2);
-
- % apply soft shrinkage
- X = softshrink(X, softshrinkparam.weights, softshrinkparam.type, nr_of_sources^2);
-
- t = (1+sqrt(1+4*talt^2))/2;
-
- Y = myplus(X,myprod((talt-1)/t,(myplus(X,myprod(-1,Xalt,nr_of_sources^2),nr_of_sources^2)),nr_of_sources^2),nr_of_sources^2);
-
- Xalt = X;
- talt = t;
-
- Resnorm(iteri+1) = norm(Res); % HS norm of residual
- iteri = iteri+1; % update iteration counter
-
-end
-
-% toc
-% disp('FISTA Algorithm finished.')
-%
-% fprintf('\n\nNumber of iterations: %i\n', iteri);
-
-nr_of_iterations = iteri;
-
-%% shift back and compute inverse Fourier transform
-
-F = cell(nr_of_sources^2, 1); % cell array containing splitted far field components
-for iters = 1 : nr_of_sources^2
- Fhelp = ifft2(X{iters})*(nxhat*nd);
- incidentDir = indices(iters,1);
- detectorPos = indices(iters,2);
- F{iters} = translOp(Fhelp,sampling,-z(:,incidentDir),-z(:,detectorPos),kappa);
- clear incidentDir detectorPos
-end
-
-
-end
-
-
-
-function Bsum = mysum(B, nr_of_sources)
-
-Bsum = B{1};
-for iterk = 2:nr_of_sources
-
- Bsum = Bsum + B{iterk};
-
-end
-end
-
-
-function Z = myplus(X, Y, nr_of_sources)
-
-Z = cell(nr_of_sources, 1);
-for iterk = 1:nr_of_sources
-
- Z{iterk} = X{iterk} + Y{iterk};
-
-end
-end
-
-
-function Z = myprod(alpha, X, nr_of_sources)
-
-Z = cell(nr_of_sources, 1);
-for iterk = 1:nr_of_sources
-
- Z{iterk} = alpha*X{iterk};
-
-end
-end
\ No newline at end of file
+function [F, nr_of_iterations, Res] = FISTA_secondOrder(G, sampling, kappa, P_OmegaC, z, kmax)
+
+% Solves the l1xl1 minimization problem for the splitting and/ or completion problem in Born
+% approximation of order 2 numerically by using the FISTA algorithm.
+%
+% INPUT: G Observed far field matrix, nxhat*nd-array.
+% sampling Structure containing information about the discretization.
+% kappa Wave number, >0.
+% P_OmegaC Projection operator corresponding to observable part Omega^C, nxhat*nd-array
+% with only 0- and 1-entries, only 1-entries for splitting only.
+% z Positions of the individual scatterers, 2*nr_of_sources-array,
+% 2*1-array for completion only.
+% kmax Maximal number of FISTA iterations for stopping criterion.
+%
+% OUTPUT: F Structure containing numerical reconstructions for solutions of
+% l1xl1 minimization problem.
+% nr_of_iterations Number of performed FISTA iterations.
+% Res Residuum-norms for all FISTA iterations, vector of length nr_of_iterations.
+%
+% *********************************************************************************************************
+
+nr_of_sources = size(z,2);
+
+% split and complete far field using iterated soft shrinkage:
+landweberparam.kmax = kmax; % max nr of iterations of iterated soft shrinkage algorithm
+landweberparam.tol = 1e-3; % tolerance to be met in the Landweber iteration (stopping criterion)
+landweberparam.omega = 1/(nr_of_sources+1); % relaxation parameter in the Landweber itertation
+softshrinkparam.mu = 1e-3; % threshold parameter in the soft shrinkage operatore
+softshrinkparam.type = 'complex'; % type of soft-shrinkage to be used ('real'/'complex')
+
+softshrinkparam.weights = landweberparam.omega * softshrinkparam.mu;
+
+% fast iteraitve thresholding algorithm:
+
+[F, nr_of_iterations, Res] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam);
+
+indOfInterest = sub2ind([nr_of_sources,nr_of_sources],1:nr_of_sources,1:nr_of_sources);
+
+F = {mysum(F, nr_of_sources^2), F{indOfInterest}};
+
+end
+
+
+function [F, nr_of_iterations, Res] = fista_secondOrder(G, sampling, kappa, P_OmegaC, z, landweberparam, softshrinkparam)
+
+% initialize some variables:
+nxhat = sampling.nxhat;
+nd = sampling.nd;
+
+nr_of_sources = size(z,2);
+
+[X,Y] = meshgrid(1:nr_of_sources, 1:nr_of_sources);
+indices = [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
+
+% initialize some more variables:
+X = cell(nr_of_sources^2,1); % columns of B correspond to the vectors b_j in the manuscript
+for iters=1:nr_of_sources^2
+ X{iters} = zeros(nxhat,nd);
+end
+KstarRes = cell(nr_of_sources^2,1);
+iteri = 1; % iteration counter
+Res = zeros(landweberparam.kmax, 1); % norm of the residual in each step of the iteration
+Res(1) = landweberparam.tol+1; % just needed to get iteration started; will be removed later
+
+Y = X; Xalt = X;
+Yhelp = Y;
+talt = 1;
+
+% disp('FISTA Algorithm starts.')
+% tic
+
+% fast iterated soft shrinkage:
+while (Res(iteri) > landweberparam.tol) && (iteri < landweberparam.kmax)
+
+ % apply operator K:
+ for iters = 1 : nr_of_sources^2
+ incidentDir = indices(iters,1);
+ detectorPos = indices(iters,2);
+ Yhelp{iters} = fft2(translOp(ifft2(Y{iters})*(nxhat*nd),sampling,-z(:,incidentDir),-z(:,detectorPos), kappa))/(nxhat*nd);
+ clear incidentDir detectorPos
+ end
+ KB = ifft2(mysum(Yhelp,nr_of_sources^2))*(nxhat*nd);
+ KB = KB .* P_OmegaC; % apply P_OmegaC
+
+ % evaluate residual:
+ Res = G - KB;
+
+ % apply operator K^*:
+ Res = Res .* P_OmegaC; % apply P_OmegaC
+ Res = fft2(Res)/(nxhat*nd);
+ for iters = 1 : nr_of_sources^2
+ incidentDir = indices(iters,1);
+ detectorPos = indices(iters,2);
+ KstarRes{iters} = fft2(translOp(ifft2(Res)*(nxhat*nd),sampling,z(:,incidentDir),z(:,detectorPos),kappa))/(nxhat*nd);
+ clear incidentDir detectorPos
+ end
+
+ X = myplus(Y, myprod(landweberparam.omega,KstarRes,nr_of_sources^2),nr_of_sources^2);
+
+ % apply soft shrinkage:
+ X = softshrink(X, softshrinkparam.weights, softshrinkparam.type, nr_of_sources^2);
+
+ t = (1+sqrt(1+4*talt^2))/2;
+
+ Y = myplus(X,myprod((talt-1)/t,(myplus(X,myprod(-1,Xalt,nr_of_sources^2),nr_of_sources^2)),nr_of_sources^2),nr_of_sources^2);
+
+ Xalt = X;
+ talt = t;
+
+ Res(iteri+1) = norm(Res); % HS norm of residual
+ iteri = iteri+1; % update iteration counter
+
+end
+
+% toc
+% disp('FISTA Algorithm finished.')
+%
+% fprintf('\n\nNumber of iterations: %i\n', iteri);
+
+nr_of_iterations = iteri;
+
+% shift back and compute inverse Fourier transform:
+
+F = cell(nr_of_sources^2, 1); % cell array containing splitted far field components
+for iters = 1 : nr_of_sources^2
+ Fhelp = ifft2(X{iters})*(nxhat*nd);
+ incidentDir = indices(iters,1);
+ detectorPos = indices(iters,2);
+ F{iters} = translOp(Fhelp,sampling,-z(:,incidentDir),-z(:,detectorPos),kappa);
+ clear incidentDir detectorPos
+end
+
+
+end
+
+
+function Bsum = mysum(B, nr_of_sources)
+
+Bsum = B{1};
+for iterk = 2:nr_of_sources
+
+ Bsum = Bsum + B{iterk};
+
+end
+end
+
+
+function Z = myplus(X, Y, nr_of_sources)
+
+Z = cell(nr_of_sources, 1);
+for iterk = 1:nr_of_sources
+
+ Z{iterk} = X{iterk} + Y{iterk};
+
+end
+end
+
+
+function Z = myprod(alpha, X, nr_of_sources)
+
+Z = cell(nr_of_sources, 1);
+for iterk = 1:nr_of_sources
+
+ Z{iterk} = alpha*X{iterk};
+
+end
+end