% fourierdemo.m
%
% Program to illustrate Fourier Series decomposition of a vibrating string.
% Can also show oscillations (normal mode decomposition) as a function of
% time.
% 
% Raghuveer Parthasarathy
% Nov. 2007
%
% Options (functions)
% 1.  Triangle
%    y(x) = kx,       x <  L/2
%           kL - kx,  x >= L/2
%    Fourier:  sine terms Bn = 4 *sin(j*pi/2)/ (n pi)^2 (odd n)
%                         Bn = 0 (even n)
% 2.  Weird rectangle function
%    y(x) = 0  for x < L/4
%         = 1  for L/4 <= x < 3L/5
%         = -1 for 3L/5 <= x < 7L/8
%         = 0 for 7L/8 <= x
% 3.  Gaussian
%    y(x) = exp(-(x-3L/4).^2 / 2 / 0.05 / 0.05);
% 
% Raghuveer Parthasarathy, 5 Nov. 2006


clear all
close all


% ----------------
disp('* Fourierdemo.m -- RP Nov. 2007')
disp(' ');
disp('Functions: ')
disp('   (1) Triangle');
disp('   (2) Weird rectangle function');
disp('   (3) Gaussian');
opt = input('Enter option:  ');
disp(' ');

N = input('Number of terms (max. n value) of Fourier Series:  '); 
    % number of terms in series

oscopt = logical(input('Display oscillations?  (1==yes):  '));

% ------------------------------------------



x = 0.0:0.005:1.0;  % x, in units of L
Nx = length(x);

% Define the functions
switch opt
    case 1
        %triangle
        y = x.*(x <= 0.5) + (1.0-x).*(x > 0.5);  % y, in units of k
    case 2
        % weird rectangle
        y = zeros(1, Nx);
        for j=1:Nx,
            if (x(j) < 0.25),
                y(j) = 0.0;
            elseif (x(j) < 0.6),
                y(j) = 1.0;
            elseif (x(j) < 0.875),
                y(j) = -1.0;
            else
                y(j) = 0.0;
            end
        end
    case 3
        % Gaussian
        y = exp(-(x-3/4).^2 / 2 / 0.05 / 0.05);
    otherwise
        disp('Error!  Invalid option for funtion.  Hit Control-C');
        pause
end


% Fourier components  (cosine amplitudes are zero)
B = zeros(1,N);
switch opt
    case 1
        % triangle
        % Analytically determine Fourier components
        for j=1:N,
            B(j) = 4 *sin(j*pi/2)/ j / j/ pi / pi;
            ycomp(j,:) = B(j)*sin(j*pi*x);
        end
    case {2 , 3}
        % weird rectangle, or Gaussian
        % computationally determine Fourier components -- I'm being lazy!
        for j=1:N
            B(j) = 2 * sum(y.*sin(j*pi*x)*mean(diff(x)));
            ycomp(j,:) = B(j)*sin(j*pi*x);
        end
end


% partial sum after N terms
Ntermy = zeros(size(y));
for j=1:N,
    Ntermy = Ntermy + ycomp(j,:);
end


% plot Fourier sum after N terms
figure(1)
plot(x, y, 'k-', 'LineWidth', 2.0)
xlabel('x / L', 'Fontweight', 'bold');
ylabel('y / k', 'Fontweight', 'bold');
title('y(x), and N term Fourier Sum', 'Fontweight', 'bold');
hold on
plot(x, Ntermy, 'b-', 'LineWidth', 2.0)
legend('y(x)', 'N term Fourier series', 'Location', 'Northwest');
legend('boxoff')
axis auto 

% Plot each component
figure(2)
plot(x, y, 'k-', 'LineWidth', 2.0)
xlabel('x / L', 'Fontweight', 'bold');
ylabel('y / k', 'Fontweight', 'bold');
title('y(x), and each Fourier component', 'Fontweight', 'bold');
hold on
for j=1:N,
    plot(x, ycomp(j,:), 'color', [0 1 0]*j/N/4 + [0 0.75 0], 'LineWidth', 2.0)
end
axis auto
a = axis;
axis([a(1) a(2) (a(3)-0.1) (a(4)+0.1)])


% Plot Fourier amplitudes
figure(3)
plot(B, 'ko-', 'LineWidth', 2.0)
xlabel('n', 'Fontweight', 'bold');
ylabel('B_n', 'Fontweight', 'bold');
title('Fourier amplitudes B_n', 'Fontweight', 'bold');
grid on





if oscopt==1,
    % Illustrate oscillations
    disp('  ');
    disp('Illustration of oscillations.');
    disp('Press enter to continue.');
    pause
    disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.');
    figure;
    pause

    w0 = 2*pi*0.2;  % 0.2 Hz "natural frequency"
    % Frequencies of the normal modes
    omega = zeros(1,N);
    for j=1:N,
        omega(j) = j*w0;
    end

    % calculate components and sum at each time interval
    k = 1;
    dt = 0.1;
    for t=0:dt:6*pi/omega(1),
        yfourier(k,:) = zeros(size(y));
        cla
        for j=1:N
            ycomptime(k,j,:) = ycomp(j,:)*cos(omega(j)*t);
            yfourier(k,:) = yfourier(k,:) + ycomp(j,:)*cos(omega(j)*t);
        end
        k = k+1;
    end
    Nt = k-1;

    % Illustration of oscillations
    axis([-0.05 1.05  min(min(yfourier)) max(max(yfourier))]);
    hold on
    grid on
    title('Oscillations of non-normal modes.', 'Fontweight', 'bold');
    xlabel('x', 'Fontweight', 'bold');
    ylabel('y (Fourier sum)', 'Fontweight', 'bold');
    dt = 0.05;
    for k = 1:Nt,
        % clear previous plot
        cla
        for j=1:N
            plot(x, squeeze(ycomptime(k,j,:)), 'color', ...
                [1 1 1]*0.7*(1.0 - (j/N)), 'Linewidth', 2.0);
        end
        plot(x, yfourier(k,:), 'r', 'linewidth', 2.0);
        if (k==1)
            disp('Press Enter to begin...');
            pause
        end
        pause(0.1);
    end

end