% Calculates a 1-d Fourier transform of a given function, and plots the
% function and its transform. 


% Parameters

Nx = 2^6;    % Number of points
Lx = 5;      % Length of domain
x = linspace(0, Lx, Nx);   % vector of grid points


% A couple of functions to get you started

f1 = @(y) (exp(-(y-Lx/2).^2/2));
   % We use '.^' to apply the operator to each element of the vector y
   % individually. Using just '^' would mean matrix power. 
   
%f1 = @(y) ((y>=Lx/2-Lx/4).*(y<=Lx/2+Lx/4));
   % Similar comment, this time about '.*' . If we use '*', this is the 
   % dot product
   % This example uses logical operators.


% Calculate the value of the function at each grid point
f = f1(x);

% Calculate its Fourier transform
[k,ft] = mhfft(x,f);


% Plot the results!

figure(1)

plot(x,f);
xlabel('x')
ylabel('f(x)')
title('Real space')

figure(2)

plot(k,abs(ft))
xlabel('k');
ylabel('f_{hat}(k)')
title('Fourier transform');
axis tight    % makes the axes the same size as the data values

% Note that we could omit the number from 'figure', and then Matlab
% would use the next available number. What happens if we omit the command
% altogether?
%
% It is also possible to put both of these graphs on one plot - 
% see the help page for 'subplot'. 

% Try typing the commands 'k' or 'f' or 'ft' in the command window 
% to see what these variables look like. 


%%%===========%%%
%   EXERCISES   %
%%%===========%%%

% 1. Play around with different functions. See how the Fourier transform
%    compares to the function itself. 
%    (If you don't like this implementation of the Fourier transform, 
%     you can write another one which has 0 as the middle value of x.)
% 2. Write down the mathematical formulation of the Fourier transform 
%    we are using. What is the relationship between Nx, Lx, dx 
%    (= x(2)-x(1)), Lk, dk (= k(2) - k(1)) ? 
% 3. Write a code to start with the Fourier transform of a function, and 
%    calculate the function itself. (Hint: you might want to start with 
%    x, compute k from this using k_of_x, and then calculate the transform 
%    on points in k.) Most of the time you will get a
%    complex-valued function. What conditions on ft(k) ensure that the
%    inverse transform is real-valued? 



