% This code is a modified version of the code for Problems 5 and related ones
%  of Sec. 5.3 of the book by C.W. Groetsch, "Inverse problems".

% The length of the integration interval is scaled to equal 1 (one).

n = input('Enter the number of data points:  n = ')
                 % # of discretization subintervals in the entire interval
h=1/n;           % discretization step size
xx=linspace(h/2,1-h/2,n);       % vector of discretization points,
                                %  located at the midpoints of the subintervals
                          
w=ones(1,n);                    % gold density that engenders the exact solution mu(x) below

% The analytical solution  mu = A*w' : 
aux1_4mu=sqrt(1+xx.^2);
aux2_4mu=sqrt(1+(xx-1).^2);     %  auxiliary expressions used in the final expression of mu
mu=(aux1_4mu + xx.*(aux2_4mu - aux1_4mu))./(aux1_4mu.*aux2_4mu);    

% Matrix A:
A=h./( (xx'*ones(1,n) - ones(n,1)*xx).^2 + 1 ).^(3/2);   % matrix A of the discretized integral equation


disp('Hit  "Enter"  to plot the measured disturbance  mu(x).')

pause

figure(101);
plot(xx,A*w','linewidth',2)
title('The measured "solution" \mu(x)','fontsize',14)
     % The reason the solution is computed as A*w' rather than simply
     %  as  mu(x)  is that I am simultaneously verifying that matrix A is set up correctly
     %  and the exact solution  mu  is also entered correctly.


w_computed = inv(A)*mu';           % gold density reconstructed from the solution  muep

disp('Hit  "Enter"  to plot the measured disturbance in gravitational force,  mu(x).')

pause

figure(102);
plot(xx,w_computed,'linewidth',2)
title('Numerically reconstructed gold density w(x)','fontsize',14)
