% Ncoupledosc.m
%
% Program to illustrate normal modes of N coupled oscillators
% Will also illustrate "redundancy" of solutions for "modes" n>N
% Loaded string
%
% Raghuveer Parthasarathy
% October 28, 2006

clear all
close all

N = round(input('Enter the number of oscillators (4):  '));
% Will consider 2*N+1 modes, though of course "unique" modes are only 1:N

w0 = 2*pi*0.2;  % 0.2 Hz "natural frequency"
C = 1;  % same amplitude for all modes

% beads on string, separated by l=1; string goes from x = 0 to N+1

% Frequencies of the modes
omega = zeros(1,2*N+1);
for j=1:2*N+1,
    omega(j) = 2*w0*sin(j*pi / 2 / (N+1));
end

% Amplitude array; dim 1 = p, dim 2 = n
% Also include endpoints (zero amplitude), so point p is at index p+1
A = zeros(N+2,2*N+1);
for n = 1:2*N+1,  % modes
   for p=1:N+2, % particles
        A(p,n) = C*sin((p-1)*n*pi/(N+1));
    end
end

% Will also plot continuous sine function along x:
x = 0:0.05:N+1;
Nx = length(x);
y_cont = zeros(2*N+1,Nx);
for n=1:2*N+1,
    y_cont(n,:) = C*sin(x*n*pi/(N+1));
end

% ------------------------
% Options
disp('coupledosc.m -- R. Parthasarathy 28Oct06');
disp('** Options **');
disp('   (1) Illustration of normal modes at t=0.');
disp('   (2) Illustration of oscillations of normal modes.');
disp('   (3) Illustration of oscillations of normal modes for "modes" n>N.');
opt = input('Enter the option number: ');

switch opt
    case 1
        disp('Illustration of normal modes at t=0.');

        % Illustration of normal modes at t=0
        % Initial conditions: \delta=0 (zero initial velocity)
        figure;
        subplot(N+1,1,1);
        for j=1:N+1,
            subplot(N+1,1,j);
            plot(x, y_cont(j,:), 'g-', 'Linewidth', 2.0);
            hold on
            plot([0:N+1], A(:,j), 'bo', 'LineWidth', 2.0);
            plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black
            grid on
            axis([-0.1 N+1+0.1 -1.2 1.2]);
            if (j==1)
                title('normal modes at t=0', 'Fontweight', 'bold');
            end
            if (j==N+1)
                xlabel('x', 'Fontweight', 'bold');
            end
            ylabel('y', 'Fontweight', 'bold');
        end


    case 2

        disp('Illustration of oscillations of normal modes.');
        disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.');
        figure;
        pause

        % Illustration of oscillations
        subplot(N,1,1);
        dt = 0.2;
        for t=0:dt:4*pi/omega(1),
            for j=1:N,
                subplot(N,1,j);
                % clear previous plot
                cla
                plot(x, y_cont(j,:)*cos(omega(j)*t), 'color', [0.7 0.7 0.7], 'Linewidth', 2.0);
                hold on
                plot([0:N+1], A(:,j)*cos(omega(j)*t), 'bo-', 'LineWidth', 2.0);
                plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black
                grid on
                axis([-0.1 N+1+0.1 -1.2 1.2]);
                if (j==1)
                    title('Oscillations of normal modes.', 'Fontweight', 'bold');
                end
                if (j==N+1)
                    xlabel('x', 'Fontweight', 'bold');
                end
                ylabel('y', 'Fontweight', 'bold');
            end
            if (t<dt)
                disp('Press Enter to begin...');
                pause
            end
            pause(0.1);
        end

    case 3

        disp('Illustration of oscillations of normal modes for "modes" n>N.');
        disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.');
        figure;
        pause
        % Illustration of oscillations -- aliased mode number (n > N)
        subplot(N+1,2,1);
        dt = 0.2;
        for t=0:dt:4*pi/omega(1),
            %    for j=1:(2*N+1),
            for j=1:(2*N+1),
                ycont_t(j,:) = y_cont(j,:)*cos(omega(j)*t);  % continuous function
                ypart_t(j,:) = A(:,j)*cos(omega(j)*t);  % discrete particle positions
            end
            for j=1:(2*N+1),
                plotno = (j<(N+2))*(2*j-1) + (j>(N+1))*((2-2*N)*j/(N-1) + 4*(N+1));
                    % which subplot element
                subplot(N+1,2,plotno);
                % clear previous plot
                cla
                plot(x, ycont_t(j,:), 'color', [0.7 0.7 0.7], 'Linewidth', 2.0);
                hold on
                plot([0:N+1], ypart_t(j,:), 'ro-', 'LineWidth', 2.0);
                plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black
                grid on
                axis([-0.1 N+1+0.1 -1.2 1.2]);
                if (j==N+1)
                    xlabel('x', 'Fontweight', 'bold');
                end
                ylabel('y', 'Fontweight', 'bold');
                fs = sprintf('Mode n=%d', j);
                title(fs, 'Fontweight', 'bold');
            end
            if (t<dt)
                disp('Press Enter to begin...');
                pause
            end
                pause(0.04);
        end
end