% Code for Problem 3 of HW 10.

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

n = 20;            % # 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

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

part = input('Which part of Problem 3 are you doing?  Enter 2 for (b) or 3 for (c) :  ')

if part == 2
    
    xx = xx';       % one needs to use a column  xx  in this part
    disp(sprintf(' Vector  xx  and matrix  A  have been set up. \n Now follow instructions for part (b).'))
    
elseif part == 3

    epsilon = 0.0001;          % amplitude of the noise that will be added to the solution mu(x)
    
    w = ones(1,n);           % gold density that engenders the exact solution mu(x) below

    % Set up the unperturbed analytical form of the "measurement" mu(x),
    %  using  aux{1,2}_4mu  as auxiliary quantities for this setup:
	aux1_4mu = sqrt(1+xx.^2);
	aux2_4mu = sqrt(1+(xx-1).^2);     
	mu =(aux1_4mu + xx.*(aux2_4mu - aux1_4mu))./(aux1_4mu.*aux2_4mu);    

    noise = epsilon*(1-2*rand(1,n));
	mu_pert = mu + noise;      % the "measurement" perturbed by a small noise of maximum magnitude epsilon
	
	w_pert = inv(A)*mu_pert';  % gold density reconstructed from the slightly perturbed measurement mu_pert
    
    disp('  You will be asked to run this part 5 times')
    disp('  for 5 different seeds of the randon number generator')
    disp('  (which are chosen automatically by Matlab).')
    disp('  Enter one number at a time, from 1 to 5:')
    linecolor = input('Select a color for plots: 1=blue, 2=green, 3=magenta, 4=black, 5=cyan: ')
    if linecolor == 1
        line_color = 'b';
    elseif linecolor == 2
        line_color = 'g';
    elseif linecolor == 3
        line_color = 'm';
    elseif linecolor == 4
        line_color = 'k';
    elseif linecolor == 5
%         line_color = 'r';
%     elseif linecolor == 6
        line_color = 'c';
    end
	
    % Plot gold density w(x):
	figure(11);
	plot(xx, w_pert, line_color,  xx, w, 'r--',  'linewidth',2)
	legend('reconstructed as A^{-1}*\mu','exact')
	title(sprintf('Numerically reconstructed \n and exact gold densities w(x)'),'fontsize',14)
	
    % Plot "measurement" mu(x):
	figure(12);
	plot(xx, mu_pert,  xx, A*w','r--', 'linewidth',2)
	legend('perturbed','unperturbed')
	title('Perturbed and unperturbed "measurement" \mu(x)','fontsize',14)
         % The reason the unperturbed "measurement" 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 (unperturbed) measurement  mu  is also entered correctly.
         
    % Plot realizations of noise that were added to the "measurements" (before density was reconstructed):
    figure(13);
    plot(xx, noise, line_color, 'linewidth',2)
	title(sprintf('Realizations of noise \n added to the "measurement" mu(x)'),'fontsize',14)
         
 else
     
     disp('You entered something other than 2 or 3. Rerun the code.')
     
 end