% animate_2D_vib.m

%

% This MATLAB m-file "animates" a vibrational mode

% for the 2-D lattice.

%

% K. Beers. MIT ChE. 10/4/2002

% load in results of calculation

function iflag_main = animate_2D_vib();

iflag_main = 0;

% read in results of normal mode analysis

load lattice_2D_vib.mat;

% display normal modes and their frequencies

disp(' ');

disp('Normal mode frequencies (rad/sec) : ');

disp('-----------------------------------');

for j=1:length(freq)

    disp([int2str(j), '  ', num2str(freq(j))]);

end

disp(' ');

% select normal mode to animate

k_show = input('Enter number of mode to visualize : ');

% Create vector of theta values for a single

% oscillation.

num_frames = 50;

theta_vect = linspace(0,2*pi,num_frames);

[sin_max,iframe_max] = max(sin(theta_vect));

if(length(iframe_max) > 1)

    iframe_max = iframe_max(1);

end

% Set amplitude of oscillation so that the

% maximum displacement is equal to 1/3 times

% the equilibirum bond length.

displace_max = max(abs(V_eig(:,k_show)));

a = l_b/3/displace_max;

% Create movie by plotting 2-D lattice positions

title_phrase = ['2-D lattice, N = ', int2str(N), ...

        ', mode # = ', int2str(k_show), ...

        ', \omega = ', num2str(freq(k_show))];

movie_fig = figure;

for iframe = 1:num_frames

    hold off;

    % compute departure vector

    delta = a*sin(theta_vect(iframe))*V_eig(:,k_show);

    q = q_min + delta;

    % plot positions of atoms

    for i=1:N

        for j=1:N

            [k,kx,ky] = get_master_index(i,j,N);

            plot(q(kx),q(ky),'o');

            hold on;

        end

    end

    axis([ (-l_b) (N+l_b) (-l_b) (N+l_b) ]);

    title(title_phrase);

    % get frame for movie

    Movie_of_mode(iframe) = getframe;

    % if maximum amplitude, make separate plot

    if(iframe == iframe_max)

        figure;

        for i=1:N

        for j=1:N

            [k,kx,ky] = get_master_index(i,j,N);

            plot(q(kx),q(ky),'o');

            hold on;

            plot(q_min(kx),q_min(ky),'x');

            x_line = [q_min(kx), q(kx)];

            y_line = [q_min(ky), q(ky)];

            plot(x_line,y_line);

        end

        end

        axis([ (-l_b) (N+l_b) (-l_b) (N+l_b) ]);

        title([title_phrase, ', max. amplitude']);

        figure(movie_fig);

    end

end

% Show several periods of oscillation, repeat

% again and again until stop.

iquit = 1;

times_repeat = 10;

FPS = num_frames/2;  % frames per second

while(iquit)

    movie(Movie_of_mode,times_repeat,FPS);

    iquit = input('Quit (0) or Show again (1) : ');

end

iflag_main = 1;

return;

% ===================================================

% get_master_index.m

%

% This function calculates the master index label of

% a site in the 2D lattice of NxN sites.

function [k,kx,ky] = get_master_index(i,j,N);

k = (i-1)*N+j;

kx = 2*k-1;

ky = 2*k;

return;