Update evaluateFarfieldSecondOrder.m

src/evaluateFarfieldSecondOrder.m

@@ -20,21 +20,17 @@ function [F] = evaluateFarfieldSecondOrder(k, sampling, objects, par, R, q)
 % ************************************************************************************
 
 kR = k * R;
-% kz = k * z;
 
 nd = sampling.nd; % number of illumination directions
-nxhat = sampling.nxhat; 
-
-xhat = sampling.xhat; % number of detector positions on the sphere
+nxhat = sampling.nxhat; % number of observation directions
+xhat = sampling.xhat;
 d = sampling.d;
 
 doub = 8; % number of doublings for discretization
 N = 2^doub; % number of discretizations in one dimension
 h = 2*kR/N; % discretization step width
 
-[X,Y] = meshgrid(-(N/2)*h:h:(N/2-1)*h,-(N/2)*h:h:(N/2-1)*h); % Discretization of supp(q), more precisely of [-R,R]^2
-% x = reshape(X,[N^2,1]) + kz(1)*ones(N^2,1);
-% y = reshape(Y,[N^2,1]) + kz(2)*ones(N^2,1);
+[X,Y] = meshgrid(-(N/2)*h:h:(N/2-1)*h,-(N/2)*h:h:(N/2-1)*h); % Discretization of [-R,R]^2
 x = reshape(X,[N^2,1]);
 y = reshape(Y,[N^2,1]);
 
@@ -44,9 +40,7 @@ Q1 = q(1) * reshape(q1,[N,N]);
 q2 = qBorn(objects{2}, par(2,:), x/k, y/k);
 Q2 = q(2) * reshape(q2,[N,N]);
 
-
-%% Kernel Khat
-
+% Kernel Khat:
 Khat = zeros(N,N);
 J0R1 = besselj(0,kR);
 J1R1 = besselj(1,kR);
@@ -69,25 +63,20 @@ for j1=1:N
         
     end
 end
-Khat = 2*kR * ifftshift(Khat);  % the factor 2R corrects an error in Vainikko's paper
-
-
-%% Calculate part of far field operator that corresponds to second order scattering first on q1 then on q2:
+Khat = 2*kR * ifftshift(Khat);
 
+% Calculate part of far field operator corresponding to second order scattering first on q1 then on q2:
 F = zeros(nxhat,nd);
 
 for iteri = 1:nd
     ui = uinc(x,y,d);
     Ui = reshape(ui(:,iteri),[N,N]);
 
-%     Vs1 = fftshift(ifft2(Khat.*fft2(ifftshift(Q1.*Ui)))); % Born approximation of us by scattering on q1
-    
     Vs1 = ifft2(ifftshift(fftshift(Khat).*fftshift(fft2(Q1.*Ui))));% fftshift(ifft2(Khat.*fft2(ifftshift(Q1.*Ui))));
     
     c = getc(N,h);
-%     Us1 = -h^2 * ToepPhi(c,Vs1);
     Us1 = h^2 * ToepPhi(c,Q1.*Ui);
-    % F as Born approximation of far field of Us1:
+    % F as first order Born approximation of far field of Us1:
     Vcol = reshape(Q2.*Us1,(N)^2,1);
     Vhelp = repmat(Vcol,1,nxhat);
     E = exp(-1i*[x,y]*[cos(xhat), sin(xhat)].');
@@ -97,7 +86,7 @@ for iteri = 1:nd
 
 end
 
-F = 2*pi/nd * F; % weights of trapecoidal rule
+F = 2*pi/nd * F; % adding weights of trapecoidal rule
     
 end
 
@@ -118,16 +107,6 @@ rot = par(5)*pi/180;
 R = [ cos(-rot) -sin(-rot) ; sin(-rot) cos(-rot) ];
 shift = repmat([par(3); par(4)],1,n);
 
-% switch object
-%     case 'circle'
-%         al = par(1);
-%         be = par(2);
-%         y_mod = R*([y1.';y2.']-shift); % 2*n array
-%         q_value = lambda*double(sqrt((y_mod(1,:)/al).^2+(y_mod(2,:)/be).^2) <= 1);
-%         % q_value = double(and(abs(y_mod(1,:)) <= al, abs(y_mod(2,:)) <= be));
-%         q_value = q_value.';
-% end
-
 [x,~,~,~] = kurve(sqrt(n), object, par);
 [x_phi,x_rad] = cart2pol(x(:,1)-par(3)*ones(sqrt(n),1),x(:,2)-par(4)*ones(sqrt(n),1));