% get_Fourier_series.m

%

% [a0,a,b,omega,iflag] = ...

%    get_Fourier_series(t_val,f_val,...

%    num_truncate,isymmetry,iplot,fun_name);

%

% This MATLAB function returns the Fourier series

% of an input function.

%

% Input :

% -------

% t_val = vector of time values at which f(t) is given

% f_val = vector of f(t) at time values in t_val

% num_truncate = the order at which to truncate the

%                Fourier series

% isymmetry = declares symmetry of f(t)

%      (-1, odd), (+1, even), (0, not specified)

% iplot = if nonzero, plot f(t) and Fourier series

% fun_name = name of function that is fitted

%

% K. Beers. MIT ChE. 8/9/2002

function [a0,a,b,omega,iflag] = ...

    get_Fourier_series(t_val,f_val,...

    num_truncate,isymmetry,iplot,fun_name);

iflag = 0;

% First, we calculate the period 2P.

num_time = length(t_val);

t_start = t_val(1);

t_end = t_val(num_time);

P = 0.5*(t_end - t_start);

% We now calculate the angular frequencies of each

% mode.

omega = zeros(num_truncate,1);

for m=1:num_truncate

    omega(m) = m*pi/P;

end

% Next, we determine the coefficient a0.

a0 = (trapz(t_val,f_val))/P;

% Then, we get the coefficients multiplying the

% cosine functions.

a = zeros(num_truncate,1);

if(isymmetry >= 0)  % if not odd function

    for m=1:num_truncate

        cos_val = cos(m*pi.*t_val./P);

        integrand = cos_val.*f_val;

        a(m) = (trapz(t_val,integrand))/P;

    end

end

% Then, we calculate the coefficients that multiply

% the sine functions.

b = zeros(num_truncate,1);

if(isymmetry <= 0)  % if not even function

    for m=1:num_truncate

        sin_val = sin(m*pi.*t_val./P);

        integrand = sin_val.*f_val;

        b(m) = (trapz(t_val,integrand))/P;

    end

end

% If desired, a plot of the function and it's truncated

% Fourier series representation is made.

if(iplot)

    f_FS = (a0/2)*ones(size(t_val));

    for m=1:num_truncate

        f_FS = f_FS + a(m).*cos(m*pi.*t_val./P) + ...

            b(m).*sin(m*pi.*t_val./P);

    end

    figure;

    plot(t_val,f_val);

    hold on;

    plot(t_val,f_FS,'-.');

    xlabel('t');

    ylabel('f(t)');

    title([fun_name]);

    legend('f(t)', 'Fourier series');

end

iflag = 1;

return;