%   This short program calculates the fourier transform of a function f(x)
%   according to the formula

%        fhat(s) = integral f(x) exp(-i*x*s) dx.

%   User must input the value a. f is then sampled over the interval
%   -a < x < a at N points. The transform is
%   computed using the fast Fourier transform of MATLAB. The result
%   (after rescaling) is computed on the interval -s0 < s < s0 where
%   N and a should be chosen such that s0 < N*pi/(2a). The program plots
%   the real and imaginary part of fhat.

%   The function f must be  provided in the mfile f.m.



a = input('Enter the value of a (sampling range)    ')

s0 = input('Enter the value of s0 (viewing range)    ')

N  = input('Enter the value of N (number of sampling points)   ')


delx = 2*a/N;
xsample = -a : delx :a-delx;
xsample(N/2+1) = eps; 

dels = pi/a;

s = -(N/2)*dels:dels:((N-1)/2)*dels;
u = f(xsample);
d = fft(u)/N;
dd = fftshift(d);

fhat  = 2*a*exp(i*a*s).*dd;

M = max(abs(fhat))+.5;
subplot(2,1,1)
 plot(s,real(fhat))
 xlabel('s axis')
 title(' real part of fhat ')
 axis([-s0 s0 -M M])
subplot(2,1,2)
 plot(s,imag(fhat))
 xlabel('s axis')
 title(' imaginary part of fhat ')
 axis([-s0 s0 -M M])