%%%%%% Problem Set 2 %%%%%%%%%%%
% Problem 2.2
% V. Anant (modification of original code by B. Goldberg)
% further modified W.M. Kaminsky 9/24/2006 
%   - change units to nanometers throughout
%   - make prettier graph (annotation, xlim, higher res analytic
%   expression)
%   - fixed mislabeling of standard deviation as variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% See the notes on the 6.728 website for information 
% on taking a discrete fourier transfom

% preamble
clear
close all
format compact

disp('====================');
disp('Problem Set 2 output');
disp('====================');
% Define an x vector
x = [-20:0.02:20].';

% Define an L vector
L = [0.5 1 2].';

% Define a matrix of Psi_L(x)
Psi(:,1) = (pi*L(1)^2)^(-1/4).*exp(-x.^2/(2*L(1)^2));
Psi(:,2) = (pi*L(2)^2)^(-1/4).*exp(-x.^2/(2*L(2)^2));
Psi(:,3) = (pi*L(3)^2)^(-1/4).*exp(-x.^2/(2*L(3)^2));

%%%% Part b %%%%%
deltax = 0.02;
N = (20 - -20)/deltax + 1;
deltaq = 2*pi/(N*deltax);
qmax = pi*(N-1)/(deltax*N);
qmin = -qmax;

% Create q vector
q = [qmin:deltaq:qmax].';

% Create phi matrix
phi = exp(i*(x*q.'));

% Calculate the fourier transform A(q)
A = phi'*Psi*deltax;

clf;
plot(q, abs(A(:,1)).^2,'x',q, abs(A(:,2)).^2,'o',q, abs(A(:,3)).^2,'+');

%%%%%%%%  Part c   %%%%%%%%%
disp('Part (c): A_L normalization');
disp('-----------------------');
fprintf('Sum should be 1 for all L''s\n\n');
Sum = [ sum(A(:,1)'*A(:,1)*deltaq/(2*pi));
        sum(A(:,2)'*A(:,2)*deltaq/(2*pi));
        sum(A(:,3)'*A(:,3)*deltaq/(2*pi)); ]
   
% Compare against the analytic expression for A(q)
qa = [-5:0.002:5].';  %% Essentially the entire nonzero part of the Fourier
                      %% transform happens between q = -5/nm to 5/nm
                      %% To get a nice Gaussian shape, we plot it out at
                      %% 10x higher resolution than in q
A_analytic(:,1) = (4*pi*L(1)^2)^(1/4)*exp((-qa.^2*L(1)^2)/2);
A_analytic(:,2) = (4*pi*L(2)^2)^(1/4)*exp((-qa.^2*L(2)^2)/2);
A_analytic(:,3) = (4*pi*L(3)^2)^(1/4)*exp((-qa.^2*L(3)^2)/2);

hold on;
plot(qa, abs(A_analytic).^2);
title('Fourier Transform of \psi(x) = [1/(\pi L^2)]^{1/4} e^{-x^2 /(2L^2)}');
xlim([-5 5])  %% Again, essentially the entire nonzero part of the Fourier transform happens between q = -5/nm to 5/nm
xlabel('q [nm^{-1}]');
ylabel('|A(q)|^2');
legend('L = 0.5nm', 'L = 1nm', 'L = 2nm');
annotation1 = annotation('textbox',...
  'Position',[0.6175 0.7014 0.245 0.06968],...
  'FitHeightToText','off',...
  'String',{'Solid lines denote analytical expression (4\piL^2)^{1/4} e^{-q^2 L^2 /2}','','x,o, and + signs denote numerical discrete Fourier transform.'});
 
hold off;

%%%% Part d %%%%%
disp('Part (d): variance of x,q');
disp('-----------------------');
mean_x = [ sum(Psi(:,1)'*(x.*Psi(:,1)));
           sum(Psi(:,2)'*(x.*Psi(:,2)));
           sum(Psi(:,3)'*(x.*Psi(:,3)))]*deltax;
meansq_x = [ sum(Psi(:,1)'*(x.^2.*Psi(:,1)));
             sum(Psi(:,2)'*(x.^2.*Psi(:,2)));
             sum(Psi(:,3)'*(x.^2.*Psi(:,3)))]*deltax;
std_x = sqrt(meansq_x - mean_x.^2);
mean_q = [ sum(A(:,1)'*(q.*A(:,1)));
           sum(A(:,2)'*(q.*A(:,2)));
           sum(A(:,3)'*(q.*A(:,3)))]*deltaq/(2*pi);
meansq_q = [ sum(A(:,1)'*(q.^2.*A(:,1)));
             sum(A(:,2)'*(q.^2.*A(:,2)));
             sum(A(:,3)'*(q.^2.*A(:,3)))]*deltaq/(2*pi);
std_q = sqrt(meansq_q - mean_q.^2);
uncertainty = abs(std_q.*std_x);

disp('Standard Deviation (Uncertainty) of x')
disp(std_x);
disp('Standard Deviation (Uncertainty) of q')
disp(std_q);
disp('Uncertainty Product')
disp(uncertainty)