function [t,y,y_p,a,b]=part_sum(fun, nth_sum,interval)
% PART_SUM computes the nth partial sum of a cosine/sine 
% Fourier series for 'fun', which is a function or expression 
% given as a string. 'nth_sum' specifies the partial sum desired, 
% and 'interval' specifies the interval on which 'fun' is defined. 
% If no 'interval' is given, the default is [0,2*pi].
%
% Francis J. Narcowich--9/28/98 (Revised 2/4/03)
n_fft=4096;
n=nth_sum;
if nargin==3,
  low=interval(1);
  high=interval(2);
else 
  low=0;
  high=2*pi;
end
t=linspace(low,high,n_fft+1);
fun=fcnchk(fun); %Yields an inline function for expressions.
y=feval(fun,t); 
f_jump=(y(1)-y(end))/2; %Jump at ends.
hat_vals=fft(y(1:end-1))-f_jump; %Correction to fft due to jump.
if low~=0 | high-low~=2*pi,
   period=high-low;
   phase=2*pi*low/(high-low);
   hat_vals=exp(-phase*i*(0:(n_fft-1))).*hat_vals;
else
   period=2*pi;
end
a=[hat_vals(1)/n_fft,2*real(hat_vals(2:end))/n_fft];
b=-2*imag(hat_vals(2:end))/n_fft;
indx=1:n;
C=[ones(1,length(t));cos(indx'*2*pi*t/period)];
S=sin(indx'*2*pi*t/period);
y_p=a(1:(n+1))*C + b(1:n)*S;